CLHEP 2.4.7.1
C++ Class Library for High Energy Physics
SphericalHarmonicExpansion.icc
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id:
3#include <sstream>
4#include <cmath>
5#include <gsl/gsl_sf_legendre.h>
6#include <complex>
7#include <cstdlib>
8#include <sstream>
9#include <stdexcept>
10namespace Genfun {
11
12 FUNCTION_OBJECT_IMP(SphericalHarmonicExpansion)
13
24
25
26 inline
28 c(new Clockwork(type,coefficients))
29 {
30
31 }
32
33
34 inline
38
39 inline
45
46 inline
48 throw std::runtime_error("Dimensionality error in SphericalHarmonicExpansion");
49 return 0;
50 }
51
52 inline
54 const unsigned int LMAX=c->coefficients.getLMax();
55 double x = a[0];
56 double phi=a[1];
57
58 // Note, the calling sequence of the GSL Special Function forces us to
59 // transpose Plm from its "natural" order.. It is addressed as P[m][l].
60
61 //double Plm[LMAX+1][LMAX+1];
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());
66 //gsl_sf_legendre_sphPlm_array (LMAX, m, x, Plm[m]);
67 }
68
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++) {
73 {
74 int LP=l-abs(m);
75 double Pn= Plm[abs(m)][LP];
76
77 if (!finite(Pn)) {
78 std::ostringstream stream;
79 stream << "Non-finite Pn(x=" << x << ")";
80 throw std::runtime_error(stream.str());
81 }
82 // Once for positive m (in all cases):
83 P+=(c->coefficients(l,m)*Pn*exp(I*(m*phi)));
84 // Once for negative m (skip if m==0);
85 if (m!=0) P+= ( (m%2 ?-1.0:1.0)*c->coefficients(l,-m)*Pn*exp(-I*(m*phi)));
86 }
87 }
88 }
89 double retVal=0;
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");
96 }
97 return retVal;
98 }
99
100 inline
104
105 inline
109
110} // end namespace Genfun
#define FUNCTION_OBJECT_IMP(classname)
Clockwork(SphericalHarmonicExpansion::Type type, const SphericalHarmonicCoefficientSet &coefficients)
SphericalHarmonicExpansion(Type type, const SphericalHarmonicCoefficientSet &coefficients)
SphericalHarmonicCoefficientSet & coefficientSet()
virtual double operator()(double argument) const override
Definition Abs.hh:14