CLHEP/GenericFunctions/SphericalHarmonicCoefficientSet.icc

0001 #include <vector>
0002 #include "CLHEP/GenericFunctions/ClebschGordanCoefficientSet.hh"
0003 #include <stdexcept>
0004 namespace Genfun {
0005
0006   class SphericalHarmonicCoefficientSet::Clockwork {
0007
0008   public:
0009
0010     std::vector<std::vector<std::complex<double> > > data;
0011
0012   };
0013
0014   inline
0015   SphericalHarmonicCoefficientSet::SphericalHarmonicCoefficientSet(unsigned int LMAX):c(new Clockwork()){
0016     for (unsigned int l=0;l<=LMAX;l++) {
0017       std::vector<std::complex<double> > theMs;
0018       for (int m=-l; m<=int(l);m++) {
0019     theMs.push_back(std::complex<double> (0.0));
0020       }
0021       c->data.push_back(theMs);
0022     }
0023   }
0024
0025   inline
0026   SphericalHarmonicCoefficientSet::~SphericalHarmonicCoefficientSet(){
0027     delete c;
0028   }
0029
0030   inline
0031   SphericalHarmonicCoefficientSet::SphericalHarmonicCoefficientSet(const SphericalHarmonicCoefficientSet & right):
0032     c(new Clockwork(*right.c))
0033   {
0034   }
0035
0036   inline
0037   unsigned int SphericalHarmonicCoefficientSet::getLMax() const {
0038     return c->data.size()-1;
0039   }
0040
0041   inline
0042   const std::complex<double> &SphericalHarmonicCoefficientSet:: operator () (unsigned int l, int m) const {
0043     return c->data[l][m+l];
0044   }
0045
0046   inline
0047   std::complex<double> & SphericalHarmonicCoefficientSet::operator () (unsigned int l, int m) {
0048     return c->data[l][m+l];
0049   }
0050
0051   inline
0052   std::ostream & operator << ( std::ostream & o, const SphericalHarmonicCoefficientSet & c)
0053   {
0054     for (unsigned int l=0;l<=c.getLMax();l++) {
0055       for (int m=-l;m<=int(l);m++) {
0056     o << "l=" << l << " m=" ;
0057     if (m==0) o << " ";
0058     if (m>0 ) o << "+";
0059     o << m << " mag: " << c(l,m) << std::endl;
0060       }
0061       o << std::endl;
0062     }
0063     return o;
0064   }
0065
0066   inline
0067   SphericalHarmonicCoefficientSet & SphericalHarmonicCoefficientSet::operator= (const SphericalHarmonicCoefficientSet & source ){
0068     if (this!=&source) {
0069       delete c;
0070       c = new Clockwork(*source.c);
0071     }
0072     return *this;
0073   }
0074
0075
0076
0077
0078   inline
0079   SphericalHarmonicCoefficientSet & SphericalHarmonicCoefficientSet::operator*= (const std::complex<double> & s ){
0080     unsigned int LMAX=getLMax();
0081     for (unsigned int l=0;l<=LMAX;l++) {
0082       for (int m=-l;m<=int(l);m++) {
0083     operator()(l,m)*=s;
0084       }
0085     }
0086     return *this;
0087   }
0088
0089
0090   inline
0091   SphericalHarmonicCoefficientSet & SphericalHarmonicCoefficientSet::operator+= (const SphericalHarmonicCoefficientSet & source ){
0092     unsigned int LMAX=getLMax();
0093     for (unsigned int l=0;l<=LMAX;l++) {
0094       for (int m=-l;m<=int(l);m++) {
0095     operator()(l,m)+=source(l,m);
0096       }
0097     }
0098     return *this;
0099   }
0100
0101
0102   inline
0103   SphericalHarmonicCoefficientSet & SphericalHarmonicCoefficientSet::operator-= (const SphericalHarmonicCoefficientSet & source ){
0104     unsigned int LMAX=getLMax();
0105     for (unsigned int l=0;l<=LMAX;l++) {
0106       for (int m=-l;m<=int(l);m++) {
0107     operator()(l,m)-=source(l,m);
0108       }
0109     }
0110     return *this;
0111   }
0112
0113
0114
0115   inline
0116   std::complex<double> dot(const SphericalHarmonicCoefficientSet &a,
0117                const SphericalHarmonicCoefficientSet &b) {
0118
0119     std::complex<double> result=0.0;
0120     if (a.getLMax()!=b.getLMax()) throw std::runtime_error ("function dot:  SphericalHarmonicCoefficientSets of different dimension");
0121
0122     for (unsigned int l=0;l<=a.getLMax();l++) {
0123       for (int m=-l;m<=int(l);m++) {
0124     result += a(l,m)*conj(b(l,m));
0125       }
0126     }
0127     return result;
0128   }
0129
0130   inline SphericalHarmonicCoefficientSet squareExpansionCoefficients(const SphericalHarmonicCoefficientSet & coefficientsA) {
0131     unsigned int LMAX=coefficientsA.getLMax();
0132     SphericalHarmonicCoefficientSet coefficientsASq(2*LMAX);
0133     static ClebschGordanCoefficientSet clebschGordan;
0134     for (unsigned int L=0;L<=2*LMAX;L++) {
0135       for (int M=-L; M<=int(L); M++) {
0136     coefficientsASq(L,M)=0.0;
0137     for (unsigned int l1=0;l1<=LMAX;l1++) {
0138       for (unsigned int l2=0;l2<=LMAX;l2++) {
0139         for (int m1=-l1;m1<=int(l1);m1++) {
0140           for (int m2=-l2;m2<=int(l2);m2++) {
0141         if (m1-m2==M) {
0142           if (((l1+l2) >= L) && abs(l1-l2) <= int(L))  {
0143             coefficientsASq(L,M) += (coefficientsA(l1,m1)*
0144                          conj(coefficientsA(l2,m2))*
0145                          (m2%2 ? -1.0:1.0) *
0146                          sqrt((2*l1+1)*(2*l2+1)/(4*M_PI*(2*L+1)))*
0147                          clebschGordan(l1,l2,0,0,L,0)*clebschGordan(l1,l2,m1,-m2,L,M));
0148           }
0149         }
0150           }
0151         }
0152       }
0153     }
0154       }
0155     }
0156     return coefficientsASq;
0157   }
0158
0159 }