Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
202 lines (171 sloc) 5.27 KB
function varargout=paul(Lmax,x0)
% [Itab,Itabo,dels,dems]=paul(Lmax,x0)
%
% Recursively evaluates a list of integrals of a single Schmidt
% semi-normalized real Legendre polynomial P_lm(x)dx from x0 to 1.
% The iteration is set up according to Paul (1978) 'Recurrence relations
% for integrals of associated Legendre functions'.
%
% INPUT:
%
% Lmax Bandwidth (maximum spherical harmonic degree)
% x0 Vector of lower bounds for the integrals
%
% OUTPUT:
%
% Itab List of integrals of LEGENDRE(l,x,'sch') from x0 to 1
% Itabo Same but without the sqrt(2-dom) factor subsumed by Matlab
%
% EXAMPLE:
%
% paul('demo1') Random Lmax, 100 values for x0, PAUL vs GAUSSLEGENDRE
% paul('demo2') Complete list comparison PAUL vs GAUSSLEGENDRE
%
% Last modified by plattner-at-princeton.edu, 05/17/2011
% Last modified by fjsimons-at-alum.mit.edu, 06/01/2011
defval('Lmax',ceil(rand*256))
defval('x0',rand(256,1)*2-1)
if ~isstr(Lmax)
% Turn the lower bounds into row
x0=x0(:)';
% first we test if Lmax is zero or positive
if(Lmax<0)
error('Lmax must be non-negative');
Itable=[];
end
% Initialize Itab
[Itab,Ptab]=deal(zeros(Lmax*(Lmax+1)/2+Lmax+1,length(x0)));
% Now the starting evaluations for the LEGENDRE(l,x,'sch')
Itab(1,:)=legendreint01(0,0,x0);
Itab(2,:)=legendreint01(1,0,x0);
% Get rid of the sqrt(2) factor
Itab(3,:)=legendreint01(1,1,x0)/sqrt(2);
% Now the Legendre functions themselves
Ptab(2:3,:)=legendre(1,x0,'sch');
% Get rid of the sqrt(2) factor
Ptab(3,:)=Ptab(3,:)/sqrt(2);
% This is a common argument
y2=(1-x0.^2);
% Now begin with the iterations
for l=2:Lmax
% This is the index for m=0 at this l, i.e. addmup(l-1)+1
lind=l*(l+1)/2+1;
% This is the index for m=0 at l-1, i.e. addmup(l-2)+1
lond=(l-1)*l/2+1;
% This is the index for m=0 at l-2, i.e. addmup(l-3)+1
lund=(l-2)*(l-1)/2+1;
% Calculate the associated Legendre functions for this l
Ptab(lind:lind+l,:)=legendre(l,x0,'sch');
% Get rid of the sqrt(2) factor
Ptab(lind+1:lind+l,:)=Ptab(lind+1:lind+l,:)/sqrt(2);
% Some common factors at this l
fux=(2*l-1)/(l+1);
fex=sqrt((2*l-1)*(2*l-2)*2*l);
% These are the cases where m~=l
for m=0:l-1
% Some common factors at this m
fax=sqrt((l +m)*(l -m));
fox=sqrt((l-1+m)*(l-1-m));
if m<l-1
Itab(lind+m,:)=Itab(lund+m,:)*(l-2)/(l+1)/fax*fox;
end
Itab(lind+m,:)=Itab(lind+m,:)+y2.*Ptab(lond+m,:)/fax*fux;
end
% This is the case where m=l
Itab(lind+l,:)=Itab(lund+l-2,:)*l*sqrt(2*l-3)/fex*fux-...
y2.*Ptab(lond+l-2,:)/fex*fux;
end
[dems,dels,mz]=addmon(Lmax);
% Now put the sqrt(2) factor in again since that's what we asked for
nonzonal=skip(1:length(dems),mz);
Itabo=Itab;
Itab(nonzonal,:)=Itab(nonzonal,:)*sqrt(2);
% Optional output
varns={Itab,Itabo,dels,dems};
varargout=varns(1:nargout);
elseif strcmp(Lmax,'demo1')
% Demo to compare the Paul iteration method to Gauss-Legendre integration
L=round(256*rand);
m=round(L*rand);
N=100;
disp(sprintf('Calculate integrals for l up to %i for %i values for x0',L,N))
x0=linspace(-1,1,N);
% Paul's (1978) method
tic
Itab=paul(L,x0);
toc
% Gauss Legendre integration of the associated Legendre polynomials
tic
%integrand=inline(sprintf(...
% ['rindeks(legendre(%i,x,''sch''),%i)'],...
% L,m+1));
integrand=@(x) rindeks(legendre(L,x,'sch'),m+1);
%[w,xgl,nsel]=gausslegendrecof(max(L,200*((m)~=0)));
[w,xgl,nsel]=gausslegendrecof(1000);
in=zeros(1,N);
for i=1:N
in(i)=gausslegendre([x0(i) 1],integrand,[w(:) xgl(:)]);
end
toc
subplot(211)
plot(x0,Itab(L*(L+1)/2+m+1,:),'r-+',x0,in,'b')
xlabel('Lower limit x0')
ylabel('I_{l,m}')
title(sprintf('l = %d, m = %d',L,m))
legend('paul','gausslegendre')
% Calculate the error
errs=Itab(L*(L+1)/2+m+1,:)-in;
difer(errs);
subplot(212)
plot(x0,errs,'k')
xlabel('Lower limit x0')
ylabel('Difference')
title(sprintf('l = %d, m = %d',L,m))
elseif strcmp(Lmax,'demo2')
L=round(rand*50);
N=10;
x0=-1:2/(N-1):1;
x0=x0(:)';
% Paul's method
tic
Itab=paul(L,x0);
toc
% Gauss Legendre integration of the associated Legendre polynomials
gleg=zeros(L*(L+1)/2+L+1,N);
[w,xgl,nsel]=gausslegendrecof(2000);
tic
for LL=0:L
for m=0:LL
% Slower...
%integrand=inline(sprintf(...
% ['rindeks(legendre(%i,x,''sch''),%i)'],...
% LL,m+1));
% Faster...
integrand=@(x) rindeks(legendre(LL,x,'sch'),m+1);
in=zeros(1,N);
for i=1:N
in(i)=gausslegendre([x0(i) 1],integrand,[w(:) xgl(:)]);
end
gleg(LL*(LL+1)/2+m+1,:)=in;
end
end
toc
% Calculate errors of everything
err=Itab-gleg;
difer(err,7)
% Setup error matrix
errs=zeros(L,L);
for LL=0:L
for m=0:LL
errs(LL+1,m+1)=sum(abs(gleg(LL*(LL+1)/2+m+1,:)...
-Itab(LL*(LL+1)/2+m+1,:)));
end
end
clf
imagefnan([0 0],[L L],errs,'kelicol',minmax(errs))
longticks(gca,2)
xlabel('order m'); ylabel('degree l');
[cb,xcb]=addcb('hor',minmax(errs),minmax(errs),'kelicol',max(errs(:)));
movev(cb,-.05)
set(xcb,'string',sprintf('sum of the absolute errors for %i lower limits',N))
end