function [Legan,DerLegan] = LegDerLeg(order,uarg)

% [Legan,DerLegan] = LegDerLeg(order,uarg)
%
% Prepares two vectors  
% Legendre and Derivative of Legendre
% corresponding to given argument and given degree
% in lexicographic ordering
arrsize=((order+1)*(order+2))/2;
Legan = zeros(arrsize,1);
DerLegan = zeros(arrsize,1);


 coefficient=1;sqr_of_2n_plus1=1.;  % numerical coefficient */    
 power_for_brackets=1.;  % variable [1-uarg**2]**(n/2) */    
 sqrt_brackets=sqrt(1-uarg*uarg);       
 bracckets=1-uarg*uarg;        
 diagonal=0;underdiagonal=0; datasize=((order+1)*(order+2))/2;        
 prev_in_column=0; now_in_column=0; next_in_column=0;line=0; column=0;    
              %% this are calculations for the initial (diagonal) P_nn and d(P_nn) values */    
              %% and also calculations of the below diagonal P_(n+1)n and d[P_(n+1)n] */    
                
for line=0:order     
   diagonal=(line*(line+3))/2;    
   Legan(1+diagonal)=coefficient*power_for_brackets;                      
   if  bracckets > 0.
       DerLegan(1+diagonal)=(-line)*coefficient*(uarg)*power_for_brackets/bracckets;
   else
       DerLegan(1+diagonal)=0;
   end
   underdiagonal=diagonal+line+1;               
   if underdiagonal < datasize
       sqr_of_2n_plus1=sqrt(2*line+1);
       Legan(1+underdiagonal)= uarg*sqr_of_2n_plus1*Legan(1+diagonal);              
       DerLegan(1+underdiagonal)= uarg*sqr_of_2n_plus1*DerLegan(1+diagonal)+sqr_of_2n_plus1*Legan(1+diagonal);                
   end                 
   if line == 0 
       coefficient=coefficient*sqrt(2);
   end
   coefficient=coefficient*(2*line+1)/sqrt((2*line+1)*(2*line+2));                      
   power_for_brackets=power_for_brackets*sqrt_brackets;
end
 
for column=0:(order-1)
  prev_in_column =(column*(column+3))/2;                       
  now_in_column=prev_in_column+column+1; 
  for line=(column+1):(order-1)
      next_in_column=now_in_column +line+1;                
      Legan(1+next_in_column)=(1/(sqrt((line+1)*(line+1)-column*column)))*((2*line+1)*uarg*Legan(1+now_in_column)-sqrt(line*line-column*column)*Legan(1+prev_in_column) );
      DerLegan(1+next_in_column)=(1/(sqrt((line+1)*(line+1)-column*column)))*((2*line+1)*uarg*DerLegan(1+now_in_column)+(2*line+1)*Legan(1+now_in_column)-sqrt(line*line-column*column)*DerLegan(1+prev_in_column) );
      prev_in_column =now_in_column;
      now_in_column=next_in_column;
  end
end

