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