Geant4  10.02.p01
G4LegendrePolynomial.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 
26 #include "G4ios.hh"
27 #include "G4LegendrePolynomial.hh"
28 #include "G4Pow.hh"
29 #include "G4Exp.hh"
30 #include "G4Log.hh"
31 
32 using namespace std;
33 
35 {
36  if(order >= fCoefficients.size()) BuildUpToOrder(order);
37  if(order >= fCoefficients.size() ||
38  i/2 >= fCoefficients[order].size() ||
39  (i%2) != order %2) return 0;
40  return fCoefficients[order][i/2];
41 }
42 
44 {
45  // Call EvalAssocLegendrePoly with m=0
46  return (EvalAssocLegendrePoly(order,0,x));
47 }
48 
50 {
51  if(l<0 || m<-l || m>l) return 0;
52  G4Pow* g4pow = G4Pow::GetInstance();
53 
54  // Use non-log factorial, which is more efficient until l and m get
55  // above 10 or so.
56  // FIXME: G4Pow doesn't check whether the argument gets too large,
57  // which is unsafe! Max is 512; VI: It is assume that Geant4 does not
58  // need higher order
59  if(m<0) return (m%2 ? -1. : 1.) * g4pow->factorial(l+m)/g4pow->factorial(l-m)
60  * EvalAssocLegendrePoly(l, -m, x);
61 
62  // hard-code the first few orders for speed
63  if(l==0) return 1;
64  if(l==1) {
65  if(m==0) return x;
66  /*m==1*/ return -sqrt(1.-x*x);
67  }
68  if(l<5) {
69  G4double x2 = x*x;
70  if(l==2) {
71  if(m==0) return 0.5*(3.*x2 - 1.);
72  if(m==1) return -3.*x*sqrt(1.-x2);
73  /*m==2*/ return 3.*(1.-x2);
74  }
75  if(l==3) {
76  if(m==0) return 0.5*(5.*x*x2 - 3.*x);
77  if(m==1) return -1.5*(5.*x2-1.)*sqrt(1.-x2);
78  if(m==2) return 15.*x*(1.-x2);
79  /*m==3*/ return -15.*(1.-x2)*sqrt(1.-x2);
80  }
81  if(l==4) {
82  if(m==0) return 0.125*(35.*x2*x2 - 30.*x2 + 3.);
83  if(m==1) return -2.5*(7.*x*x2-3.*x)*sqrt(1.-x2);
84  if(m==2) return 7.5*(7.*x2-1.)*(1.-x2);
85  if(m==3) return -105.*x*(1.-x2)*sqrt(1.-x2);
86  /*m==4*/ return 105.*(1. - 2.*x2 + x2*x2);
87  }
88  }
89 
90  // Easy special cases
91  // FIXME: G4Pow doesn't check whether the argument gets too large, which is unsafe! Max is 512.
92  if(m==l) return (l%2 ? -1. : 1.) *
93  G4Exp(g4pow->logfactorial(2*l) - g4pow->logfactorial(l)) *
94  G4Exp(G4Log((1.-x*x)*0.25)*0.5*G4double(l));
95  if(m==l-1) return x*(2.*G4double(m)+1.)*EvalAssocLegendrePoly(m,m,x);
96 
97  // Otherwise calculate recursively
98  return (x*G4double(2*l-1)*EvalAssocLegendrePoly(l-1,m,x) -
99  (G4double(l+m-1))*EvalAssocLegendrePoly(l-2,m,x))/G4double(l-m);
100 }
101 
103 {
104  if(orderMax > 30) {
105  G4cout << "G4LegendrePolynomial::GetCoefficient(): "
106  << "I refuse to make a Legendre Polynomial of order "
107  << orderMax << G4endl;
108  return;
109  }
110  while(fCoefficients.size() < orderMax+1) { /* Loop checking, 30-Oct-2015, G.Folger */
111  size_t order = fCoefficients.size();
112  fCoefficients.resize(order+1);
113  if(order <= 1) fCoefficients[order].push_back(1.);
114  else {
115  for(size_t iCoeff = 0; iCoeff < order+1; ++iCoeff) {
116  if((order % 2) == (iCoeff % 2)) {
117  G4double coeff = 0;
118  if(iCoeff <= order-2) coeff -= fCoefficients[order-2][iCoeff/2]*G4double(order-1);
119  if(iCoeff > 0) coeff += fCoefficients[order-1][(iCoeff-1)/2]*G4double(2*order-1);
120  coeff /= G4double(order);
121  fCoefficients[order].push_back(coeff);
122  }
123  }
124  }
125  }
126 }
127 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double EvalLegendrePoly(G4int order, G4double x)
Definition: G4Pow.hh:56
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double factorial(G4int Z) const
Definition: G4Pow.hh:264
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double EvalAssocLegendrePoly(G4int l, G4int m, G4double x)
G4double logfactorial(G4int Z) const
Definition: G4Pow.hh:269
const G4double x[NPOINTSGL]
#define G4endl
Definition: G4ios.hh:61
static const double m
Definition: G4SIunits.hh:128
double G4double
Definition: G4Types.hh:76
void BuildUpToOrder(size_t order)
G4double GetCoefficient(size_t i, size_t order)