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: geant4-09-02-patch-03 $ >> 29 // 30 30 31 #include "G4ChebyshevApproximation.hh" 31 #include "G4ChebyshevApproximation.hh" 32 #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