Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
161 lines (140 sloc) 4.14 KB
function varargout=legendreint01(l,m,x0)
% val=legendreint01(l,m,x0)
%
% Evaluates integrals of a Schmidt semi-normalized real Legendre
% polynomial P_lm(x)dx between x0 to 1 for l=0 or 1, analytically.
%
% INPUT:
%
% l Angular degree of the polynomial, l>=0
% m Angular order of the polynomial, 0<=m<=l
% x0 Lower bound(s) for the integral(s)
%
% OUTPUT:
%
% val Integral(s) of LEGENDRE(N,x,'sch') or LIBBRECHT(N,x,'sch')
%
% SEE ALSO:
%
% LIBBRECHT, PAUL, LEGENREPRODINT
%
% EXAMPLE:
%
% legendreint01('demo1') Compare analytic to numerical Legendre functions
% legendreint01('demo2') Compare Gauss-Legendre integration to analytic
% legendreint01('demo3') Evaluates Gauss-Legendre integration errors
%
% Last modified by plattner-at-princeton.edu, 05/17/2011
% Last modified by fjsimons-at-alum.mit.edu, 06/01/2011
% Calculates the integral of a single Legendre function between x0 and 1
if ~isstr(l)
if l<0 || m<0 || m>l
error('Bad choice of l and m');
end
if m>1
error('Use PAUL or LEGENDREPRODINT instead.')
end
if l==0
% l=0 m=0, the integral of 1
val=(1-x0);
elseif l==1
if m==0
% l=1 m=0, the integral of x
val=(1-x0.^2)/2;
else
% l=1 m=1, the integral of sqrt(1-x.^2)
% See Research Book 9 page 79 for these alternative expressions
%val=( +acos(x0)/2-(x0.*sqrt(1-x0.^2))/2);
%val= coscos([],1,1,[asin(x0(:)) repmat(pi/2,length(x0),1)])';
%val=-sinsin([],1,1,[acos(x0(:)) repmat(0,length(x0),1)])';
val=(pi/4-asin(x0)/2-(x0.*sqrt(1-x0.^2))/2);
end
end
% Optional output
varns={val};
varargout=varns(1:nargout);
elseif strcmp(l,'demo1')
x=linspace(-1,1,100);
clf
subplot(311)
plot(x,ones(size(x)),'k')
hold on
plot(x,legendre(0,x,'sch'),'o')
plot(x,libbrecht(0,x,'sch'),'+')
subplot(312)
plot(x,x,'k')
hold on
plot(x,rindeks(legendre(1,x,'sch'),1),'o')
plot(x,rindeks(libbrecht(1,x,'sch'),1),'+')
subplot(313)
plot(x,sqrt(1-x.^2),'k')
hold on
plot(x,rindeks(legendre(1,x,'sch'),2),'o')
plot(x,rindeks(libbrecht(1,x,'sch'),2),'+')
elseif strcmp(l,'demo2')
x=linspace(-1,1,100);
ngl=5;
integrand1=inline('ones(size(x))');
integrand2=inline('x');
integrand3=inline('sqrt(1-x.^2)');
[gl1,gl2,gl3]=deal(zeros(size(x)));
for i=1:length(x)
% Here by GL integration
gl1(i)=gausslegendre([x(i) 1],integrand1,ngl);
gl2(i)=gausslegendre([x(i) 1],integrand2,ngl);
gl3(i)=gausslegendre([x(i) 1],integrand3,ngl);
end
% Now the code itself
val1=legendreint01(0,0,x);
val2=legendreint01(1,0,x);
val3=legendreint01(1,1,x);
% And another piece of code at a random point
randi=ceil(rand*length(x));
opti={'automatic','dumb','gl','paul'};
% Note that the choice may be overridden inside LEGENDREPRODINT
rando=ceil(rand*length(opti));
% Use various algorithms explicitly at random
vol1=legendreprodint(0,0,0,0,x(randi),opti{rando});
vol2=legendreprodint(1,0,0,0,x(randi),opti{rando});
vol3=legendreprodint(1,1,0,0,x(randi),opti{rando});
clf
subplot(311)
plot(x,val1,'k')
hold on
plot(x,gl1,'b+')
plot(x(randi),vol1,'rv','MarkerF','r','markers',6)
subplot(312)
plot(x,val2,'k')
hold on
plot(x,gl2,'b+')
plot(x(randi),vol2,'rv','MarkerF','r','markers',6)
subplot(313)
plot(x,val3,'k')
hold on
plot(x,gl3,'b+')
plot(x(randi),vol3,'rv','MarkerF','r','markers',6)
elseif strcmp(l,'demo3')
x=-1;
ngl=linspace(1,1000,25);
integrand1=inline('ones(size(x))');
integrand2=inline('x');
integrand3=inline('sqrt(1-x.^2)');
[gl1,gl2,gl3]=deal(zeros(size(x)));
for i=1:length(ngl)
% Here by GL integration
gl1(i)=gausslegendre([x 1],integrand1,ngl(i));
gl2(i)=gausslegendre([x 1],integrand2,ngl(i));
gl3(i)=gausslegendre([x 1],integrand3,ngl(i));
end
% Now the code itself
val1=legendreint01(0,0,x);
val2=legendreint01(1,0,x);
val3=legendreint01(1,1,x);
clf
subplot(311)
semilogy(ngl,abs(gl1-val1),'k')
subplot(312)
semilogy(ngl,abs(gl2-val2),'k')
subplot(313)
semilogy(ngl,abs(gl3-val3),'k')
end