54 const unsigned int LMAX=c->coefficients.getLMax();
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());
83 P+=(c->coefficients(l,m)*Pn*exp(I*(m*phi)));
85 if (m!=0) P+= ( (m%2 ?-1.0:1.0)*c->coefficients(l,-m)*Pn*exp(-I*(m*phi)));
90 if (c->type==
MAGSQ)
return norm(P);
91 if (c->type==
MAG)
return abs(P);
92 if (c->type==
REAL)
return real(P);
93 if (c->type==
IMAG)
return imag(P);
94 if (!finite(retVal)) {
95 throw std::runtime_error(
"Non-finite return value in SphericalHarmonicExpansion");