function [YL,DYLdt,DYLdp] = spharmrd(L,z,phi)

%   Return a matrix YL containing the real spherical harmonics
%   Y_{L,-L}:Y_{L,L} of degree L at the points (z, phi),
%   with (z,phi) representing points on the sphere S^2 in cylindrical
%   coordinates.
%   Optionally return matrices DYLdt and DYLdp containing
%   partial derivatives with respect to theta and phi
%   where z == cos(theta).
%
% -- Input arguments --
% L      : Degree, Integer >= 0
% z      : Row vector of altitudes of points where functions are to be evaluated
% phi    : Row vector of longitudes of points where functions are to be evaluated
% -- Output arguments --
% YL     : Matrix of size (2*L+1)*prod(size(z)) giving the function values at z
% DYLdt  : Matrix of size (2*L+1)*prod(size(z)) giving theta derivatives at z
% DYLdp  : Matrix of size (2*L+1)*prod(size(z)) giving phi derivatives at z
%        : If called with only one output argument the derivatives will not be calculated
%
%   Layout of the matrix YL is:
%   YL(1,:)   = Y_{L,-L}(z,phi) = scaled P_L^L(z)*sin(L*phi)
%   ...
%   YL(L,:)   = Y_{L,-1}(z,phi) = scaled P_L^1(z)*sin(phi)
%   YL(L+1,:) = Y_{L,0} (z,phi) = scaled P_L(z)
%   YL(L+2,:) = Y_{L,1} (z,phi) = scaled P_L^1(z)*cos(phi)
%   ...
%   YL(N,:)   = Y_{L,L} (z,phi) = scaled P_L^1(z)*cos(L*phi)
%
%   The matrices DYLdt and DYLdp have appropriate partial
%   derivatives in the corresponding positions.
%
%   See http://mathworld.wolfram.com/SphericalHarmonic.html
%   for definition.
%   See also http://mathworld.wolfram.com/LegendrePolynomial.html
%   References: 
%   A,A&R: Andrews, Askey & Roy, "Special Functions",
%   Cambridge, 1999, pp302-303, 456-457.
%   A&S: Abramowitz & Stegun, "Handbook of Mathematical Functions",
%   pp333-334,778
%
%   [YL,DYLdt,DYLdp] = spharmrd(L,z,phi)

%   $Revision: 2.2  $ Paul Leopardi 2003-08-22
%   Improve comments and tidy up.
%   $Revision: 2.1  $ Paul Leopardi 2003-08-21
%   Factor PL calculation back out into plkd
%   $Revision: 2.0  $ Paul Leopardi 2003-08-20
%   Thoroughly debug
%   $Revision: 1.5  $ Paul Leopardi 2002-08-02
%   Incorporate diff_PL
%   $Revision: 1.4  $ Paul Leopardi 2002-08-01
%   Optimized spharmrd based on spharmvr and diff_PL
%   $Revision: 1.3  $ Paul Leopardi 2002-07-11
%   Optimized spharmvr based on spharmvz
%   $Revision: 1.2  $ Paul Leopardi 2002-07-10
%   Optimized spharmpl based on spharmx
%   $Revision: 1.1  $ Paul Leopardi 2002-07-09
%   Harmonize with spharmxc
%   $Revision: 1.0  $ Paul Leopardi 2002-07-02
%   Adapted by Paul Leopardi 2002-07-02
%   for UNSW School of Mathematics
%   Based on SPHARM2:
%   Denise L. Chen  9-1-93
%   Copyright 1984-2000 The MathWorks, Inc.
%   $Revision: 5.10 $  $Date: 2000/06/01 03:46:27 $

if nargin == 0
  help spharmrd
  return;
end

CalcDt = nargout > 1;
CalcDp = nargout > 2;

% Ensure input consists of row vectors
%
if size(z,2) == 1
    z = z';
end
if size(phi,2) == 1
    phi = phi';
end

% Get associated Legendre functions and derivatives.
%
if CalcDt
    [PL,DPLdt] = plkd(L,z);
else
    PL = plkd(L,z);
end

k = 1:L;
N = 2*L+1;
lz = prod(size(z));

YL = zeros(N,lz);
scaling = sqrt(N/(4*pi));
YL(1+L,:)        = scaling * PL(1,:);
if CalcDt
    DYLdt = zeros(N,lz);
    DYLdt(1+L,:) = scaling * DPLdt(1,:);
end
if CalcDp
    DYLdp = zeros(N,lz);
end
if L > 0
    cp = (cumprod(L:2*L).*[1,1,cumprod(L-1:-1:1)])';
    scaling = diag(sqrt(N./(2*pi*cp(k+1))));
    scaledPL  = scaling * PL(1+k,:);
    kkp = k'*phi;

    % sin comes first, in reverse order
    %
    sk = L:-1:1;

    % cos comes after P_L, in forward order
    %
    ck = L+2:N;

    YL(sk,:)        = scaledPL  .* sin(kkp);
    YL(ck,:)        = scaledPL  .* cos(kkp);
    
    % If required, use Leibniz rul to calculate derivatives
    %
    if CalcDt
        scaledDPLdt = scaling * DPLdt(1+k,:);
        DYLdt(sk,:) = scaledDPLdt .* sin(kkp);
        DYLdt(ck,:) = scaledDPLdt .* cos(kkp);
    end
    if CalcDp
        DYLdp(sk,:) = scaledPL  .* ( diag(k) * cos(kkp));
        DYLdp(ck,:) = scaledPL  .* (-diag(k) * sin(kkp));
    end
end
