| Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of
>> 2 // the GEANT4 collaboration.
1 // 3 //
2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work
3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement,
4 // * << 6 // and all its terms.
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 // 7 //
26 // G4ChebyshevApproximation class implementati << 8 // $Id: G4ChebyshevApproximation.cc,v 1.2 1999/11/16 17:31:09 gcosmo Exp $
>> 9 // GEANT4 tag $Name: geant4-03-01 $
27 // 10 //
28 // Author: V.Grichine, 24.04.97 <<
29 // ------------------------------------------- <<
30 11
31 #include "G4ChebyshevApproximation.hh" 12 #include "G4ChebyshevApproximation.hh"
32 #include "G4PhysicalConstants.hh" <<
33 13
34 // Constructor for initialisation of the class << 14
35 // It creates the array fChebyshevCof[0,...,fN << 15 // Constructor for initialisation of the class data members. It creates the array
36 // which consists of Chebyshev coefficients de << 16 // fChebyshevCof[0,...,fNumber-1], fNumber = n ; which consists of Chebyshev
37 // pointed by pFunction. The values a and b fi << 17 // coefficients describing the function pointed by pFunction. The values a and b
38 // of the Chebyshev approximation. << 18 // fixe the interval of validity of Chebyshev approximation.
39 << 19
40 G4ChebyshevApproximation::G4ChebyshevApproxima << 20
41 << 21 G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
42 : fFunction(pFunction) << 22 G4int n,
43 , fNumber(n) << 23 G4double a,
44 , fChebyshevCof(new G4double[fNumber]) << 24 G4double b )
45 , fMean(0.5 * (b + a)) <<
46 , fDiff(0.5 * (b - a)) <<
47 { 25 {
48 G4int i = 0, j = 0; << 26 G4int i, j ;
49 G4double rootSum = 0.0, cofj = 0.0; << 27 G4double rootSum, cof, cofj, weight ;
50 auto tempFunction = new G4double[fNumber]; << 28
51 G4double weight = 2.0 / fNumber; << 29 fFunction = pFunction ;
52 G4double cof = 0.5 * weight * pi; << 30 fNumber = n ;
53 << 31 fDiff = 0.5*(b-a) ;
54 for(i = 0; i < fNumber; ++i) << 32 fMean = 0.5*(b+a) ;
55 { << 33 fChebyshevCof = new G4double[fNumber] ;
56 rootSum = std::cos(cof * (i + 0.5) << 34
57 tempFunction[i] = fFunction(rootSum * fDif << 35 G4double* tempFunction = new G4double[fNumber] ;
58 } << 36
59 for(j = 0; j < fNumber; ++j) << 37 weight = 2.0/fNumber ;
60 { << 38 cof = 0.5*weight*pi ; // pi/n
61 cofj = cof * j; << 39
62 rootSum = 0.0; << 40 for (i=0;i<fNumber;i++)
63 << 41 {
64 for(i = 0; i < fNumber; ++i) << 42 rootSum = cos(cof*(i+0.5)) ;
65 { << 43 tempFunction[i]= fFunction(rootSum*fDiff+fMean) ;
66 rootSum += tempFunction[i] * std::cos(co << 44 }
67 } << 45 for (j=0;j<fNumber;j++)
68 fChebyshevCof[j] = weight * rootSum; << 46 {
69 } << 47 cofj = cof*j ;
70 delete[] tempFunction; << 48 rootSum = 0.0 ;
>> 49
>> 50 for (i=0;i<fNumber;i++)
>> 51 {
>> 52 rootSum += tempFunction[i]*cos(cofj*(i+0.5)) ;
>> 53 }
>> 54 fChebyshevCof[j] = weight*rootSum ;
>> 55 }
>> 56 delete[] tempFunction ;
71 } 57 }
72 58
73 // ------------------------------------------- 59 // --------------------------------------------------------------------
74 // 60 //
75 // Constructor for creation of Chebyshev coeff << 61 // Constructor for creation of Chebyshev coefficients for m-derivative
76 // from pFunction. The value of mx ! MUST BE ! << 62 // from pFunction. The value of m ! MUST BE ! < n , because the result
77 // array of fChebyshevCof will be of (nx-mx) s << 63 // array of fChebyshevCof will be of (n-m) size. The values a and b
78 // fix the interval of validity of the Chebysh << 64 // fixe the interval of validity of Chebyshev approximation.
79 << 65
80 G4ChebyshevApproximation::G4ChebyshevApproxima << 66
81 << 67 G4ChebyshevApproximation::
82 << 68 G4ChebyshevApproximation( function pFunction,
83 : fFunction(pFunction) << 69 G4int n,
84 , fNumber(nx) << 70 G4int m,
85 , fChebyshevCof(new G4double[fNumber]) << 71 G4double a,
86 , fMean(0.5 * (b + a)) << 72 G4double b )
87 , fDiff(0.5 * (b - a)) <<
88 { 73 {
89 if(nx <= mx) << 74 if(n <= m)
90 { << 75 {
91 G4Exception("G4ChebyshevApproximation::G4C << 76 G4Exception
92 "InvalidCall", FatalException, << 77 ("Invalid arguments in G4ChebyshevApproximation::G4ChebyshevApproximation") ;
93 } << 78 }
94 G4int i = 0, j = 0; << 79 G4int i, j ;
95 G4double rootSum = 0.0, cofj = 0.0; << 80 G4double rootSum, cof, cofj, weight ;
96 auto tempFunction = new G4double[fNumber]; << 81
97 G4double weight = 2.0 / fNumber; << 82 fFunction = pFunction ;
98 G4double cof = 0.5 * weight * pi; << 83 fNumber = n ;
99 << 84 fDiff = 0.5*(b-a) ;
100 for(i = 0; i < fNumber; ++i) << 85 fMean = 0.5*(b+a) ;
101 { << 86 fChebyshevCof = new G4double[fNumber] ;
102 rootSum = std::cos(cof * (i + 0.5) << 87
103 tempFunction[i] = fFunction(rootSum * fDif << 88 G4double* tempFunction = new G4double[fNumber] ;
104 } << 89
105 for(j = 0; j < fNumber; ++j) << 90 weight = 2.0/fNumber ;
106 { << 91 cof = 0.5*weight*pi ; // pi/n
107 cofj = cof * j; << 92
108 rootSum = 0.0; << 93 for (i=0;i<fNumber;i++)
109 << 94 {
110 for(i = 0; i < fNumber; ++i) << 95 rootSum = cos(cof*(i+0.5)) ;
111 { << 96 tempFunction[i] = fFunction(rootSum*fDiff+fMean) ;
112 rootSum += tempFunction[i] * std::cos(co << 97 }
113 } << 98 for (j=0;j<fNumber;j++)
114 fChebyshevCof[j] = weight * rootSum; // c << 99 {
115 } << 100 cofj = cof*j ;
116 // Chebyshev coefficients for (mx)-derivativ << 101 rootSum = 0.0 ;
117 << 102
118 for(i = 1; i <= mx; ++i) << 103 for (i=0;i<fNumber;i++)
119 { << 104 {
120 DerivativeChebyshevCof(tempFunction); << 105 rootSum += tempFunction[i]*cos(cofj*(i+0.5)) ;
121 fNumber--; << 106 }
122 for(j = 0; j < fNumber; ++j) << 107 fChebyshevCof[j] = weight*rootSum ; // corresponds to pFunction
123 { << 108 }
124 fChebyshevCof[j] = tempFunction[j]; // << 109 // Chebyshev coefficients for (m)-derivative of pFunction
125 } << 110
126 } << 111 for(i=1;i<=m;i++)
127 delete[] tempFunction; // delete of dynamic << 112 {
>> 113 DerivativeChebyshevCof(tempFunction) ;
>> 114 fNumber-- ;
>> 115 for(j=0;j<fNumber;j++)
>> 116 {
>> 117 fChebyshevCof[j] = tempFunction[j] ; // corresponds to (i)-derivative
>> 118 }
>> 119 }
>> 120 delete[] tempFunction ; // delete of dynamically allocated tempFunction
128 } 121 }
129 122
130 // ------------------------------------------- 123 // ------------------------------------------------------
131 // 124 //
132 // Constructor for creation of Chebyshev coeff 125 // Constructor for creation of Chebyshev coefficients for integral
133 // from pFunction. 126 // from pFunction.
134 << 127
135 G4ChebyshevApproximation::G4ChebyshevApproxima << 128 G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
136 << 129 G4double a,
137 << 130 G4double b,
138 : fFunction(pFunction) << 131 G4int n )
139 , fNumber(n) <<
140 , fChebyshevCof(new G4double[fNumber]) <<
141 , fMean(0.5 * (b + a)) <<
142 , fDiff(0.5 * (b - a)) <<
143 { 132 {
144 G4int i = 0, j = 0; << 133 G4int i,j;
145 G4double rootSum = 0.0, cofj = 0.0; << 134 G4double rootSum, cof, cofj, weight ;
146 auto tempFunction = new G4double[fNumber]; << 135
147 G4double weight = 2.0 / fNumber; << 136 fFunction = pFunction ;
148 G4double cof = 0.5 * weight * pi; << 137 fNumber = n ;
149 << 138 fDiff = 0.5*(b-a) ;
150 for(i = 0; i < fNumber; ++i) << 139 fMean = 0.5*(b+a) ;
151 { << 140 fChebyshevCof = new G4double[fNumber] ;
152 rootSum = std::cos(cof * (i + 0.5) << 141
153 tempFunction[i] = fFunction(rootSum * fDif << 142 G4double* tempFunction = new G4double[fNumber] ;
154 } << 143
155 for(j = 0; j < fNumber; ++j) << 144 weight = 2.0/fNumber ;
156 { << 145 cof = 0.5*weight*pi ; // pi/n
157 cofj = cof * j; << 146
158 rootSum = 0.0; << 147 for (i=0;i<fNumber;i++)
159 << 148 {
160 for(i = 0; i < fNumber; ++i) << 149 rootSum = cos(cof*(i+0.5)) ;
161 { << 150 tempFunction[i]= fFunction(rootSum*fDiff+fMean) ;
162 rootSum += tempFunction[i] * std::cos(co << 151 }
163 } << 152 for (j=0;j<fNumber;j++)
164 fChebyshevCof[j] = weight * rootSum; // c << 153 {
165 } << 154 cofj = cof*j ;
166 // Chebyshev coefficients for integral of pF << 155 rootSum = 0.0 ;
167 << 156
168 IntegralChebyshevCof(tempFunction); << 157 for (i=0;i<fNumber;i++)
169 for(j = 0; j < fNumber; ++j) << 158 {
170 { << 159 rootSum += tempFunction[i]*cos(cofj*(i+0.5)) ;
171 fChebyshevCof[j] = tempFunction[j]; // co << 160 }
172 } << 161 fChebyshevCof[j] = weight*rootSum ; // corresponds to pFunction
173 delete[] tempFunction; // delete of dynamic << 162 }
>> 163 // Chebyshev coefficients for integral of pFunction
>> 164
>> 165 IntegralChebyshevCof(tempFunction) ;
>> 166 for(j=0;j<fNumber;j++)
>> 167 {
>> 168 fChebyshevCof[j] = tempFunction[j] ; // corresponds to integral
>> 169 }
>> 170 delete[] tempFunction ; // delete of dynamically allocated tempFunction
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
195 FatalException, "Argument out << 196 ("Argument out of range in G4ChebyshevApproximation::GetChebyshevCof") ;
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 G4int i;
208 xReduced2 = 0.0; << 210 G4double evaluate = 0.0, evaluate2 = 0.0, temp, xReduced, xReduced2 ;
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("Invalid argument in G4ChebyshevApproximation::ChebyshevEvaluation");
213 "InvalidCall", FatalException, << 215 }
214 } << 216 xReduced = (x-fMean)/fDiff ;
215 xReduced = (x - fMean) / fDiff; << 217 xReduced2 = 2.0*xReduced ;
216 xReduced2 = 2.0 * xReduced; << 218 for (i=fNumber-1;i>=1;i--)
217 for(G4int i = fNumber - 1; i >= 1; --i) << 219 {
218 { << 220 temp = evaluate ;
219 temp = evaluate; << 221 evaluate = xReduced2*evaluate - evaluate2 + fChebyshevCof[i] ;
220 evaluate = xReduced2 * evaluate - evaluat << 222 evaluate2 = temp ;
221 evaluate2 = temp; << 223 }
222 } << 224 return xReduced*evaluate - evaluate2 + 0.5*fChebyshevCof[0] ;
223 return xReduced * evaluate - evaluate2 + 0.5 <<
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 G4int i ;
234 derCof[fNumber - 1] = 0.0; << 236 G4double cof = 1.0/fDiff ;
235 derCof[fNumber - 2] = 2 * (fNumber - 1) * fC << 237 derCof[fNumber-1] = 0.0 ;
236 for(G4int i = fNumber - 3; i >= 0; --i) << 238 derCof[fNumber-2] = 2*(fNumber-1)*fChebyshevCof[fNumber-1] ;
237 { << 239 for(i=fNumber-3;i>=0;i--)
238 derCof[i] = derCof[i + 2] + 2 * (i + 1) * << 240 {
239 } << 241 derCof[i] = derCof[i+2] + 2*(i+1)*fChebyshevCof[i+1] ;
240 for(G4int j = 0; j < fNumber; ++j) << 242 }
241 { << 243 for(i=0;i<fNumber;i++)
242 derCof[j] *= cof; << 244 {
243 } << 245 derCof[i] *= cof ;
>> 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 G4int i ;
259 for(G4int i = 1; i < fNumber - 1; ++i) << 262 G4double cof = 0.5*fDiff, sum = 0.0, factor = 1.0 ;
260 { << 263 for(i=1;i<fNumber-1;i++)
261 integralCof[i] = cof * (fChebyshevCof[i - << 264 {
262 sum += factor * integralCof[i]; << 265 integralCof[i] = cof*(fChebyshevCof[i-1] - fChebyshevCof[i+1])/i ;
263 factor = -factor; << 266 sum += factor*integralCof[i] ;
264 } << 267 factor = -factor ;
265 integralCof[fNumber - 1] = cof * fChebyshevC << 268 }
266 sum += factor * integralCof[fNumber - 1]; << 269 integralCof[fNumber-1] = cof*fChebyshevCof[fNumber-2]/(fNumber-1) ;
267 integralCof[0] = 2.0 * sum; // set the cons << 270 sum += factor*integralCof[fNumber-1] ;
268 } << 271 integralCof[0] = 2.0*sum ; // set the constant of integration
>> 272 }
269 273