% function [pd]=mlpmns()
%This program is a direct conversion of the corresponding Fortran program in
%S. Zhang & J. Jin "Computation of Special Functions" (Wiley, 1996).

%       ========================================================
%       Purpose: This program computes the associated Legendre
%                functions Pmn(x) and their derivatives Pmn'(x)
%                for a given order using subroutine LPMNS
%       Input :  x --- Argument of Pmn(x)
%                m --- Order of Pmn(x),  m = 0,1,2,...,n
%                n --- Degree of Pmn(x), n = 0,1,2,...,N
%       Output:  PM(n) --- Pmn(x)
%                PD(n) --- Pmn'(x)
%       Examples:
%                m = 1,  N = 5,  x = .5
%                n        Pmn(x)           Pmn'(x)
%               -------------------------------------
%                0    .00000000D+00    .00000000D+00
%                1    .86602540D+00   -.57735027D+00
%                2    .12990381D+01    .17320508D+01
%                3    .32475953D+00    .62786842D+01
%                4   -.13531647D+01    .57735027D+01
%                5   -.19282597D+01   -.43977853D+01

%                m = 2,  N = 6,  x = 2.5
%                n        Pmn(x)           Pmn'(x)
%               -------------------------------------
%                0    .00000000D+00    .00000000D+00
%                1    .00000000D+00    .00000000D+00
%                2    .15750000D+02    .15000000D+02
%                3    .19687500D+03    .26625000D+03
%                4    .16832813D+04    .29812500D+04
%                5    .12230859D+05    .26876719D+05
%                6    .81141416D+05    .21319512D+06
%       =======================================================

% pm=[];pd=[];
% 'please enter m, n, and x ',
%        READ(*,*)M,N,X
% m=1;
% n=5;
% x=.5;
% m,n,x,
% [m,n,x,pm,pd]=lpmns(m,n,x,pm,pd);
% ,
% '  n        pmn(x)           pmn''(x)    ',
% ' -------------------------------------',
% for  j=0:n;
% j,pm(j+1),pd(j+1),
% end;
      %format(1x,i3,2d17.8);
      %format(1x,'m =',i2,',  ','n =',i2,',  ','x =',f5.1);



function [m,n,x,pm,pd]=lpmns(m,n,x,pm,pd);

%       ========================================================
%       Purpose: Compute associated Legendre functions Pmn(x)
%                and Pmn'(x) for a given order
%       Input :  x --- Argument of Pmn(x)
%                m --- Order of Pmn(x),  m = 0,1,2,...,n
%                n --- Degree of Pmn(x), n = 0,1,2,...,N
%       Output:  PM(n) --- Pmn(x)
%                PD(n) --- Pmn'(x)
%       ========================================================




for  k=0:n;
pm(k+1)=0.0d0;
pd(k+1)=0.0d0;
end;
if (abs(x) == 1.0d0);
for  k=0:n;
if (m == 0);
pm(k+1)=1.0d0;
pd(k+1)=0.5d0.*k.*(k+1.0);
if (x < 0.0);
pm(k+1)=(-1).^k.*pm(k+1);
pd(k+1)=(-1).^(k+1).*pd(k+1);
end;
elseif (m == 1);
pd(k+1)=1.0d+300;
elseif (m == 2);
pd(k+1)=-0.25d0.*(k+2.0).*(k+1.0).*k.*(k-1.0);
if (x < 0.0) pd(k+1)=(-1).^(k+1).*pd(k+1); end;
end;
end;
return;
end;
x0=abs(1.0d0-x.*x);
pm0=1.0d0;
pmk=pm0;
for  k=1:m;
pmk=(2.0d0.*k-1.0d0).*sqrt(x0).*pm0;
pm0=pmk;
end;
pm1=(2.0d0.*m+1.0d0).*x.*pm0;
pm(m+1)=pmk;
pm(m+1+1)=pm1;
for  k=m+2:n;
pm2=((2.0d0.*k-1.0d0).*x.*pm1-(k+m-1.0d0).*pmk)./(k-m);
pm(k+1)=pm2;
pmk=pm1;
pm1=pm2;
end;
pd(0+1)=((1.0d0-m).*pm(1+1)-x.*pm(0+1))./(x.*x-1.0);
for  k=1:n;
pd(k+1)=(k.*x.*pm(k+1)-(k+m).*pm(k-1+1))./(x.*x-1.0d0);
end;


for  j=0:n;
s = j;
ss = pm(j+1);
sss = pd(j+1);
end;