csdms-contrib/slepian_alpha

Switch branches/tags
Nothing to show
Fetching contributors…
Cannot retrieve contributors at this time
106 lines (100 sloc) 2.69 KB
 function [PpL,PLm1,PL,PLp1]=legendrediff(L,x,norma) % [PpL,PLm1,PL,PLp1]=legendrediff(L,x,norma) % % Computes the derivative of the Legendre polynomial with respect to its % argument x=cos(theta). Compare LIBBRECHT, where it is wrt acos(x). % % INPUT: % % L Degree of spherical harmonic < 256 % x Evaluation point(s) % norma 'sch' Schmidt-normalized [default] % 'fnr' Fully normalized real % % OUTPUT: % % PpL Schmidt-normalized derivative of P_L(x), m=0 % PLm1 Schmidt-normalized polynomial P_{L-1}(x), m=0 % PL Schmidt-normalized polynomial P_L(x), m=0 % PLp1 Schmidt-normalized polynomial P_{L+1}(x), m=0 % % See Wolfram under Legendre-Gauss Quadrature. % % EXAMPLE: % % legendrediff('demo1') % Comparison with LIBBRECHT and DIFF % % SEE ALSO: LIBBRECHT, YLM % % Last modified by fjsimons-at-alum.mit.edu, 05/17/2011 defval('L','demo1') defval('norma','sch') if ~isstr(L) switch norma case 'sch' fac1=1; fac2=1; fac3=1; case 'fnr' fac1=sqrt(2*(L-1)+1); fac2=sqrt(2*L+1); fac3=sqrt(2*(L+1)+1); otherwise error('Specify valid normalization') end % Must build in what it means to be -1 PLm1=rindeks(legendre(L-1,x,'sch'),1); PL=rindeks(legendre(L,x,'sch'),1); warning off PpL=L*(x(:).*PL(:)-PLm1(:))./(x(:).^2-1); warning on % From Boyd (2001) PpL(x==1)=L*(L+1)/2; PpL(x==-1)=(-1)^(L-1)*L*(L+1)/2; PpL=fac2*PpL(:)'; PLm1=PLm1*fac1; PL=PL*fac2; if nargout==4 PLp1=fac3*rindeks(legendre(L+1,x,'sch'),1); end elseif strcmp(L,'demo1') clf theta=linspace(0,pi,500); x=cos(theta); more off deg=round(rand*20);deg=9 m=0; p=rindeks(legendre(deg,x,'sch'),1); [pp,jk,p2]=legendrediff(deg,x,'sch'); [pl,dpl]=libbrecht(deg,x,'sch',[]); pl=rindeks(pl,1); dpl=rindeks(dpl,1); subplot(211) plot(x,p,'b-','LineW',2); hold on plot(x,pl,'y-','LineW',1); if max(abs(p2(:)-p(:)))>0; error('Something wrong?'); end title(sprintf(... 'Legendre functions (Schmidt); l=%i (m= %i)',... deg,m),'FontS',12) grid on; axis tight ; openup(gca,6); nolabels(gca,1) yl=ylabel('N_l^m\timesP_l^m(cos\theta)'); movev(gca,-.1) l1=legend('LEGENDRE','LIBBRECHT'); longticks(gca,2) subplot(212) plot(x,pp,'r-','LineW',2); hold on % Watch out since x is not equally spaced. plot(x(2:end)-indeks(diff(x),1)/2,... diff(p)./diff(x),'y') warning off plot(x,-dpl./sin(theta),'k--','LineW',1); warning on yl=ylabel('dP_l^m(cos(\theta))/dcos\theta'); xl=xlabel('cos(\theta)'); axis tight ; grid on; openup(gca,6) l=legend('LEGENDREDIFF','DIFF','LIBBRECHT'); longticks(gca,2) fig2print(gcf,'portrait'); id; axes(l) figdisp end