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