1 /*
2 * $Id: legndr.i,v 1.1 2005-09-18 22:06:00 dhmunro Exp $
3 * Compute Legendre polynomials and associated Legendre functions.
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 /* After Numerical Recipes, by Press, et. al., section 6.6. */
12
legndr(l,m,x)13 func legndr(l,m, x)
14 /* DOCUMENT legndr(l,m, x)
15 return the associated Legendre function Plm(x). The X may
16 be an array (-1<=x<=1), but L and M (0<=M<=L) must be scalar
17 values. For m=0, these are the Legendre polynomials Pl(x).
18 Relation of Plm(x) to Pl(x):
19 Plm(x) = (-1)^m (1-x^2)^(m/2) d^m/dx^m(Pl(x))
20 Relation of Plm(x) to spherical harmonics Ylm:
21 Ylm(theta,phi)= sqrt((2*l+1)(l-m)!/(4*pi*(l+m)!)) *
22 Plm(cos(theta)) * exp(1i*m*phi)
23 SEE ALSO: ylm_coef
24 */
25 {
26 m+= 0;
27 l+= 0;
28 if (structof(m)!=long || structof(l)!=long ||
29 m<0 || m>l) error, "l, m integers with 0<=m<=l for Plm";
30
31 if (m > 0) {
32 /* sqrt blows up if x out of range */
33 pmm= ((m%2? -1:1) * exp(sum(log(indgen(1:2*m:2))))) *
34 (sqrt((1.-x)*(1.+x))^m);
35 } else {
36 if (anyof(abs(x)>1.0)) error, "-1<=x<=1 for Plm";
37 pmm= array(1.0, dimsof(x));
38 }
39
40 if (l==m) return pmm;
41
42 pmmp1= (2*m+1)*x*pmm;
43 if (l==m+1) return pmmp1;
44
45 for (ll=m+2 ; ; ++ll) {
46 c1= double(ll-m);
47 c= (ll+m-1)/c1;
48 c1= (2*ll-1)/c1;
49 pll= c1*x*pmmp1 - c*pmm;
50 if (ll==l) return pll;
51 pmm= pmmp1;
52 pmmp1= pll;
53 }
54 }
55
ylm_coef(l,m)56 func ylm_coef(l,m)
57 /* DOCUMENT ylm_coef(l,m)
58 return sqrt((2*l+1)(l-m)!/(4*pi*(l+m)!)), the normalization
59 coefficient for spherical harmonic Ylm with respect to the
60 associated Legendre function Plm. In this implementation,
61 0<=m<=l; use symmetry for m<0, or use sines and cosines
62 instead of complex exponentials. Unlike Plm, array L and M
63 arguments are permissible here.
64 WARNING: These get combinitorially small with large L and M;
65 probably Plm is simultaneously blowing up and should be
66 normalized directly in legndr if what you want is Ylm. But
67 I don't feel like working all that out -- if you need large
68 L and M results, you should probably be working with some
69 sort of asymptotic form anyway...
70 SEE ALSO: legndr
71 */
72 {
73 m+= 0;
74 l+= 0;
75 if (structof(m)!=long || structof(l)!=long ||
76 anyof(m<0) || anyof(m>l))
77 error, "l, m integers with 0<=m<=l; use symmetry for m<0";
78
79 lm= l-m;
80 not_scalar= dimsof(lm)(1);
81 if (not_scalar) result= array(1.0, dimsof(lm));
82 else result= [1.0];
83 iactive= indgen(numberof(lm));
84 for (list=where(lm) ; numberof(list) ; list=where(lm)) {
85 lm= lm(list);
86 iactive= iactive(list);
87 result(iactive)*= lm;
88 lm-= 1;
89 }
90
91 lm= l+m;
92 iactive= indgen(numberof(lm));
93 for (list=where(lm) ; numberof(list) ; list=where(lm)) {
94 lm= lm(list);
95 iactive= iactive(list);
96 result(iactive)/= lm;
97 lm-= 1;
98 }
99
100 result= sqrt(((2*l+1)/(4.*pi)) * result);
101 return not_scalar? result : result(1);
102 }
103