function plm = legendre(l,amu,tflag)
%LEGENDRE Associated Legendre function.
%   P = LEGENDRE(N,X) computes the associated Legendre functions 
%   of degree N and order m = 0, 1, ..., N, evaluated for each element
%   of X.  N must be a scalar integer less than 256 and X must contain
%   real values between -1 <= X <= 1.  
%
%   If X is a vector, P is an (N+1)-by-L matrix, where L = length(X).
%   The P(m+1,i) entry corresponds to the associated Legendre function 
%   of degree N and order m evaluated at X(i). 
%
%   In general, the returned array has one more dimension than X.
%   Each element P(m+1,i,j,k,...) contains the associated Legendre
%   function of degree N and order m evaluated at X(i,j,k,...).
%
%   The associated Legendre functions are defined as:
% 
%       P(n,m;x) = (-1)^m * (1-x^2)^(m/2) * (d/dx)^m { P(n,x) },
%
%   where P(n,x) is the Legendre polynomial of degree n. Note that
%   the first row of P is the Legendre polynomial evaluated at X 
%   (the M == 0 case).
%
%   SP = LEGENDRE(N,X,'sch') computes the Schmidt semi-normalized
%   associated Legendre functions SP(n,m;x). These functions are 
%   related to the (un-normalized) associated Legendre functions 
%   P(n,m;x) by:
%               
%   SP(n,m;x) = P(n,x), m = 0
%             = (-1)^m * sqrt(2*(n-m)!/(n+m)!) * P(n,m;x), m > 0
%
%   Examples: 
%     1. legendre(2, 0.0:0.1:0.2) returns the matrix:
% 
%              |    x = 0           x = 0.1         x = 0.2
%        ------|---------------------------------------------
%        m = 0 |   -0.5000         -0.4850         -0.4400
%        m = 1 |         0         -0.2985         -0.5879
%        m = 2 |    3.0000          2.9700          2.8800  
%
%     2. X = rand(2,4,5); N = 2; 
%        P = legendre(N,X); 
%
%     so that size(P) is 3x2x4x5 and 
%     P(:,1,2,3) is the same as legendre(N,X(1,2,3)). 

%   Acknowledgment:
%
%   This program is based on a Fortran program by Robert L. Parker,
%   Scripps Institution of Oceanography, Institute for Geophysics and 
%   Planetary Physics, UCSD. February 1993.
%
%   MATLAB version: D. L. Chen, 5/24/93, 1/20/95.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.15 $  $Date: 1998/08/24 18:29:03 $
%
%   Reference:
%     [1] M. Abramowitz and I.A. Stegun, "Handbook of Mathematical
%         Functions", Dover Publications, 1965, Ch. 8.
%     [2] J. A. Jacobs, "Geomagnetism", Academic Press, 1987, Ch.4.
%
%   Note on Algorithm:
%
%   LEGENDRE uses a three-term backward recursion relationship in M.
%   This recursion is on a version of the Schmidt semi-normalized 
%   Associated Legendre functions SPc(n,m;x), which are complex 
%   spherical harmonics. These functions are related to the standard 
%   Abramowitz & Stegun functions P(n,m;x) by
%
%       P(n,m;x) = sqrt((n+m)!/(n-m)!) * SPc(n,m;x)
%   
%   They are related to the Schmidt form given previously by:
%
%       SP(n,m;x) = SPc(n,0;x), m = 0
%                 = (-1)^m * sqrt(2) * SPc(n,m;x), m > 0

if nargin < 3
    tflag = 'unnorm';
end;

% Check degree constraints
lmx = 256;   % initialize lmx
if (prod(size(l)) > 1) | (l >= lmx),
    error(sprintf('N must be a positive scalar < %d.',lmx));
elseif l <= 0 
    plm = ones(size(amu));
    return
end

% Convert amu to a row vector
siz = size(amu); 
amu = amu(:)';

% Initialize parameters
nmx=2*lmx;
% The following three are parameters from Parker's Fortran code
tol = 0.100000001686238353e-15;
alogt = -0.368413619995117188e+02;
tstart = 0.999999993922529029e-08;

rootn=sqrt(0:nmx);
s=sqrt(1.0 - amu.^2);
plm = zeros(l+1,length(amu));
tlm = zeros(l+3,length(amu));

% Calculate TWOCOT, separating out the amu = -1,+1 cases first
twocot = amu;       % initialize dimensions of TWOCOT

kk = find(amu == -1);
if ~isempty(kk)
    twocot(kk) = Inf;
end;

kk = find(amu == 1);
if ~isempty(kk)
    twocot(kk) = -Inf;
end;

kk = find(s~=0);
if ~isempty(kk)
    twocot(kk)  = -2.0*amu(kk)./s(kk);
end;

% Small sine: prevent underflow By starting at m=m1-1<l

% First, replace s_copy = 0 with s_copy = NaN to avoid log(0) error messages.
% This is okay since we want nonzero s anyway.
scopy = s;
kk = find(s == 0);
scopy(kk) = NaN;
ind = find(l*log(scopy) < alogt);

if ~isempty(ind)
% Approx solution of x*ln(x)=y 
    v = 9.2-alogt./(l*s(ind));
    w = 1./log(v);
    m1=1+l*s(ind).*v.*w.*(1.0058+ w.*(3.819 - w*12.173));
    m1=min(l+1, floor(m1));

% column-by-column recursion
    for k = 1: length(m1)
        mm1 = m1(k);
        col = ind(k);
        tlm(mm1:l+1,col)= zeros(size(mm1:l+1))';

        % Start recursion with proper sign
        tlm(mm1,col)  = sign(rem(mm1,2)-0.5)*tstart;
        if amu(col) < 0, tlm(mm1,col)=sign(rem(l+1,2)-0.5)*tstart; end;

        % Recur from m1 to m=0, carrying normalizing 
        sumsq=tol;
        for m= mm1-2:-1:0,
            tlm(m+1,col)=((m+1)*twocot(col)*tlm(m+2,col)- ...
            rootn(l+m+3)*rootn(l-m)* tlm(m+3,col))/ ...
            (rootn(l+m+2)*rootn(l-m+1));
            sumsq=tlm(m+1,col)^2 + sumsq;
        end;
        scale=1.0/sqrt(2*sumsq - tlm(1,col)^2);
        tlm(1:mm1+1,col)=scale*tlm(1:mm1+1,col);
    end     % FOR loop

end % small sine IF loop

% Zero out the columns that have been calculated above, saving S
tmps = s(ind);
s(ind) = zeros(size(ind));
nind = find(s);

% Regular recursion from m=l to 0
if ~isempty(nind)

    % Produce normalization constant
    c=1.0;
    twol=2.0*l;
    for x= 2.0: 2: twol
        c=c - c/x;
    end;
    
    tlm(l+1,nind)= sqrt(c)*(-s(nind)).^l;
    tlm(l,nind)  = tlm(l+1,nind).*twocot(nind).*l./sqrt(twol);

    % Recur downwards to m=0
    for m= l-2: -1: 0,
        tlm(m+1,nind)=((m+1)*twocot(nind).*tlm(m+2,nind) ...
            - rootn(l+m+3)*rootn(l-m).* ...
            tlm(m+3,nind))/(rootn(l+m+2)*rootn(l-m+1));
    end;

end % IF loop

% Calculate the Schmidt semi-normalized functions
s(ind) = tmps;  % restore s
plm(1:l,:)=tlm(1:l,:);
plm(l+1,:)=tlm(l+1,:);

% Polar argument   (mu=+-1)
s0 = find(s == 0);
if ~isempty(s0)
    plm(1,s0)=amu(s0).^l;
end

if strcmp(tflag,'sch')
    % Calculate the standard Schmidt semi-normalized functions from these
    % functions. For m = 1,...,l, multiply by (-1)^m*sqrt(2)
    row1 = plm(1,:);
    plm = sqrt(2)*plm;
    plm(1,:) = row1;    % restore first row
    const1 = 1;
    for r = 2: l+1;
        const1 = -const1;
        plm(r,:) = const1*plm(r,:);
    end;

else
    % Calculate the standard A&S functions (i.e., unnormlized case) from 
    % these functions by multiplying each row by: sqrt((l+m)!/(l-m)!)
    const1 = 1;
    for r = 2: l+1
        m = r-1;

        % Find const for row r
        const = exp(0.5 * (sum(log(2:l+m)) - sum(log(2:l-m)) ) );

        plm(r,:) = const*plm(r,:);
    end;

end;

% Pack into output array
if length(siz) == 2 & min(siz) == 1 
    plm = reshape(plm,l+1,max(siz)); % Row or column vector case
else
    plm = reshape(plm,[l+1 siz]);
end

