| 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 // G4ChebyshevApproximation class implementati <<
27 // 26 //
28 // Author: V.Grichine, 24.04.97 << 27 //
29 // ------------------------------------------- <<
30 28
31 #include "G4ChebyshevApproximation.hh" 29 #include "G4ChebyshevApproximation.hh"
32 #include "G4PhysicalConstants.hh" 30 #include "G4PhysicalConstants.hh"
33 31
34 // Constructor for initialisation of the class 32 // Constructor for initialisation of the class data members.
35 // It creates the array fChebyshevCof[0,...,fN 33 // It creates the array fChebyshevCof[0,...,fNumber-1], fNumber = n ;
36 // which consists of Chebyshev coefficients de 34 // which consists of Chebyshev coefficients describing the function
37 // pointed by pFunction. The values a and b fi 35 // pointed by pFunction. The values a and b fix the interval of validity
38 // of the Chebyshev approximation. << 36 // of the Chebyshev approximation.
39 37
40 G4ChebyshevApproximation::G4ChebyshevApproxima << 38 G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
41 << 39 G4int n,
42 : fFunction(pFunction) << 40 G4double a,
43 , fNumber(n) << 41 G4double b )
44 , fChebyshevCof(new G4double[fNumber]) << 42 : fFunction(pFunction), fNumber(n),
45 , fMean(0.5 * (b + a)) << 43 fChebyshevCof(new G4double[fNumber]),
46 , fDiff(0.5 * (b - a)) << 44 fMean(0.5*(b+a)), fDiff(0.5*(b-a))
47 { << 45 {
48 G4int i = 0, j = 0; << 46 G4int i=0, j=0 ;
49 G4double rootSum = 0.0, cofj = 0.0; << 47 G4double rootSum=0.0, cofj=0.0 ;
50 auto tempFunction = new G4double[fNumber]; << 48 G4double* tempFunction = new G4double[fNumber] ;
51 G4double weight = 2.0 / fNumber; << 49 G4double weight = 2.0/fNumber ;
52 G4double cof = 0.5 * weight * pi; << 50 G4double cof = 0.5*weight*pi ; // pi/n
53 << 51
54 for(i = 0; i < fNumber; ++i) << 52 for (i=0;i<fNumber;i++)
55 { << 53 {
56 rootSum = std::cos(cof * (i + 0.5) << 54 rootSum = std::cos(cof*(i+0.5)) ;
57 tempFunction[i] = fFunction(rootSum * fDif << 55 tempFunction[i]= fFunction(rootSum*fDiff+fMean) ;
58 } << 56 }
59 for(j = 0; j < fNumber; ++j) << 57 for (j=0;j<fNumber;j++)
60 { << 58 {
61 cofj = cof * j; << 59 cofj = cof*j ;
62 rootSum = 0.0; << 60 rootSum = 0.0 ;
63 << 61
64 for(i = 0; i < fNumber; ++i) << 62 for (i=0;i<fNumber;i++)
65 { << 63 {
66 rootSum += tempFunction[i] * std::cos(co << 64 rootSum += tempFunction[i]*std::cos(cofj*(i+0.5)) ;
67 } << 65 }
68 fChebyshevCof[j] = weight * rootSum; << 66 fChebyshevCof[j] = weight*rootSum ;
69 } << 67 }
70 delete[] tempFunction; << 68 delete[] tempFunction ;
71 } 69 }
72 70
73 // ------------------------------------------- 71 // --------------------------------------------------------------------
74 // 72 //
75 // Constructor for creation of Chebyshev coeff 73 // Constructor for creation of Chebyshev coefficients for mx-derivative
76 // from pFunction. The value of mx ! MUST BE ! 74 // from pFunction. The value of mx ! MUST BE ! < nx , because the result
77 // array of fChebyshevCof will be of (nx-mx) s 75 // array of fChebyshevCof will be of (nx-mx) size. The values a and b
78 // fix the interval of validity of the Chebysh << 76 // fix the interval of validity of the Chebyshev approximation.
79 77
80 G4ChebyshevApproximation::G4ChebyshevApproxima << 78 G4ChebyshevApproximation::
81 << 79 G4ChebyshevApproximation( function pFunction,
82 << 80 G4int nx, G4int mx,
83 : fFunction(pFunction) << 81 G4double a, G4double b )
84 , fNumber(nx) << 82 : fFunction(pFunction), fNumber(nx),
85 , fChebyshevCof(new G4double[fNumber]) << 83 fChebyshevCof(new G4double[fNumber]),
86 , fMean(0.5 * (b + a)) << 84 fMean(0.5*(b+a)), fDiff(0.5*(b-a))
87 , fDiff(0.5 * (b - a)) << 85 {
88 { << 86 if(nx <= mx)
89 if(nx <= mx) << 87 {
90 { << 88 G4Exception("G4ChebyshevApproximation::G4ChebyshevApproximation()",
91 G4Exception("G4ChebyshevApproximation::G4C << 89 "InvalidCall", FatalException, "Invalid arguments !") ;
92 "InvalidCall", FatalException, << 90 }
93 } << 91 G4int i=0, j=0 ;
94 G4int i = 0, j = 0; << 92 G4double rootSum = 0.0, cofj=0.0;
95 G4double rootSum = 0.0, cofj = 0.0; << 93 G4double* tempFunction = new G4double[fNumber] ;
96 auto tempFunction = new G4double[fNumber]; << 94 G4double weight = 2.0/fNumber ;
97 G4double weight = 2.0 / fNumber; << 95 G4double cof = 0.5*weight*pi ; // pi/nx
98 G4double cof = 0.5 * weight * pi; << 96
99 << 97 for (i=0;i<fNumber;i++)
100 for(i = 0; i < fNumber; ++i) << 98 {
101 { << 99 rootSum = std::cos(cof*(i+0.5)) ;
102 rootSum = std::cos(cof * (i + 0.5) << 100 tempFunction[i] = fFunction(rootSum*fDiff+fMean) ;
103 tempFunction[i] = fFunction(rootSum * fDif << 101 }
104 } << 102 for (j=0;j<fNumber;j++)
105 for(j = 0; j < fNumber; ++j) << 103 {
106 { << 104 cofj = cof*j ;
107 cofj = cof * j; << 105 rootSum = 0.0 ;
108 rootSum = 0.0; << 106
109 << 107 for (i=0;i<fNumber;i++)
110 for(i = 0; i < fNumber; ++i) << 108 {
111 { << 109 rootSum += tempFunction[i]*std::cos(cofj*(i+0.5)) ;
112 rootSum += tempFunction[i] * std::cos(co << 110 }
113 } << 111 fChebyshevCof[j] = weight*rootSum ; // corresponds to pFunction
114 fChebyshevCof[j] = weight * rootSum; // c << 112 }
115 } << 113 // Chebyshev coefficients for (mx)-derivative of pFunction
116 // Chebyshev coefficients for (mx)-derivativ << 114
117 << 115 for(i=1;i<=mx;i++)
118 for(i = 1; i <= mx; ++i) << 116 {
119 { << 117 DerivativeChebyshevCof(tempFunction) ;
120 DerivativeChebyshevCof(tempFunction); << 118 fNumber-- ;
121 fNumber--; << 119 for(j=0;j<fNumber;j++)
122 for(j = 0; j < fNumber; ++j) << 120 {
123 { << 121 fChebyshevCof[j] = tempFunction[j] ; // corresponds to (i)-derivative
124 fChebyshevCof[j] = tempFunction[j]; // << 122 }
125 } << 123 }
126 } << 124 delete[] tempFunction ; // delete of dynamically allocated tempFunction
127 delete[] tempFunction; // delete of dynamic <<
128 } 125 }
129 126
130 // ------------------------------------------- 127 // ------------------------------------------------------
131 // 128 //
132 // Constructor for creation of Chebyshev coeff 129 // Constructor for creation of Chebyshev coefficients for integral
133 // from pFunction. 130 // from pFunction.
134 131
135 G4ChebyshevApproximation::G4ChebyshevApproxima << 132 G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
136 << 133 G4double a,
137 << 134 G4double b,
138 : fFunction(pFunction) << 135 G4int n )
139 , fNumber(n) << 136 : fFunction(pFunction), fNumber(n),
140 , fChebyshevCof(new G4double[fNumber]) << 137 fChebyshevCof(new G4double[fNumber]),
141 , fMean(0.5 * (b + a)) << 138 fMean(0.5*(b+a)), fDiff(0.5*(b-a))
142 , fDiff(0.5 * (b - a)) << 139 {
143 { << 140 G4int i=0, j=0;
144 G4int i = 0, j = 0; << 141 G4double rootSum=0.0, cofj=0.0;
145 G4double rootSum = 0.0, cofj = 0.0; << 142 G4double* tempFunction = new G4double[fNumber] ;
146 auto tempFunction = new G4double[fNumber]; << 143 G4double weight = 2.0/fNumber;
147 G4double weight = 2.0 / fNumber; << 144 G4double cof = 0.5*weight*pi ; // pi/n
148 G4double cof = 0.5 * weight * pi; << 145
149 << 146 for (i=0;i<fNumber;i++)
150 for(i = 0; i < fNumber; ++i) << 147 {
151 { << 148 rootSum = std::cos(cof*(i+0.5)) ;
152 rootSum = std::cos(cof * (i + 0.5) << 149 tempFunction[i]= fFunction(rootSum*fDiff+fMean) ;
153 tempFunction[i] = fFunction(rootSum * fDif << 150 }
154 } << 151 for (j=0;j<fNumber;j++)
155 for(j = 0; j < fNumber; ++j) << 152 {
156 { << 153 cofj = cof*j ;
157 cofj = cof * j; << 154 rootSum = 0.0 ;
158 rootSum = 0.0; << 155
159 << 156 for (i=0;i<fNumber;i++)
160 for(i = 0; i < fNumber; ++i) << 157 {
161 { << 158 rootSum += tempFunction[i]*std::cos(cofj*(i+0.5)) ;
162 rootSum += tempFunction[i] * std::cos(co << 159 }
163 } << 160 fChebyshevCof[j] = weight*rootSum ; // corresponds to pFunction
164 fChebyshevCof[j] = weight * rootSum; // c << 161 }
165 } << 162 // Chebyshev coefficients for integral of pFunction
166 // Chebyshev coefficients for integral of pF << 163
167 << 164 IntegralChebyshevCof(tempFunction) ;
168 IntegralChebyshevCof(tempFunction); << 165 for(j=0;j<fNumber;j++)
169 for(j = 0; j < fNumber; ++j) << 166 {
170 { << 167 fChebyshevCof[j] = tempFunction[j] ; // corresponds to integral
171 fChebyshevCof[j] = tempFunction[j]; // co << 168 }
172 } << 169 delete[] tempFunction ; // delete of dynamically allocated tempFunction
173 delete[] tempFunction; // delete of dynamic <<
174 } 170 }
175 171
>> 172
>> 173
176 // ------------------------------------------- 174 // ---------------------------------------------------------------
177 // 175 //
178 // Destructor deletes the array of Chebyshev c 176 // Destructor deletes the array of Chebyshev coefficients
179 177
180 G4ChebyshevApproximation::~G4ChebyshevApproxim 178 G4ChebyshevApproximation::~G4ChebyshevApproximation()
181 { 179 {
182 delete[] fChebyshevCof; << 180 delete[] fChebyshevCof ;
183 } 181 }
184 182
185 // ------------------------------------------- 183 // ---------------------------------------------------------------
186 // 184 //
187 // Access function for Chebyshev coefficients 185 // Access function for Chebyshev coefficients
188 // 186 //
189 187
190 G4double G4ChebyshevApproximation::GetChebyshe << 188
>> 189 G4double
>> 190 G4ChebyshevApproximation::GetChebyshevCof(G4int number) const
191 { 191 {
192 if(number < 0 && number >= fNumber) << 192 if(number < 0 && number >= fNumber)
193 { << 193 {
194 G4Exception("G4ChebyshevApproximation::Get << 194 G4Exception("G4ChebyshevApproximation::GetChebyshevCof()",
195 FatalException, "Argument out << 195 "InvalidCall", FatalException, "Argument out of range !") ;
196 } << 196 }
197 return fChebyshevCof[number]; << 197 return fChebyshevCof[number] ;
198 } 198 }
199 199
200 // ------------------------------------------- 200 // --------------------------------------------------------------
201 // 201 //
202 // Evaluate the value of fFunction at the poin 202 // Evaluate the value of fFunction at the point x via the Chebyshev coefficients
203 // fChebyshevCof[0,...,fNumber-1] 203 // fChebyshevCof[0,...,fNumber-1]
204 204
205 G4double G4ChebyshevApproximation::ChebyshevEv << 205 G4double
>> 206 G4ChebyshevApproximation::ChebyshevEvaluation(G4double x) const
206 { 207 {
207 G4double evaluate = 0.0, evaluate2 = 0.0, te << 208 G4double evaluate = 0.0, evaluate2 = 0.0, temp = 0.0,
208 xReduced2 = 0.0; << 209 xReduced = 0.0, xReduced2 = 0.0 ;
209 210
210 if((x - fMean + fDiff) * (x - fMean - fDiff) << 211 if ((x-fMean+fDiff)*(x-fMean-fDiff) > 0.0)
211 { << 212 {
212 G4Exception("G4ChebyshevApproximation::Che << 213 G4Exception("G4ChebyshevApproximation::ChebyshevEvaluation()",
213 "InvalidCall", FatalException, << 214 "InvalidCall", FatalException, "Invalid argument !") ;
214 } << 215 }
215 xReduced = (x - fMean) / fDiff; << 216 xReduced = (x-fMean)/fDiff ;
216 xReduced2 = 2.0 * xReduced; << 217 xReduced2 = 2.0*xReduced ;
217 for(G4int i = fNumber - 1; i >= 1; --i) << 218 for (G4int i=fNumber-1;i>=1;i--)
218 { << 219 {
219 temp = evaluate; << 220 temp = evaluate ;
220 evaluate = xReduced2 * evaluate - evaluat << 221 evaluate = xReduced2*evaluate - evaluate2 + fChebyshevCof[i] ;
221 evaluate2 = temp; << 222 evaluate2 = temp ;
222 } << 223 }
223 return xReduced * evaluate - evaluate2 + 0.5 << 224 return xReduced*evaluate - evaluate2 + 0.5*fChebyshevCof[0] ;
224 } 225 }
225 226
226 // ------------------------------------------- 227 // ------------------------------------------------------------------
227 // 228 //
228 // Returns the array derCof[0,...,fNumber-2], << 229 // Returns the array derCof[0,...,fNumber-2], the Chebyshev coefficients of the
229 // derivative of the function whose coefficien 230 // derivative of the function whose coefficients are fChebyshevCof
230 231
231 void G4ChebyshevApproximation::DerivativeCheby << 232 void
>> 233 G4ChebyshevApproximation::DerivativeChebyshevCof(G4double derCof[]) const
232 { 234 {
233 G4double cof = 1.0 / fDiff; << 235 G4double cof = 1.0/fDiff ;
234 derCof[fNumber - 1] = 0.0; << 236 derCof[fNumber-1] = 0.0 ;
235 derCof[fNumber - 2] = 2 * (fNumber - 1) * fC << 237 derCof[fNumber-2] = 2*(fNumber-1)*fChebyshevCof[fNumber-1] ;
236 for(G4int i = fNumber - 3; i >= 0; --i) << 238 for(G4int i=fNumber-3;i>=0;i--)
237 { << 239 {
238 derCof[i] = derCof[i + 2] + 2 * (i + 1) * << 240 derCof[i] = derCof[i+2] + 2*(i+1)*fChebyshevCof[i+1] ;
239 } << 241 }
240 for(G4int j = 0; j < fNumber; ++j) << 242 for(G4int j=0;j<fNumber;j++)
241 { << 243 {
242 derCof[j] *= cof; << 244 derCof[j] *= cof ;
243 } << 245 }
244 } 246 }
245 247
246 // ------------------------------------------- 248 // ------------------------------------------------------------------------
247 // 249 //
248 // This function produces the array integralCo 250 // This function produces the array integralCof[0,...,fNumber-1] , the Chebyshev
249 // coefficients of the integral of the functio << 251 // coefficients of the integral of the function whose coefficients are
250 // fChebyshevCof[]. The constant of integratio << 252 // fChebyshevCof[]. The constant of integration is set so that the integral
251 // vanishes at the point (fMean - fDiff), i.e. << 253 // vanishes at the point (fMean - fDiff), i.e. at the begining of the interval of
252 // of validity (we start the integration from << 254 // validity (we start the integration from this point).
253 // << 255 //
254 << 256
255 void G4ChebyshevApproximation::IntegralChebysh << 257 void
256 G4double integralCof[]) const << 258 G4ChebyshevApproximation::IntegralChebyshevCof(G4double integralCof[]) const
257 { << 259 {
258 G4double cof = 0.5 * fDiff, sum = 0.0, facto << 260 G4double cof = 0.5*fDiff, sum = 0.0, factor = 1.0 ;
259 for(G4int i = 1; i < fNumber - 1; ++i) << 261 for(G4int i=1;i<fNumber-1;i++)
260 { << 262 {
261 integralCof[i] = cof * (fChebyshevCof[i - << 263 integralCof[i] = cof*(fChebyshevCof[i-1] - fChebyshevCof[i+1])/i ;
262 sum += factor * integralCof[i]; << 264 sum += factor*integralCof[i] ;
263 factor = -factor; << 265 factor = -factor ;
264 } << 266 }
265 integralCof[fNumber - 1] = cof * fChebyshevC << 267 integralCof[fNumber-1] = cof*fChebyshevCof[fNumber-2]/(fNumber-1) ;
266 sum += factor * integralCof[fNumber - 1]; << 268 sum += factor*integralCof[fNumber-1] ;
267 integralCof[0] = 2.0 * sum; // set the cons << 269 integralCof[0] = 2.0*sum ; // set the constant of integration
268 } << 270 }
269 271