62 std::vector< std::vector<double> > Plm(LMAX+1);
63 for (
int m=0;m<=int(LMAX);m++) {
64 Plm[m].resize(LMAX+1);
65 gsl_sf_legendre_sphPlm_array (LMAX, m, x, &*Plm[m].begin());
69 std::complex<double> P=0.0;
70 std::complex<double> I(0,1.0);
71 for (
unsigned int l=0;l<=LMAX;l++) {
72 for (
int m=0;m<=int(l);m++) {
75 double Pn= Plm[abs(m)][LP];
78 std::ostringstream stream;
79 stream <<
"Non-finite Pn(x=" << x <<
")";
80 throw std::runtime_error(stream.str());
85 if (m!=0) P+= ( (m%2 ?-1.0:1.0)*c->
coefficients(l,-m)*Pn*exp(-I*(m*phi)));
91 if (c->
type==
MAG)
return abs(P);
94 if (!finite(retVal)) {
95 throw std::runtime_error(
"Non-finite return value in SphericalHarmonicExpansion");