| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 25
26 #include "G4ios.hh" 26 #include "G4ios.hh"
27 #include "G4LegendrePolynomial.hh" 27 #include "G4LegendrePolynomial.hh"
28 #include "G4Pow.hh" 28 #include "G4Pow.hh"
29 #include "G4Exp.hh" 29 #include "G4Exp.hh"
30 #include "G4Log.hh" 30 #include "G4Log.hh"
31 31
32 using namespace std; 32 using namespace std;
33 33
34 G4double G4LegendrePolynomial::GetCoefficient( << 34 G4double G4LegendrePolynomial::GetCoefficient(size_t i, size_t order)
35 { 35 {
36 if(order >= fCoefficients.size()) BuildUpToO 36 if(order >= fCoefficients.size()) BuildUpToOrder(order);
37 if(order >= fCoefficients.size() || 37 if(order >= fCoefficients.size() ||
38 i/2 >= fCoefficients[order].size() || 38 i/2 >= fCoefficients[order].size() ||
39 (i%2) != order %2) return 0; 39 (i%2) != order %2) return 0;
40 return fCoefficients[order][i/2]; 40 return fCoefficients[order][i/2];
41 } 41 }
42 42
43 G4double G4LegendrePolynomial::EvalLegendrePol 43 G4double G4LegendrePolynomial::EvalLegendrePoly(G4int order, G4double x)
44 { 44 {
45 // Call EvalAssocLegendrePoly with m=0 45 // Call EvalAssocLegendrePoly with m=0
46 return (EvalAssocLegendrePoly(order,0,x)); 46 return (EvalAssocLegendrePoly(order,0,x));
47 } 47 }
48 48
49 G4double G4LegendrePolynomial::EvalAssocLegend << 49 G4double G4LegendrePolynomial::EvalAssocLegendrePoly(G4int l, G4int m, G4double x,
>> 50 map<G4int, map<G4int, G4double> >* cache)
50 { 51 {
51 // Invalid calls --> 0. (Keeping for backwar << 52 // Calculate P_l^m(x).
52 if (l<0 || m<-l || m>l) return 0.0; << 53 // If cache ptr is non-null, use cache[l][m] if it exists, otherwise compute
>> 54 // P_l^m(x) and cache it in that position. The cache speeds up calculations
>> 55 // where many P_l^m computations are need at the same value of x.
>> 56
>> 57 if(l<0 || m<-l || m>l) return 0;
>> 58 G4Pow* g4pow = G4Pow::GetInstance();
>> 59
>> 60 // Use non-log factorial for low l, m: it is more efficient until
>> 61 // l and m get above 10 or so.
>> 62 // FIXME: G4Pow doesn't check whether the argument gets too large,
>> 63 // which is unsafe! Max is 512; VI: It is assume that Geant4 does not
>> 64 // need higher order
>> 65 if(m<0) {
>> 66 G4double value = (m%2 ? -1. : 1.) * EvalAssocLegendrePoly(l, -m, x);
>> 67 if(l < 10) return value * g4pow->factorial(l+m)/g4pow->factorial(l-m);
>> 68 else { return value * G4Exp(g4pow->logfactorial(l+m) - g4pow->logfactorial(l-m));
>> 69 }
>> 70 }
53 71
54 // New check: g4pow doesn't handle integer a << 72 // hard-code the first few orders for speed
55 // For us that means: << 73 if(l==0) return 1;
56 if ((l+m) > 512 || (l-m) > 512 || (2*m) > 51 << 74 if(l==1) {
>> 75 if(m==0){return x;}
>> 76 /*m==1*/ return -sqrt(1.-x*x);
>> 77 }
>> 78 if(l<5) {
>> 79 G4double x2 = x*x;
>> 80 if(l==2) {
>> 81 if(m==0){return 0.5*(3.*x2 - 1.);}
>> 82 if(m==1){return -3.*x*sqrt(1.-x2);}
>> 83 /*m==2*/ return 3.*(1.-x2);
>> 84 }
>> 85 if(l==3) {
>> 86 if(m==0){return 0.5*(5.*x*x2 - 3.*x);}
>> 87 if(m==1){return -1.5*(5.*x2-1.)*sqrt(1.-x2);}
>> 88 if(m==2){return 15.*x*(1.-x2);}
>> 89 /*m==3*/ return -15.*(1.-x2)*sqrt(1.-x2);
>> 90 }
>> 91 if(l==4) {
>> 92 if(m==0){return 0.125*(35.*x2*x2 - 30.*x2 + 3.);}
>> 93 if(m==1){return -2.5*(7.*x*x2-3.*x)*sqrt(1.-x2);}
>> 94 if(m==2){return 7.5*(7.*x2-1.)*(1.-x2);}
>> 95 if(m==3){return -105.*x*(1.-x2)*sqrt(1.-x2);}
>> 96 /*m==4*/ return 105.*(1. - 2.*x2 + x2*x2);
>> 97 }
>> 98 }
57 99
58 G4Pow* g4pow = G4Pow::GetInstance(); << 100 // Easy special cases
59 G4double x2 = x*x; << 101 // FIXME: G4Pow doesn't check whether the argument gets too large, which is unsafe! Max is 512.
>> 102 if(m==l) return (l%2 ? -1. : 1.) *
>> 103 G4Exp(g4pow->logfactorial(2*l) - g4pow->logfactorial(l)) *
>> 104 G4Exp(G4Log((1.-x*x)*0.25)*0.5*G4double(l));
>> 105 if(m==l-1) return x*(2.*G4double(m)+1.)*EvalAssocLegendrePoly(m,m,x);
>> 106
>> 107 // See if we have this value cached.
>> 108 if(cache != NULL && cache->count(l) > 0 && (*cache)[l].count(m) > 0) {
>> 109 return (*cache)[l][m];
>> 110 }
60 111
61 // hard-code the first few orders for speed << 112 // Otherwise calculate recursively
62 switch (l) { << 113 G4double value = (x*G4double(2*l-1)*EvalAssocLegendrePoly(l-1,m,x) -
63 case 0 : << 114 (G4double(l+m-1))*EvalAssocLegendrePoly(l-2,m,x))/G4double(l-m);
64 return 1; << 115
65 case 1 : << 116 // If we are working with a cache, cache this value.
66 switch (m) { << 117 if(cache != NULL) {
67 case -1 : return 0.5 * sqrt(1.-x2); << 118 (*cache)[l][m] = value;
68 case 0 : return x; <<
69 case 1 : return -sqrt(1.-x2); <<
70 }; <<
71 break; <<
72 case 2 : <<
73 switch (m) { <<
74 case -2 : return 0.125 * (1.0 - x2); <<
75 case -1 : return 0.5 * x * sqrt(1.0 - <<
76 case 0 : return 0.5*(3.*x2 - 1.); <<
77 case 1 : return -3.*x*sqrt(1.-x2); <<
78 case 2 : return 3.*(1.-x2); <<
79 }; <<
80 break; <<
81 case 3 : <<
82 switch (m) { <<
83 case -3 : return (1.0/48.0) * (1.0 - x <<
84 case -2 : return 0.125 * x * (1.0 - x2 <<
85 case -1 : return 0.125 * (5.0 * x2 - 1 <<
86 case 0 : return 0.5*(5.*x*x2 - 3.*x); <<
87 case 1 : return -1.5*(5.*x2-1.)*sqrt(1 <<
88 case 2 : return 15.*x*(1.-x2); <<
89 case 3 : return -15.*(1.-x2)*sqrt(1.-x <<
90 }; <<
91 break; <<
92 case 4 : <<
93 switch (m) { <<
94 case -4 : return (105.0/40320.0)*(1. - <<
95 case -3 : return (105.0/5040.0)*x*(1.- <<
96 case -2 : return (15.0/720.0)*(7.*x2-1 <<
97 case -1 : return 0.125*(7.*x*x2-3.*x)* <<
98 case 0 : return 0.125*(35.*x2*x2 - 30. <<
99 case 1 : return -2.5*(7.*x*x2-3.*x)*sq <<
100 case 2 : return 7.5*(7.*x2-1.)*(1.-x2) <<
101 case 3 : return -105.*x*(1.-x2)*sqrt(1 <<
102 case 4 : return 105.*(1. - 2.*x2 + x2* <<
103 }; <<
104 break; <<
105 }; <<
106 <<
107 // if m<0, compute P[l,-m,x] and correct <<
108 if (m < 0) <<
109 { <<
110 G4double complementary_value = EvalAssocLe <<
111 return complementary_value * (m%2==0 ? 1.0 <<
112 } <<
113 <<
114 // Iteratively walk up from P[m,m,x] to P[l, <<
115 <<
116 // prime the pump: P[l<m,m,x] = 0 <<
117 G4double previous = 0.0; <<
118 <<
119 // prime the pump: P[m,m,x] <<
120 G4double current; <<
121 if (m == 0) current = 1.0; <<
122 else if (m == 1) current = -sqrt((1.0 - (x2) <<
123 else { <<
124 current = (m%2==0 ? 1.0 : -1.0) * <<
125 G4Exp(g4pow->logfactorial(2*m) - g4pow- <<
126 G4Exp(G4Log((1.0-(x2))*0.25)*0.5*G4doub <<
127 } <<
128 <<
129 // Work up to P[l,m,x] <<
130 for(G4int i=m+1; i<=l; i++) <<
131 { <<
132 G4double next = (-(G4double(i+m-1))*previo <<
133 previous = current; <<
134 current = next; <<
135 } 119 }
136 << 120 return value;
137 return current; <<
138 } 121 }
139 122
140 void G4LegendrePolynomial::BuildUpToOrder(std: << 123 void G4LegendrePolynomial::BuildUpToOrder(size_t orderMax)
141 { 124 {
142 if(orderMax > 30) { 125 if(orderMax > 30) {
143 G4cout << "G4LegendrePolynomial::GetCoeffi 126 G4cout << "G4LegendrePolynomial::GetCoefficient(): "
144 << "I refuse to make a Legendre Pol 127 << "I refuse to make a Legendre Polynomial of order "
145 << orderMax << G4endl; 128 << orderMax << G4endl;
146 return; 129 return;
147 } 130 }
148 while(fCoefficients.size() < orderMax+1) { 131 while(fCoefficients.size() < orderMax+1) { /* Loop checking, 30-Oct-2015, G.Folger */
149 std::size_t order = fCoefficients.size(); << 132 size_t order = fCoefficients.size();
150 fCoefficients.resize(order+1); 133 fCoefficients.resize(order+1);
151 if(order <= 1) fCoefficients[order].push_b 134 if(order <= 1) fCoefficients[order].push_back(1.);
152 else { 135 else {
153 for(std::size_t iCoeff = 0; iCoeff < ord << 136 for(size_t iCoeff = 0; iCoeff < order+1; ++iCoeff) {
154 if((order % 2) == (iCoeff % 2)) { 137 if((order % 2) == (iCoeff % 2)) {
155 G4double coeff = 0; 138 G4double coeff = 0;
156 if(iCoeff <= order-2) coeff -= fCoef 139 if(iCoeff <= order-2) coeff -= fCoefficients[order-2][iCoeff/2]*G4double(order-1);
157 if(iCoeff > 0) coeff += fCoefficient 140 if(iCoeff > 0) coeff += fCoefficients[order-1][(iCoeff-1)/2]*G4double(2*order-1);
158 coeff /= G4double(order); 141 coeff /= G4double(order);
159 fCoefficients[order].push_back(coeff 142 fCoefficients[order].push_back(coeff);
160 } 143 }
161 } 144 }
162 } 145 }
163 } 146 }
164 } 147 }
165 148
166 149