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