| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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
34 G4double G4LegendrePolynomial::GetCoefficient(
35 {
36 if(order >= fCoefficients.size()) BuildUpToO
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
43 G4double G4LegendrePolynomial::EvalLegendrePol
44 {
45 // Call EvalAssocLegendrePoly with m=0
46 return (EvalAssocLegendrePoly(order,0,x));
47 }
48
49 G4double G4LegendrePolynomial::EvalAssocLegend
50 {
51 // Invalid calls --> 0. (Keeping for backwar
52 if (l<0 || m<-l || m>l) return 0.0;
53
54 // New check: g4pow doesn't handle integer a
55 // For us that means:
56 if ((l+m) > 512 || (l-m) > 512 || (2*m) > 51
57
58 G4Pow* g4pow = G4Pow::GetInstance();
59 G4double x2 = x*x;
60
61 // hard-code the first few orders for speed
62 switch (l) {
63 case 0 :
64 return 1;
65 case 1 :
66 switch (m) {
67 case -1 : return 0.5 * sqrt(1.-x2);
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 }
136
137 return current;
138 }
139
140 void G4LegendrePolynomial::BuildUpToOrder(std:
141 {
142 if(orderMax > 30) {
143 G4cout << "G4LegendrePolynomial::GetCoeffi
144 << "I refuse to make a Legendre Pol
145 << orderMax << G4endl;
146 return;
147 }
148 while(fCoefficients.size() < orderMax+1) {
149 std::size_t order = fCoefficients.size();
150 fCoefficients.resize(order+1);
151 if(order <= 1) fCoefficients[order].push_b
152 else {
153 for(std::size_t iCoeff = 0; iCoeff < ord
154 if((order % 2) == (iCoeff % 2)) {
155 G4double coeff = 0;
156 if(iCoeff <= order-2) coeff -= fCoef
157 if(iCoeff > 0) coeff += fCoefficient
158 coeff /= G4double(order);
159 fCoefficients[order].push_back(coeff
160 }
161 }
162 }
163 }
164 }
165
166