function [PL,DPLdt] = plkd(L,z)

%   Return a matrix PL containing the associated Legendre functions
%   P_L^0:P_L^L of degree L at the points z
%   Optionally return a matrix DPLdt containing the
%   derivative with respect to theta where z == cos(theta).
%
% -- Input arguments --
% L      : Degree, Integer >= 0
% z      : Row vector of points in [-1, 1] where functions are to be evaluated
% -- Output arguments --
% PL     : Matrix of size (L+1)*prod(size(z)) giving the function values at z
% DPLdt  : Matrix of size (L+1)*prod(size(z)) giving derivatives at z
%        : If called with only one output argument the derivative will not be calculated
%
%   See http://mathworld.wolfram.com/LegendrePolynomial.html
%   References:
%   A&S: Abramowitz & Stegun, "Handbook of Mathematical Functions".
%   PL: Paul Leopardi, "Computing derivatives of associated Legendre functions using Matlab".
%   S: Szego, "Orthogonal Polynomials".
%
%   [PL,DPLdt] = plkd(L,z)

%   $Revision: 1.2  $ Paul Leopardi 2003-08-26
%   Clean up comments
%   $Revision: 1.1  $ Paul Leopardi 2003-08-22
%   Simplify expressions and remove the calls to ultra
%   $Revision: 1.0  $ Paul Leopardi 2003-08-20
%   Based on spharmrd
%   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 plkd
  return;
end

CalcDt = nargout > 1;
if size(z,2) == 1
    z = z';
end
k = 1:L;
lz = prod(size(z));
PL = legendre(L,z);
if CalcDt
    DPLdt = zeros(L+1,lz);
    if L > 0
        pole = abs(z) == 1;
        if ~isempty(z(pole))

            % The following special cases for the derivatives at -1 and at 1
            % are defined in terms of limits. 
            % The formula for order 1 is PL Theorem 4.
            % It can be derived from S 4.21.7, p63, S 4.1.1, p58 and S 4.1.3, p59.
            % For all other orders, the derivative is 0, as per PL Theorems 3 and 5.
            %
            DPLdt(2,pole) = -L*(L+1)/2*PL(1,pole);
        end
        if ~isempty(z(~pole))
            oz = ones(size(z));
            ok = ones(size(k))';
            sint = sqrt(1-z.*z);

            % The following general formula for the derivative is PL Theorem 1.
            % It can be derived from A&S 8.5.2, p334, assuming -1 < z < 1.
            %
            DPLdt(1+k,~pole) = (((k-L-1).*(k+L))'*oz(~pole)).*PL(k,~pole) - ...
                (k'*oz(~pole)).*PL(1+k,~pole).*(ok*(z(~pole)./sint(~pole)));

            % The following special case is PL Theorem 2.
            %
            DPLdt(1,~pole) = PL(2,~pole);
        end
    end
end
