function [schmidt] = schmit(n,m,theta)
%function [schmidt] = schmit(n,m,theta);
%
%  Returns Schmidt normalized associated Legendre polynomial.
%  Input parameters:
%    Argument of polynomial is theta, in degrees.  Degree and 
%    order of polynomial are n and m, respectively, 0 <= m <= n.
%    NOTE: theta may be a scalar or vector
%  Output:
%    schmidt is the same size as theta.
%    NaN is returned if the function exits with an error

%  Written in MATLAB by Jeff Freymueller
%  Modified from Press, W.H., Flannery, 
%  B.P., Teukolsky, S.A., and Vetterling, W.T., "Numerical 
%  Recipes, The Art of Scientific Computing,"  Cambridge Univ. 
%  Press, 1986. 

%  This function is vectorized by making use of the .* operator
%  in MATLAB. a.*b gives (a1*b1, a2*b2, ..., an*bn), it multiplies
%  the each element of a with the corresponding element of b.

schmidt = NaN*ones(size(theta));

if (m<0 | m>n)
   disp('SCHMIT:  Bad argument detected');
   disp('Must have m >= 0, m < n');
   return;
end 

d2rad = pi/180;
x=cos(theta*d2rad);

%
% Recursive calculation of pmm for degree and order m
% 
pmm = ones(size(x));
if (m > 0)
   somx2 = sqrt((1.0-x).*(1.0+x));
   fact = 1.0;
   for i = 1:m,
      pmm = -fact*pmm.*somx2;
      fact = fact + 2.0;
   end
end

%
% Calculate the value for degree n
%  Special cases are n = m, and n = m+1
%

if (n == m)
   schmidt = pmm;
else
   pmmp1 = (2*m+1)*x.*pmm;
   if (n == m+1)
      schmidt=pmmp1;
   else
      for nn = m+2:n
         pnn = ((2*nn-1)*x.*pmmp1 - (nn+m-1)*pmm)/(nn-m);
         pmm = pmmp1;
         pmmp1 = pnn;
      end
      schmidt = pnn;
   end
end

%
% Schmidt normalization
%

if (m > 0)
   xnorm = sqrt(2*factorial(n-m)/factorial(n+m));
   schmidt = xnorm*schmidt;    
end
         
return;
