Homec4science
Diffusion COSOlver (new_coulomb) Browse

cosolver/matlab/legendrePlm.m
41e7a23b16aenew_coulomb

Last Change
Raw File

legendrePlm.m
View Options

function [p, err]=legendrePlm(l,m,x)
% Calculate the Associated Legendre polynomials. Speed is obtained by
% direct calculation of polynomial coefficients rather than recursion.
% Polynomial coefficients can increase in magnitude very quickly with
% polynomial degree, leading to decreased accuracy (estimated by err return
% parameter). If you need higher degrees for the polynomials, use
% recursion-based algorithms.
% If called with two return values, the second is a scalar, estimating the
% error, based on worst case scenario and the significant digits. Error is
% based on the largest polynomial coefficient and machine error, "eps"
% maxcf is used to store the largest polynomial coefficient If there are
% only two input arguments, the second argument is x, and m defaults to 0;
% then legendreP is the Legendre polynomial of degree m
if(nargin==2);
x=m;
clear m;
m=0;
end;
%Some basic error checking of parameters
% Definition states that if |m|>l, the polynomial is 0; so set to 0 rather
% than return an error. Some algorithms depend on this behaviour. Note that
% the original definition requires 0<=m<=l
if(abs(m)>abs(l)), p=zeros(size(x)); err=0; return; end;
% Definition for l<0
if(l<0), l=-l-1; end;
% For m<0, polynomials are proportional to those with m>0, cfnm is the
% proportionality coefficient
cfnm=1;
if(m<0),
m=-m;
cfnm=(-1)^m*factorial(l-m)/factorial(l+m);
end;
% Calculate coef of maximum degree in x from the explicit analytical
% formula
cl=(-1)^m*cfnm*factorial(2*l)/((2^l)*factorial(l)*factorial(l-m));
maxcf=abs(cl);
%fprintf('Coef is %.16f\n',cl);
px=l-m;
% Power of x changes from one term to the next by 2. Also needed for
% sqrt(1-x^2).
x2=x.*x;
% Calculate efficiently P_l^m (x)/sqrt(1-x^2)^(m/2) - that is, only the
% polynomial part. At least one coefficient is guaranteed to exist - there
% is no null Legendre polynomial.
p=cl*ones(size(x));
for j=l-1:-1:0,
% Check the exponent of x for current coefficient, px. If it is 0 or 1,
% just exit the loop
if(px<2), break; end;
% If current exponent is >=2, there is a "next" coefficient; multiply p
% by x2 and add it. Calculate the current coefficient
cl=-(j+j+2-l-m)*(j+j+1-l-m)/(2*(j+j+1)*(l-j))*cl;
if(maxcf<abs(cl)), maxcf=abs(cl); end;
%fprintf('Coef is %.16f\n',cl);
% ...and add to the polynomial
p=p.*x2 + cl;
% Decrease the exponent of x - this is the exponent of x corresponding
% to the newly added coefficient
px=px-2;
end;
% Estimate the error
err=maxcf*eps;
%fprintf('Coef is %.16f, err %.16f\n',maxcf, err);
% Now we're done adding coefficients. However, if the exponent of x
% corresponding to the last added coefficient is 1 (polynomial is odd),
% multiply the polynomial by x
if(px==1), p=p.*x; end;
% All that's left is to multiply the whole thing with sqrt(1-x^2)^(m/2). No
% further calculations are needed if m=0.
if(m==0), return; end;
x2=1-x2;
%First, multiply by the integer part of m/2
for j=1:floor(m/2), p=p.*x2; end;
%If m is odd, there is an additional factor sqrt(1-x^2)
if(m~=2*floor(m/2)), p=p.*sqrt(x2); end;
% Finally, the polynomials are not defined for |x|>1. If you do not need
% this behaviour, comment the following line
p(abs(x)>1)=NaN;
% Done.
c4science · Help