CLHEP/GenericFunctions/LegendreExpansion.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 <stdexcept>
0009 namespace Genfun {
0010
0011   FUNCTION_OBJECT_IMP(LegendreExpansion)
0012
0013   class LegendreExpansion::Clockwork {
0014
0015   public:
0016
0017     Clockwork(LegendreExpansion::Type type, const LegendreCoefficientSet & coefficients):type(type),coefficients(coefficients){}
0018
0019     LegendreExpansion::Type type;
0020     LegendreCoefficientSet coefficients;
0021
0022   };
0023
0024
0025   inline
0026   LegendreExpansion::LegendreExpansion(Type type, const LegendreCoefficientSet & coefficients):
0027     c(new Clockwork(type,coefficients))
0028   {
0029
0030   }
0031
0032
0033   inline
0034   LegendreExpansion::~LegendreExpansion() {
0035     delete c;
0036   }
0037
0038   inline
0039   LegendreExpansion::LegendreExpansion(const LegendreExpansion & right):
0040     AbsFunction(),
0041     c(new Clockwork(right.c->type,right.c->coefficients))
0042   {
0043   }
0044
0045   inline
0046   double LegendreExpansion::operator() (double x) const {
0047
0048     int N=c->coefficients.getLMax();
0049     std::vector<double> Pk(N+1);
0050     gsl_sf_legendre_Pl_array(N, x, &Pk[0]);
0051     unsigned int n=N;
0052     std::complex<double> P=0.0;
0053     std::complex<double> I(0,1.0);
0054     while (1) {
0055       if (n==0) {
0056     P+=c->coefficients(n)*Pk[n];
0057     break;
0058       }
0059       else {
0060     P+=c->coefficients(n)*Pk[n];
0061     n--;
0062       }
0063     }
0064
0065     double retVal=0;
0066     if (c->type==MAGSQ) return norm(P);
0067     if (c->type==MAG)   return abs(P);
0068     if (c->type==REAL)  return real(P);
0069     if (c->type==IMAG)  return imag(P);
0070     if (!finite(retVal)) {
0071       throw std::runtime_error("Non-finite return value in LegendreExpansion");
0072     }
0073     return retVal;
0074   }
0075
0076   inline
0077   const LegendreCoefficientSet & LegendreExpansion::coefficientSet() const {
0078     return c->coefficients;
0079   }
0080
0081 } // end namespace Genfun