CLHEP/GenericFunctions/SphericalHarmonicExpansion.icc

0001 // -*- C++ -*-
0002 // $Id:
0003 #include <sstream>
0004 #include <cmath>
0005 #include <gsl/gsl_sf_legendre.h>
0006 #include <complex>
0007 #include <cstdlib>
0008 #include <sstream>
0009 #include <stdexcept>
0010 namespace Genfun {
0011
0012   FUNCTION_OBJECT_IMP(SphericalHarmonicExpansion)
0013
0014   class SphericalHarmonicExpansion::Clockwork {
0015
0016   public:
0017
0018     Clockwork(SphericalHarmonicExpansion::Type type, const SphericalHarmonicCoefficientSet & coefficients):type(type),coefficients(coefficients){}
0019
0020     SphericalHarmonicExpansion::Type type;
0021     SphericalHarmonicCoefficientSet coefficients;
0022
0023   };
0024
0025
0026   inline
0027   SphericalHarmonicExpansion::SphericalHarmonicExpansion(Type type, const SphericalHarmonicCoefficientSet & coefficients):
0028     c(new Clockwork(type,coefficients))
0029   {
0030
0031   }
0032
0033
0034   inline
0035   SphericalHarmonicExpansion::~SphericalHarmonicExpansion() {
0036     delete c;
0037   }
0038
0039   inline
0040   SphericalHarmonicExpansion::SphericalHarmonicExpansion(const SphericalHarmonicExpansion & right):
0041     AbsFunction(),
0042     c(new Clockwork(right.c->type,right.c->coefficients))
0043   {
0044   }
0045
0046   inline
0047   double SphericalHarmonicExpansion::operator() (double ) const {
0048     throw std::runtime_error("Dimensionality error in SphericalHarmonicExpansion");
0049     return 0;
0050   }
0051
0052   inline
0053   double SphericalHarmonicExpansion::operator() (const Argument & a ) const {
0054     const unsigned int LMAX=c->coefficients.getLMax();
0055     double x = a[0];
0056     double phi=a[1];
0057
0058     // Note, the calling sequence of the GSL Special Function forces us to
0059     // transpose Plm from its "natural" order.. It is addressed as P[m][l].
0060
0061     //double Plm[LMAX+1][LMAX+1];
0062     std::vector< std::vector<double> > Plm(LMAX+1);
0063     for (int m=0;m<=int(LMAX);m++) {
0064       Plm[m].resize(LMAX+1);
0065       gsl_sf_legendre_sphPlm_array (LMAX, m, x, &*Plm[m].begin());
0066       //gsl_sf_legendre_sphPlm_array (LMAX, m, x, Plm[m]);
0067     }
0068
0069     std::complex<double> P=0.0;
0070     std::complex<double> I(0,1.0);
0071     for (unsigned int l=0;l<=LMAX;l++) {
0072       for (int m=0;m<=int(l);m++) {
0073     {
0074       int LP=l-abs(m);
0075       double Pn= Plm[abs(m)][LP];
0076
0077       if (!finite(Pn)) {
0078         std::ostringstream stream;
0079         stream << "Non-finite Pn(x=" << x << ")";
0080         throw std::runtime_error(stream.str());
0081       }
0082       // Once for positive m (in all cases):
0083       P+=(c->coefficients(l,m)*Pn*exp(I*(m*phi)));
0084       // Once for negative m (skip if m==0);
0085       if (m!=0) P+= ( (m%2 ?-1.0:1.0)*c->coefficients(l,-m)*Pn*exp(-I*(m*phi)));
0086     }
0087       }
0088     }
0089     double retVal=0;
0090     if (c->type==MAGSQ) return norm(P);
0091     if (c->type==MAG)   return abs(P);
0092     if (c->type==REAL)  return real(P);
0093     if (c->type==IMAG)  return imag(P);
0094     if (!finite(retVal)) {
0095       throw std::runtime_error("Non-finite return value in SphericalHarmonicExpansion");
0096     }
0097     return retVal;
0098   }
0099
0100   inline
0101   const SphericalHarmonicCoefficientSet & SphericalHarmonicExpansion::coefficientSet() const {
0102     return c->coefficients;
0103   }
0104
0105   inline
0106   SphericalHarmonicCoefficientSet & SphericalHarmonicExpansion::coefficientSet() {
0107     return c->coefficients;
0108   }
0109
0110 } // end namespace Genfun