52 |
|
ComplexType SphericalHarmonic::getValueAt(RealType costheta, RealType phi) { |
53 |
|
|
54 |
|
RealType p; |
55 |
– |
ComplexType phase; |
56 |
– |
ComplexType I(0.0, 1.0); |
55 |
|
|
56 |
|
// associated Legendre polynomial |
57 |
< |
p = Legendre(L, M, costheta); |
58 |
< |
|
59 |
< |
phase = exp(I * (ComplexType)M * (ComplexType)phi); |
60 |
< |
|
63 |
< |
return coefficient * phase * (ComplexType)p; |
57 |
> |
p = Ptilde(L, M, costheta); |
58 |
> |
ComplexType phase(0.0, (RealType)M * phi); |
59 |
> |
|
60 |
> |
return exp(phase) * (ComplexType)p; |
61 |
|
|
62 |
|
} |
66 |
– |
|
67 |
– |
//---------------------------------------------------------------------------// |
63 |
|
// |
64 |
< |
// RealType LegendreP (int l, int m, RealType x); |
64 |
> |
// Routine to calculate the associated Legendre polynomials for m>=0 |
65 |
|
// |
66 |
< |
// Computes the value of the associated Legendre polynomial P_lm (x) |
67 |
< |
// of order l at a given point. |
68 |
< |
// |
69 |
< |
// Input: |
70 |
< |
// l = degree of the polynomial >= 0 |
71 |
< |
// m = parameter satisfying 0 <= m <= l, |
72 |
< |
// x = point in which the computation is performed, range -1 <= x <= 1. |
73 |
< |
// Returns: |
79 |
< |
// value of the polynomial in x |
80 |
< |
// |
81 |
< |
//---------------------------------------------------------------------------// |
82 |
< |
RealType SphericalHarmonic::LegendreP (int l, int m, RealType x) { |
83 |
< |
// check parameters |
84 |
< |
if (m < 0 || m > l || fabs(x) > 1.0) { |
85 |
< |
printf("LegendreP got a bad argument: l = %d\tm = %d\tx = %lf\n", l, m, x); |
66 |
> |
RealType SphericalHarmonic::LegendreP(int l,int m, RealType x) { |
67 |
> |
|
68 |
> |
RealType temp1, temp2, temp3, temp4, result; |
69 |
> |
RealType temp5; |
70 |
> |
int i, ll; |
71 |
> |
|
72 |
> |
if (fabs(x) > 1.0) { |
73 |
> |
printf("LegendreP: x out of range: l = %d\tm = %d\tx = %lf\n", l, m, x); |
74 |
|
return std::numeric_limits <RealType>:: quiet_NaN(); |
75 |
|
} |
76 |
|
|
77 |
< |
RealType pmm = 1.0; |
78 |
< |
if (m > 0) { |
79 |
< |
RealType h = sqrt((1.0-x)*(1.0+x)), |
92 |
< |
f = 1.0; |
93 |
< |
for (int i = 1; i <= m; i++) { |
94 |
< |
pmm *= -f * h; |
95 |
< |
f += 2.0; |
96 |
< |
} |
77 |
> |
if (m>l) { |
78 |
> |
printf("LegendreP: m > l: l = %d\tm = %d\tx = %lf\n", l, m, x); |
79 |
> |
return std::numeric_limits <RealType>:: quiet_NaN(); |
80 |
|
} |
81 |
< |
if (l == m) |
82 |
< |
return pmm; |
83 |
< |
else { |
84 |
< |
RealType pmmp1 = x * (2 * m + 1) * pmm; |
85 |
< |
if (l == (m+1)) |
86 |
< |
return pmmp1; |
87 |
< |
else { |
88 |
< |
RealType pll = 0.0; |
89 |
< |
for (int ll = m+2; ll <= l; ll++) { |
90 |
< |
pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m); |
91 |
< |
pmm = pmmp1; |
92 |
< |
pmmp1 = pll; |
81 |
> |
|
82 |
> |
if (m<0) { |
83 |
> |
printf("LegendreP: m < 0: l = %d\tm = %d\tx = %lf\n", l, m, x); |
84 |
> |
return std::numeric_limits <RealType>:: quiet_NaN(); |
85 |
> |
} else { |
86 |
> |
temp3=1.0; |
87 |
> |
|
88 |
> |
if (m>0) { |
89 |
> |
temp1=sqrt(1.0-pow(x,2)); |
90 |
> |
temp5 = 1.0; |
91 |
> |
for (i=1;i<=m;++i) { |
92 |
> |
temp3 *= -temp5*temp1; |
93 |
> |
temp5 += 2.0; |
94 |
|
} |
111 |
– |
return pll; |
95 |
|
} |
96 |
+ |
if (l==m) { |
97 |
+ |
result = temp3; |
98 |
+ |
} else { |
99 |
+ |
temp4=x*(2.*m+1.)*temp3; |
100 |
+ |
if (l==(m+1)) { |
101 |
+ |
result = temp4; |
102 |
+ |
} else { |
103 |
+ |
for (ll=(m+2);ll<=l;++ll) { |
104 |
+ |
temp2 = (x*(2.*ll-1.)*temp4-(ll+m-1.)*temp3)/(RealType)(ll-m); |
105 |
+ |
temp3=temp4; |
106 |
+ |
temp4=temp2; |
107 |
+ |
} |
108 |
+ |
result = temp2; |
109 |
+ |
} |
110 |
+ |
} |
111 |
|
} |
112 |
+ |
return result; |
113 |
|
} |
114 |
|
|
115 |
+ |
|
116 |
|
// |
117 |
|
// Routine to calculate the associated Legendre polynomials for all m... |
118 |
|
// |
124 |
|
} else if (m >= 0) { |
125 |
|
result = LegendreP(l,m,x); |
126 |
|
} else { |
127 |
+ |
//result = mpow(-m)*LegendreP(l,-m,x); |
128 |
|
result = mpow(-m)*Fact(l+m)/Fact(l-m)*LegendreP(l, -m, x); |
129 |
|
} |
130 |
|
result *=mpow(m); |
131 |
|
return result; |
132 |
|
} |
133 |
|
// |
134 |
+ |
// Routine to calculate the normalized associated Legendre polynomials... |
135 |
+ |
// |
136 |
+ |
RealType SphericalHarmonic::Ptilde(int l,int m, RealType x){ |
137 |
+ |
|
138 |
+ |
RealType result; |
139 |
+ |
if (m>l || m<-l) { |
140 |
+ |
result = 0.; |
141 |
+ |
} else { |
142 |
+ |
RealType y=(RealType)(2.*l+1.)*Fact(l-m)/Fact(l+m); |
143 |
+ |
result = sqrt(y) * Legendre(l,m,x); |
144 |
+ |
} |
145 |
+ |
return result; |
146 |
+ |
} |
147 |
+ |
// |
148 |
|
// mpow returns (-1)**m |
149 |
|
// |
150 |
|
RealType SphericalHarmonic::mpow(int m) { |