| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4ChebyshevApproximation class implementati <<
27 // 23 //
28 // Author: V.Grichine, 24.04.97 << 24 // $Id: G4ChebyshevApproximation.cc,v 1.3 2001/07/11 10:00:41 gunter Exp $
29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-04-00 $
>> 26 //
30 27
31 #include "G4ChebyshevApproximation.hh" 28 #include "G4ChebyshevApproximation.hh"
32 #include "G4PhysicalConstants.hh" <<
33 29
34 // Constructor for initialisation of the class << 30
35 // It creates the array fChebyshevCof[0,...,fN << 31 // Constructor for initialisation of the class data members. It creates the array
36 // which consists of Chebyshev coefficients de << 32 // fChebyshevCof[0,...,fNumber-1], fNumber = n ; which consists of Chebyshev
37 // pointed by pFunction. The values a and b fi << 33 // coefficients describing the function pointed by pFunction. The values a and b
38 // of the Chebyshev approximation. << 34 // fixe the interval of validity of Chebyshev approximation.
39 << 35
40 G4ChebyshevApproximation::G4ChebyshevApproxima << 36
41 << 37 G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
42 : fFunction(pFunction) << 38 G4int n,
43 , fNumber(n) << 39 G4double a,
44 , fChebyshevCof(new G4double[fNumber]) << 40 G4double b )
45 , fMean(0.5 * (b + a)) <<
46 , fDiff(0.5 * (b - a)) <<
47 { 41 {
48 G4int i = 0, j = 0; << 42 G4int i, j ;
49 G4double rootSum = 0.0, cofj = 0.0; << 43 G4double rootSum, cof, cofj, weight ;
50 auto tempFunction = new G4double[fNumber]; << 44
51 G4double weight = 2.0 / fNumber; << 45 fFunction = pFunction ;
52 G4double cof = 0.5 * weight * pi; << 46 fNumber = n ;
53 << 47 fDiff = 0.5*(b-a) ;
54 for(i = 0; i < fNumber; ++i) << 48 fMean = 0.5*(b+a) ;
55 { << 49 fChebyshevCof = new G4double[fNumber] ;
56 rootSum = std::cos(cof * (i + 0.5) << 50
57 tempFunction[i] = fFunction(rootSum * fDif << 51 G4double* tempFunction = new G4double[fNumber] ;
58 } << 52
59 for(j = 0; j < fNumber; ++j) << 53 weight = 2.0/fNumber ;
60 { << 54 cof = 0.5*weight*pi ; // pi/n
61 cofj = cof * j; << 55
62 rootSum = 0.0; << 56 for (i=0;i<fNumber;i++)
63 << 57 {
64 for(i = 0; i < fNumber; ++i) << 58 rootSum = cos(cof*(i+0.5)) ;
65 { << 59 tempFunction[i]= fFunction(rootSum*fDiff+fMean) ;
66 rootSum += tempFunction[i] * std::cos(co << 60 }
67 } << 61 for (j=0;j<fNumber;j++)
68 fChebyshevCof[j] = weight * rootSum; << 62 {
69 } << 63 cofj = cof*j ;
70 delete[] tempFunction; << 64 rootSum = 0.0 ;
>> 65
>> 66 for (i=0;i<fNumber;i++)
>> 67 {
>> 68 rootSum += tempFunction[i]*cos(cofj*(i+0.5)) ;
>> 69 }
>> 70 fChebyshevCof[j] = weight*rootSum ;
>> 71 }
>> 72 delete[] tempFunction ;
71 } 73 }
72 74
73 // ------------------------------------------- 75 // --------------------------------------------------------------------
74 // 76 //
75 // Constructor for creation of Chebyshev coeff << 77 // Constructor for creation of Chebyshev coefficients for m-derivative
76 // from pFunction. The value of mx ! MUST BE ! << 78 // from pFunction. The value of m ! MUST BE ! < n , because the result
77 // array of fChebyshevCof will be of (nx-mx) s << 79 // array of fChebyshevCof will be of (n-m) size. The values a and b
78 // fix the interval of validity of the Chebysh << 80 // fixe the interval of validity of Chebyshev approximation.
79 << 81
80 G4ChebyshevApproximation::G4ChebyshevApproxima << 82
81 << 83 G4ChebyshevApproximation::
82 << 84 G4ChebyshevApproximation( function pFunction,
83 : fFunction(pFunction) << 85 G4int n,
84 , fNumber(nx) << 86 G4int m,
85 , fChebyshevCof(new G4double[fNumber]) << 87 G4double a,
86 , fMean(0.5 * (b + a)) << 88 G4double b )
87 , fDiff(0.5 * (b - a)) <<
88 { 89 {
89 if(nx <= mx) << 90 if(n <= m)
90 { << 91 {
91 G4Exception("G4ChebyshevApproximation::G4C << 92 G4Exception
92 "InvalidCall", FatalException, << 93 ("Invalid arguments in G4ChebyshevApproximation::G4ChebyshevApproximation") ;
93 } << 94 }
94 G4int i = 0, j = 0; << 95 G4int i, j ;
95 G4double rootSum = 0.0, cofj = 0.0; << 96 G4double rootSum, cof, cofj, weight ;
96 auto tempFunction = new G4double[fNumber]; << 97
97 G4double weight = 2.0 / fNumber; << 98 fFunction = pFunction ;
98 G4double cof = 0.5 * weight * pi; << 99 fNumber = n ;
99 << 100 fDiff = 0.5*(b-a) ;
100 for(i = 0; i < fNumber; ++i) << 101 fMean = 0.5*(b+a) ;
101 { << 102 fChebyshevCof = new G4double[fNumber] ;
102 rootSum = std::cos(cof * (i + 0.5) << 103
103 tempFunction[i] = fFunction(rootSum * fDif << 104 G4double* tempFunction = new G4double[fNumber] ;
104 } << 105
105 for(j = 0; j < fNumber; ++j) << 106 weight = 2.0/fNumber ;
106 { << 107 cof = 0.5*weight*pi ; // pi/n
107 cofj = cof * j; << 108
108 rootSum = 0.0; << 109 for (i=0;i<fNumber;i++)
109 << 110 {
110 for(i = 0; i < fNumber; ++i) << 111 rootSum = cos(cof*(i+0.5)) ;
111 { << 112 tempFunction[i] = fFunction(rootSum*fDiff+fMean) ;
112 rootSum += tempFunction[i] * std::cos(co << 113 }
113 } << 114 for (j=0;j<fNumber;j++)
114 fChebyshevCof[j] = weight * rootSum; // c << 115 {
115 } << 116 cofj = cof*j ;
116 // Chebyshev coefficients for (mx)-derivativ << 117 rootSum = 0.0 ;
117 << 118
118 for(i = 1; i <= mx; ++i) << 119 for (i=0;i<fNumber;i++)
119 { << 120 {
120 DerivativeChebyshevCof(tempFunction); << 121 rootSum += tempFunction[i]*cos(cofj*(i+0.5)) ;
121 fNumber--; << 122 }
122 for(j = 0; j < fNumber; ++j) << 123 fChebyshevCof[j] = weight*rootSum ; // corresponds to pFunction
123 { << 124 }
124 fChebyshevCof[j] = tempFunction[j]; // << 125 // Chebyshev coefficients for (m)-derivative of pFunction
125 } << 126
126 } << 127 for(i=1;i<=m;i++)
127 delete[] tempFunction; // delete of dynamic << 128 {
>> 129 DerivativeChebyshevCof(tempFunction) ;
>> 130 fNumber-- ;
>> 131 for(j=0;j<fNumber;j++)
>> 132 {
>> 133 fChebyshevCof[j] = tempFunction[j] ; // corresponds to (i)-derivative
>> 134 }
>> 135 }
>> 136 delete[] tempFunction ; // delete of dynamically allocated tempFunction
128 } 137 }
129 138
130 // ------------------------------------------- 139 // ------------------------------------------------------
131 // 140 //
132 // Constructor for creation of Chebyshev coeff 141 // Constructor for creation of Chebyshev coefficients for integral
133 // from pFunction. 142 // from pFunction.
134 << 143
135 G4ChebyshevApproximation::G4ChebyshevApproxima << 144 G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
136 << 145 G4double a,
137 << 146 G4double b,
138 : fFunction(pFunction) << 147 G4int n )
139 , fNumber(n) <<
140 , fChebyshevCof(new G4double[fNumber]) <<
141 , fMean(0.5 * (b + a)) <<
142 , fDiff(0.5 * (b - a)) <<
143 { 148 {
144 G4int i = 0, j = 0; << 149 G4int i,j;
145 G4double rootSum = 0.0, cofj = 0.0; << 150 G4double rootSum, cof, cofj, weight ;
146 auto tempFunction = new G4double[fNumber]; << 151
147 G4double weight = 2.0 / fNumber; << 152 fFunction = pFunction ;
148 G4double cof = 0.5 * weight * pi; << 153 fNumber = n ;
149 << 154 fDiff = 0.5*(b-a) ;
150 for(i = 0; i < fNumber; ++i) << 155 fMean = 0.5*(b+a) ;
151 { << 156 fChebyshevCof = new G4double[fNumber] ;
152 rootSum = std::cos(cof * (i + 0.5) << 157
153 tempFunction[i] = fFunction(rootSum * fDif << 158 G4double* tempFunction = new G4double[fNumber] ;
154 } << 159
155 for(j = 0; j < fNumber; ++j) << 160 weight = 2.0/fNumber ;
156 { << 161 cof = 0.5*weight*pi ; // pi/n
157 cofj = cof * j; << 162
158 rootSum = 0.0; << 163 for (i=0;i<fNumber;i++)
159 << 164 {
160 for(i = 0; i < fNumber; ++i) << 165 rootSum = cos(cof*(i+0.5)) ;
161 { << 166 tempFunction[i]= fFunction(rootSum*fDiff+fMean) ;
162 rootSum += tempFunction[i] * std::cos(co << 167 }
163 } << 168 for (j=0;j<fNumber;j++)
164 fChebyshevCof[j] = weight * rootSum; // c << 169 {
165 } << 170 cofj = cof*j ;
166 // Chebyshev coefficients for integral of pF << 171 rootSum = 0.0 ;
167 << 172
168 IntegralChebyshevCof(tempFunction); << 173 for (i=0;i<fNumber;i++)
169 for(j = 0; j < fNumber; ++j) << 174 {
170 { << 175 rootSum += tempFunction[i]*cos(cofj*(i+0.5)) ;
171 fChebyshevCof[j] = tempFunction[j]; // co << 176 }
172 } << 177 fChebyshevCof[j] = weight*rootSum ; // corresponds to pFunction
173 delete[] tempFunction; // delete of dynamic << 178 }
>> 179 // Chebyshev coefficients for integral of pFunction
>> 180
>> 181 IntegralChebyshevCof(tempFunction) ;
>> 182 for(j=0;j<fNumber;j++)
>> 183 {
>> 184 fChebyshevCof[j] = tempFunction[j] ; // corresponds to integral
>> 185 }
>> 186 delete[] tempFunction ; // delete of dynamically allocated tempFunction
174 } 187 }
175 188
>> 189
>> 190
176 // ------------------------------------------- 191 // ---------------------------------------------------------------
177 // 192 //
178 // Destructor deletes the array of Chebyshev c 193 // Destructor deletes the array of Chebyshev coefficients
179 194
180 G4ChebyshevApproximation::~G4ChebyshevApproxim 195 G4ChebyshevApproximation::~G4ChebyshevApproximation()
181 { 196 {
182 delete[] fChebyshevCof; << 197 delete[] fChebyshevCof ;
183 } 198 }
184 199
185 // ------------------------------------------- 200 // ---------------------------------------------------------------
186 // 201 //
187 // Access function for Chebyshev coefficients 202 // Access function for Chebyshev coefficients
188 // 203 //
189 204
190 G4double G4ChebyshevApproximation::GetChebyshe << 205
>> 206 G4double
>> 207 G4ChebyshevApproximation::GetChebyshevCof(G4int number) const
191 { 208 {
192 if(number < 0 && number >= fNumber) << 209 if(number < 0 && number >= fNumber)
193 { << 210 {
194 G4Exception("G4ChebyshevApproximation::Get << 211 G4Exception
195 FatalException, "Argument out << 212 ("Argument out of range in G4ChebyshevApproximation::GetChebyshevCof") ;
196 } << 213 }
197 return fChebyshevCof[number]; << 214 return fChebyshevCof[number] ;
198 } 215 }
199 216
200 // ------------------------------------------- 217 // --------------------------------------------------------------
201 // 218 //
202 // Evaluate the value of fFunction at the poin 219 // Evaluate the value of fFunction at the point x via the Chebyshev coefficients
203 // fChebyshevCof[0,...,fNumber-1] 220 // fChebyshevCof[0,...,fNumber-1]
204 221
205 G4double G4ChebyshevApproximation::ChebyshevEv << 222 G4double
>> 223 G4ChebyshevApproximation::ChebyshevEvaluation(G4double x) const
206 { 224 {
207 G4double evaluate = 0.0, evaluate2 = 0.0, te << 225 G4int i;
208 xReduced2 = 0.0; << 226 G4double evaluate = 0.0, evaluate2 = 0.0, temp, xReduced, xReduced2 ;
209 227
210 if((x - fMean + fDiff) * (x - fMean - fDiff) << 228 if ((x-fMean+fDiff)*(x-fMean-fDiff) > 0.0)
211 { << 229 {
212 G4Exception("G4ChebyshevApproximation::Che << 230 G4Exception("Invalid argument in G4ChebyshevApproximation::ChebyshevEvaluation");
213 "InvalidCall", FatalException, << 231 }
214 } << 232 xReduced = (x-fMean)/fDiff ;
215 xReduced = (x - fMean) / fDiff; << 233 xReduced2 = 2.0*xReduced ;
216 xReduced2 = 2.0 * xReduced; << 234 for (i=fNumber-1;i>=1;i--)
217 for(G4int i = fNumber - 1; i >= 1; --i) << 235 {
218 { << 236 temp = evaluate ;
219 temp = evaluate; << 237 evaluate = xReduced2*evaluate - evaluate2 + fChebyshevCof[i] ;
220 evaluate = xReduced2 * evaluate - evaluat << 238 evaluate2 = temp ;
221 evaluate2 = temp; << 239 }
222 } << 240 return xReduced*evaluate - evaluate2 + 0.5*fChebyshevCof[0] ;
223 return xReduced * evaluate - evaluate2 + 0.5 <<
224 } 241 }
225 242
226 // ------------------------------------------- 243 // ------------------------------------------------------------------
227 // 244 //
228 // Returns the array derCof[0,...,fNumber-2], << 245 // Returns the array derCof[0,...,fNumber-2], the Chebyshev coefficients of the
229 // derivative of the function whose coefficien 246 // derivative of the function whose coefficients are fChebyshevCof
230 247
231 void G4ChebyshevApproximation::DerivativeCheby << 248 void
>> 249 G4ChebyshevApproximation::DerivativeChebyshevCof(G4double derCof[]) const
232 { 250 {
233 G4double cof = 1.0 / fDiff; << 251 G4int i ;
234 derCof[fNumber - 1] = 0.0; << 252 G4double cof = 1.0/fDiff ;
235 derCof[fNumber - 2] = 2 * (fNumber - 1) * fC << 253 derCof[fNumber-1] = 0.0 ;
236 for(G4int i = fNumber - 3; i >= 0; --i) << 254 derCof[fNumber-2] = 2*(fNumber-1)*fChebyshevCof[fNumber-1] ;
237 { << 255 for(i=fNumber-3;i>=0;i--)
238 derCof[i] = derCof[i + 2] + 2 * (i + 1) * << 256 {
239 } << 257 derCof[i] = derCof[i+2] + 2*(i+1)*fChebyshevCof[i+1] ;
240 for(G4int j = 0; j < fNumber; ++j) << 258 }
241 { << 259 for(i=0;i<fNumber;i++)
242 derCof[j] *= cof; << 260 {
243 } << 261 derCof[i] *= cof ;
>> 262 }
244 } 263 }
245 264
246 // ------------------------------------------- 265 // ------------------------------------------------------------------------
247 // 266 //
248 // This function produces the array integralCo 267 // This function produces the array integralCof[0,...,fNumber-1] , the Chebyshev
249 // coefficients of the integral of the functio << 268 // coefficients of the integral of the function whose coefficients are
250 // fChebyshevCof[]. The constant of integratio << 269 // fChebyshevCof[]. The constant of integration is set so that the integral
251 // vanishes at the point (fMean - fDiff), i.e. << 270 // vanishes at the point (fMean - fDiff), i.e. at the begining of the interval of
252 // of validity (we start the integration from << 271 // validity (we start the integration from this point).
253 // << 272 //
254 << 273
255 void G4ChebyshevApproximation::IntegralChebysh << 274 void
256 G4double integralCof[]) const << 275 G4ChebyshevApproximation::IntegralChebyshevCof(G4double integralCof[]) const
257 { 276 {
258 G4double cof = 0.5 * fDiff, sum = 0.0, facto << 277 G4int i ;
259 for(G4int i = 1; i < fNumber - 1; ++i) << 278 G4double cof = 0.5*fDiff, sum = 0.0, factor = 1.0 ;
260 { << 279 for(i=1;i<fNumber-1;i++)
261 integralCof[i] = cof * (fChebyshevCof[i - << 280 {
262 sum += factor * integralCof[i]; << 281 integralCof[i] = cof*(fChebyshevCof[i-1] - fChebyshevCof[i+1])/i ;
263 factor = -factor; << 282 sum += factor*integralCof[i] ;
264 } << 283 factor = -factor ;
265 integralCof[fNumber - 1] = cof * fChebyshevC << 284 }
266 sum += factor * integralCof[fNumber - 1]; << 285 integralCof[fNumber-1] = cof*fChebyshevCof[fNumber-2]/(fNumber-1) ;
267 integralCof[0] = 2.0 * sum; // set the cons << 286 sum += factor*integralCof[fNumber-1] ;
268 } << 287 integralCof[0] = 2.0*sum ; // set the constant of integration
>> 288 }
269 289