/usr/include/CLHEP/GenericFunctions/SphericalHarmonicExpansion.icc is in libclhep-dev 2.1.4.1+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 | // -*- C++ -*-
// $Id:
#include <sstream>
#include <cmath>
#include <gsl/gsl_sf_legendre.h>
#include <complex>
#include <cstdlib>
#include <sstream>
#include <stdexcept>
namespace Genfun {
FUNCTION_OBJECT_IMP(SphericalHarmonicExpansion)
class SphericalHarmonicExpansion::Clockwork {
public:
Clockwork(SphericalHarmonicExpansion::Type type, const SphericalHarmonicCoefficientSet & coefficients):type(type),coefficients(coefficients){}
SphericalHarmonicExpansion::Type type;
SphericalHarmonicCoefficientSet coefficients;
};
inline
SphericalHarmonicExpansion::SphericalHarmonicExpansion(Type type, const SphericalHarmonicCoefficientSet & coefficients):
c(new Clockwork(type,coefficients))
{
}
inline
SphericalHarmonicExpansion::~SphericalHarmonicExpansion() {
delete c;
}
inline
SphericalHarmonicExpansion::SphericalHarmonicExpansion(const SphericalHarmonicExpansion & right):
AbsFunction(),
c(new Clockwork(right.c->type,right.c->coefficients))
{
}
inline
double SphericalHarmonicExpansion::operator() (double ) const {
throw std::runtime_error("Dimensionality error in SphericalHarmonicExpansion");
return 0;
}
inline
double SphericalHarmonicExpansion::operator() (const Argument & a ) const {
const unsigned int LMAX=c->coefficients.getLMax();
double x = a[0];
double phi=a[1];
// Note, the calling sequence of the GSL Special Function forces us to
// transpose Plm from its "natural" order.. It is addressed as P[m][l].
//double Plm[LMAX+1][LMAX+1];
std::vector< std::vector<double> > Plm(LMAX+1);
for (int m=0;m<=int(LMAX);m++) {
Plm[m].resize(LMAX+1);
gsl_sf_legendre_sphPlm_array (LMAX, m, x, &*Plm[m].begin());
//gsl_sf_legendre_sphPlm_array (LMAX, m, x, Plm[m]);
}
std::complex<double> P=0.0;
std::complex<double> I(0,1.0);
for (unsigned int l=0;l<=LMAX;l++) {
for (int m=0;m<=int(l);m++) {
{
int LP=l-abs(m);
double Pn= Plm[abs(m)][LP];
if (!finite(Pn)) {
std::ostringstream stream;
stream << "Non-finite Pn(x=" << x << ")";
throw std::runtime_error(stream.str());
}
// Once for positive m (in all cases):
P+=(c->coefficients(l,m)*Pn*exp(I*(m*phi)));
// Once for negative m (skip if m==0);
if (m!=0) P+= ( (m%2 ?-1.0:1.0)*c->coefficients(l,-m)*Pn*exp(-I*(m*phi)));
}
}
}
double retVal=0;
if (c->type==MAGSQ) return norm(P);
if (c->type==MAG) return abs(P);
if (c->type==REAL) return real(P);
if (c->type==IMAG) return imag(P);
if (!finite(retVal)) {
throw std::runtime_error("Non-finite return value in SphericalHarmonicExpansion");
}
return retVal;
}
inline
const SphericalHarmonicCoefficientSet & SphericalHarmonicExpansion::coefficientSet() const {
return c->coefficients;
}
inline
SphericalHarmonicCoefficientSet & SphericalHarmonicExpansion::coefficientSet() {
return c->coefficients;
}
} // end namespace Genfun
|