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 // G4DataInterpolation class implementation 26 // G4DataInterpolation class implementation 27 // 27 // 28 // Author: V.Grichine, 03.04.1997 28 // Author: V.Grichine, 03.04.1997 29 // ------------------------------------------- 29 // -------------------------------------------------------------------- 30 30 31 #include "G4DataInterpolation.hh" 31 #include "G4DataInterpolation.hh" 32 32 33 ////////////////////////////////////////////// 33 ////////////////////////////////////////////////////////////////////////////// 34 // 34 // 35 // Constructor for initializing of fArgument, 35 // Constructor for initializing of fArgument, fFunction and fNumber 36 // data members 36 // data members 37 37 38 G4DataInterpolation::G4DataInterpolation(G4dou 38 G4DataInterpolation::G4DataInterpolation(G4double pX[], G4double pY[], 39 G4int 39 G4int number) 40 : fArgument(new G4double[number]) 40 : fArgument(new G4double[number]) 41 , fFunction(new G4double[number]) 41 , fFunction(new G4double[number]) 42 , fNumber(number) 42 , fNumber(number) 43 { 43 { 44 for(G4int i = 0; i < fNumber; ++i) 44 for(G4int i = 0; i < fNumber; ++i) 45 { 45 { 46 fArgument[i] = pX[i]; 46 fArgument[i] = pX[i]; 47 fFunction[i] = pY[i]; 47 fFunction[i] = pY[i]; 48 } 48 } 49 } 49 } 50 50 51 ////////////////////////////////////////////// 51 //////////////////////////////////////////////////////////////////////////// 52 // 52 // 53 // Constructor for cubic spline interpolation. 53 // Constructor for cubic spline interpolation. It creates the array 54 // fSecondDerivative[0,...fNumber-1] which is 54 // fSecondDerivative[0,...fNumber-1] which is used in this interpolation by 55 // the function 55 // the function 56 56 57 G4DataInterpolation::G4DataInterpolation(G4dou 57 G4DataInterpolation::G4DataInterpolation(G4double pX[], G4double pY[], 58 G4int 58 G4int number, G4double pFirstDerStart, 59 G4dou 59 G4double pFirstDerFinish) 60 : fArgument(new G4double[number]) 60 : fArgument(new G4double[number]) 61 , fFunction(new G4double[number]) 61 , fFunction(new G4double[number]) 62 , fSecondDerivative(new G4double[number]) 62 , fSecondDerivative(new G4double[number]) 63 , fNumber(number) 63 , fNumber(number) 64 { 64 { 65 G4int i = 0; 65 G4int i = 0; 66 G4double p = 0.0, qn = 0.0, sig = 0.0, un = 66 G4double p = 0.0, qn = 0.0, sig = 0.0, un = 0.0; 67 const G4double maxDerivative = 0.99e30; 67 const G4double maxDerivative = 0.99e30; 68 auto u = new G4double[fNumb << 68 G4double* u = new G4double[fNumber - 1]; 69 69 70 for(i = 0; i < fNumber; ++i) 70 for(i = 0; i < fNumber; ++i) 71 { 71 { 72 fArgument[i] = pX[i]; 72 fArgument[i] = pX[i]; 73 fFunction[i] = pY[i]; 73 fFunction[i] = pY[i]; 74 } 74 } 75 if(pFirstDerStart > maxDerivative) 75 if(pFirstDerStart > maxDerivative) 76 { 76 { 77 fSecondDerivative[0] = 0.0; 77 fSecondDerivative[0] = 0.0; 78 u[0] = 0.0; 78 u[0] = 0.0; 79 } 79 } 80 else 80 else 81 { 81 { 82 fSecondDerivative[0] = -0.5; 82 fSecondDerivative[0] = -0.5; 83 u[0] = (3.0 / (fArgument[1 83 u[0] = (3.0 / (fArgument[1] - fArgument[0])) * 84 ((fFunction[1] - fFunction[0]) / (f 84 ((fFunction[1] - fFunction[0]) / (fArgument[1] - fArgument[0]) - 85 pFirstDerStart); 85 pFirstDerStart); 86 } 86 } 87 87 88 // Decomposition loop for tridiagonal algori 88 // Decomposition loop for tridiagonal algorithm. fSecondDerivative[i] 89 // and u[i] are used for temporary storage o 89 // and u[i] are used for temporary storage of the decomposed factors. 90 90 91 for(i = 1; i < fNumber - 1; ++i) 91 for(i = 1; i < fNumber - 1; ++i) 92 { 92 { 93 sig = 93 sig = 94 (fArgument[i] - fArgument[i - 1]) / (fAr 94 (fArgument[i] - fArgument[i - 1]) / (fArgument[i + 1] - fArgument[i - 1]); 95 p = sig * fSecondDeriva 95 p = sig * fSecondDerivative[i - 1] + 2.0; 96 fSecondDerivative[i] = (sig - 1.0) / p; 96 fSecondDerivative[i] = (sig - 1.0) / p; 97 u[i] = 97 u[i] = 98 (fFunction[i + 1] - fFunction[i]) / (fAr 98 (fFunction[i + 1] - fFunction[i]) / (fArgument[i + 1] - fArgument[i]) - 99 (fFunction[i] - fFunction[i - 1]) / (fAr 99 (fFunction[i] - fFunction[i - 1]) / (fArgument[i] - fArgument[i - 1]); 100 u[i] = 100 u[i] = 101 (6.0 * u[i] / (fArgument[i + 1] - fArgum 101 (6.0 * u[i] / (fArgument[i + 1] - fArgument[i - 1]) - sig * u[i - 1]) / p; 102 } 102 } 103 if(pFirstDerFinish > maxDerivative) 103 if(pFirstDerFinish > maxDerivative) 104 { 104 { 105 qn = 0.0; 105 qn = 0.0; 106 un = 0.0; 106 un = 0.0; 107 } 107 } 108 else 108 else 109 { 109 { 110 qn = 0.5; 110 qn = 0.5; 111 un = 111 un = 112 (3.0 / (fArgument[fNumber - 1] - fArgume 112 (3.0 / (fArgument[fNumber - 1] - fArgument[fNumber - 2])) * 113 (pFirstDerFinish - (fFunction[fNumber - 113 (pFirstDerFinish - (fFunction[fNumber - 1] - fFunction[fNumber - 2]) / 114 (fArgument[fNumber 114 (fArgument[fNumber - 1] - fArgument[fNumber - 2])); 115 } 115 } 116 fSecondDerivative[fNumber - 1] = 116 fSecondDerivative[fNumber - 1] = 117 (un - qn * u[fNumber - 2]) / (qn * fSecond 117 (un - qn * u[fNumber - 2]) / (qn * fSecondDerivative[fNumber - 2] + 1.0); 118 118 119 // The backsubstitution loop for the triagon 119 // The backsubstitution loop for the triagonal algorithm of solving 120 // a linear system of equations. 120 // a linear system of equations. 121 121 122 for(G4int k = fNumber - 2; k >= 0; --k) 122 for(G4int k = fNumber - 2; k >= 0; --k) 123 { 123 { 124 fSecondDerivative[k] = 124 fSecondDerivative[k] = 125 fSecondDerivative[k] * fSecondDerivative 125 fSecondDerivative[k] * fSecondDerivative[k + 1] + u[k]; 126 } 126 } 127 delete[] u; 127 delete[] u; 128 } 128 } 129 129 130 ////////////////////////////////////////////// 130 ///////////////////////////////////////////////////////////////////////////// 131 // 131 // 132 // Destructor deletes dynamically created arra 132 // Destructor deletes dynamically created arrays for data members: fArgument, 133 // fFunction and fSecondDerivative, all have d 133 // fFunction and fSecondDerivative, all have dimension of fNumber 134 134 135 G4DataInterpolation::~G4DataInterpolation() 135 G4DataInterpolation::~G4DataInterpolation() 136 { 136 { 137 delete[] fArgument; 137 delete[] fArgument; 138 delete[] fFunction; 138 delete[] fFunction; 139 delete[] fSecondDerivative; 139 delete[] fSecondDerivative; 140 } 140 } 141 141 142 ////////////////////////////////////////////// 142 ///////////////////////////////////////////////////////////////////////////// 143 // 143 // 144 // This function returns the value P(pX), wher 144 // This function returns the value P(pX), where P(x) is polynom of fNumber-1 145 // degree such that P(fArgument[i]) = fFunctio 145 // degree such that P(fArgument[i]) = fFunction[i], for i = 0, ..., fNumber-1. 146 // This is Lagrange's form of interpolation an 146 // This is Lagrange's form of interpolation and it is based on Neville's 147 // algorithm 147 // algorithm 148 148 149 G4double G4DataInterpolation::PolynomInterpola 149 G4double G4DataInterpolation::PolynomInterpolation(G4double pX, 150 150 G4double& deltaY) const 151 { 151 { 152 G4int i = 0, j = 1, k = 0; 152 G4int i = 0, j = 1, k = 0; 153 G4double mult = 0.0, difi = 0.0, deltaLow = 153 G4double mult = 0.0, difi = 0.0, deltaLow = 0.0, deltaUp = 0.0, cd = 0.0, 154 y = 0.0; 154 y = 0.0; 155 auto c = new G4double[fNumber]; << 155 G4double* c = new G4double[fNumber]; 156 auto d = new G4double[fNumber]; << 156 G4double* d = new G4double[fNumber]; 157 G4double diff = std::fabs(pX - fArgument[0]) 157 G4double diff = std::fabs(pX - fArgument[0]); 158 for(i = 0; i < fNumber; ++i) 158 for(i = 0; i < fNumber; ++i) 159 { 159 { 160 difi = std::fabs(pX - fArgument[i]); 160 difi = std::fabs(pX - fArgument[i]); 161 if(difi < diff) 161 if(difi < diff) 162 { 162 { 163 k = i; 163 k = i; 164 diff = difi; 164 diff = difi; 165 } 165 } 166 c[i] = fFunction[i]; 166 c[i] = fFunction[i]; 167 d[i] = fFunction[i]; 167 d[i] = fFunction[i]; 168 } 168 } 169 y = fFunction[k--]; 169 y = fFunction[k--]; 170 for(j = 1; j < fNumber; ++j) 170 for(j = 1; j < fNumber; ++j) 171 { 171 { 172 for(i = 0; i < fNumber - j; ++i) 172 for(i = 0; i < fNumber - j; ++i) 173 { 173 { 174 deltaLow = fArgument[i] - pX; 174 deltaLow = fArgument[i] - pX; 175 deltaUp = fArgument[i + j] - pX; 175 deltaUp = fArgument[i + j] - pX; 176 cd = c[i + 1] - d[i]; 176 cd = c[i + 1] - d[i]; 177 mult = deltaLow - deltaUp; 177 mult = deltaLow - deltaUp; 178 if(!(mult != 0.0)) 178 if(!(mult != 0.0)) 179 { 179 { 180 G4Exception("G4DataInterpolation::Poly 180 G4Exception("G4DataInterpolation::PolynomInterpolation()", "Error", 181 FatalException, "Coinciden 181 FatalException, "Coincident nodes !"); 182 } 182 } 183 mult = cd / mult; 183 mult = cd / mult; 184 d[i] = deltaUp * mult; 184 d[i] = deltaUp * mult; 185 c[i] = deltaLow * mult; 185 c[i] = deltaLow * mult; 186 } 186 } 187 y += (deltaY = (2 * k < (fNumber - j - 1) 187 y += (deltaY = (2 * k < (fNumber - j - 1) ? c[k + 1] : d[k--])); 188 } 188 } 189 delete[] c; 189 delete[] c; 190 delete[] d; 190 delete[] d; 191 191 192 return y; 192 return y; 193 } 193 } 194 194 195 ////////////////////////////////////////////// 195 //////////////////////////////////////////////////////////////////////////// 196 // 196 // 197 // Given arrays fArgument[0,..,fNumber-1] and 197 // Given arrays fArgument[0,..,fNumber-1] and fFunction[0,..,fNumber-1], this 198 // function calculates an array of coefficient 198 // function calculates an array of coefficients. The coefficients don't provide 199 // usually (fNumber>10) better accuracy for po 199 // usually (fNumber>10) better accuracy for polynom interpolation, as compared 200 // with PolynomInterpolation function. They co 200 // with PolynomInterpolation function. They could be used instead for derivate 201 // calculations and some other applications. 201 // calculations and some other applications. 202 202 203 void G4DataInterpolation::PolIntCoefficient(G4 203 void G4DataInterpolation::PolIntCoefficient(G4double cof[]) const 204 { 204 { 205 G4int i = 0, j = 0; 205 G4int i = 0, j = 0; 206 G4double factor; 206 G4double factor; 207 G4double reducedY = 0.0, mult = 1.0; 207 G4double reducedY = 0.0, mult = 1.0; 208 auto tempArgument = new G4double[fNumber]; << 208 G4double* tempArgument = new G4double[fNumber]; 209 209 210 for(i = 0; i < fNumber; ++i) 210 for(i = 0; i < fNumber; ++i) 211 { 211 { 212 tempArgument[i] = cof[i] = 0.0; 212 tempArgument[i] = cof[i] = 0.0; 213 } 213 } 214 tempArgument[fNumber - 1] = -fArgument[0]; 214 tempArgument[fNumber - 1] = -fArgument[0]; 215 215 216 for(i = 1; i < fNumber; ++i) 216 for(i = 1; i < fNumber; ++i) 217 { 217 { 218 for(j = fNumber - 1 - i; j < fNumber - 1; 218 for(j = fNumber - 1 - i; j < fNumber - 1; ++j) 219 { 219 { 220 tempArgument[j] -= fArgument[i] * tempAr 220 tempArgument[j] -= fArgument[i] * tempArgument[j + 1]; 221 } 221 } 222 tempArgument[fNumber - 1] -= fArgument[i]; 222 tempArgument[fNumber - 1] -= fArgument[i]; 223 } 223 } 224 for(i = 0; i < fNumber; ++i) 224 for(i = 0; i < fNumber; ++i) 225 { 225 { 226 factor = fNumber; 226 factor = fNumber; 227 for(j = fNumber - 1; j >= 1; --j) 227 for(j = fNumber - 1; j >= 1; --j) 228 { 228 { 229 factor = j * tempArgument[j] + factor * 229 factor = j * tempArgument[j] + factor * fArgument[i]; 230 } 230 } 231 reducedY = fFunction[i] / factor; 231 reducedY = fFunction[i] / factor; 232 mult = 1.0; 232 mult = 1.0; 233 for(j = fNumber - 1; j >= 0; --j) 233 for(j = fNumber - 1; j >= 0; --j) 234 { 234 { 235 cof[j] += mult * reducedY; 235 cof[j] += mult * reducedY; 236 mult = tempArgument[j] + mult * fArgumen 236 mult = tempArgument[j] + mult * fArgument[i]; 237 } 237 } 238 } 238 } 239 delete[] tempArgument; 239 delete[] tempArgument; 240 } 240 } 241 241 242 ////////////////////////////////////////////// 242 ///////////////////////////////////////////////////////////////////////////// 243 // 243 // 244 // The function returns diagonal rational func 244 // The function returns diagonal rational function (Bulirsch and Stoer 245 // algorithm of Neville type) Pn(x)/Qm(x) wher 245 // algorithm of Neville type) Pn(x)/Qm(x) where P and Q are polynoms. 246 // Tests showed the method is not stable and h 246 // Tests showed the method is not stable and hasn't advantage if compared 247 // with polynomial interpolation ?! 247 // with polynomial interpolation ?! 248 248 249 G4double G4DataInterpolation::RationalPolInter 249 G4double G4DataInterpolation::RationalPolInterpolation(G4double pX, 250 250 G4double& deltaY) const 251 { 251 { 252 G4int i = 0, j = 1, k = 0; 252 G4int i = 0, j = 1, k = 0; 253 const G4double tolerance = 1.6e-24; 253 const G4double tolerance = 1.6e-24; 254 G4double mult = 0.0, difi = 0.0, cd = 0.0, y 254 G4double mult = 0.0, difi = 0.0, cd = 0.0, y = 0.0, cof = 0.0; 255 auto c = new G4double[fNumber]; << 255 G4double* c = new G4double[fNumber]; 256 auto d = new G4double[fNumber]; << 256 G4double* d = new G4double[fNumber]; 257 G4double diff = std::fabs(pX - fArgument[0]) 257 G4double diff = std::fabs(pX - fArgument[0]); 258 for(i = 0; i < fNumber; ++i) 258 for(i = 0; i < fNumber; ++i) 259 { 259 { 260 difi = std::fabs(pX - fArgument[i]); 260 difi = std::fabs(pX - fArgument[i]); 261 if(!(difi != 0.0)) 261 if(!(difi != 0.0)) 262 { 262 { 263 y = fFunction[i]; 263 y = fFunction[i]; 264 deltaY = 0.0; 264 deltaY = 0.0; 265 delete[] c; 265 delete[] c; 266 delete[] d; 266 delete[] d; 267 return y; 267 return y; 268 } 268 } 269 if(difi < diff) << 269 else if(difi < diff) 270 { 270 { 271 k = i; 271 k = i; 272 diff = difi; 272 diff = difi; 273 } 273 } 274 c[i] = fFunction[i]; 274 c[i] = fFunction[i]; 275 d[i] = fFunction[i] + tolerance; // to pr 275 d[i] = fFunction[i] + tolerance; // to prevent rare zero/zero cases 276 } 276 } 277 y = fFunction[k--]; 277 y = fFunction[k--]; 278 for(j = 1; j < fNumber; ++j) 278 for(j = 1; j < fNumber; ++j) 279 { 279 { 280 for(i = 0; i < fNumber - j; ++i) 280 for(i = 0; i < fNumber - j; ++i) 281 { 281 { 282 cd = c[i + 1] - d[i]; 282 cd = c[i + 1] - d[i]; 283 difi = fArgument[i + j] - pX; 283 difi = fArgument[i + j] - pX; 284 cof = (fArgument[i] - pX) * d[i] / difi 284 cof = (fArgument[i] - pX) * d[i] / difi; 285 mult = cof - c[i + 1]; 285 mult = cof - c[i + 1]; 286 if(!(mult != 0.0)) // function to be in 286 if(!(mult != 0.0)) // function to be interpolated has pole at pX 287 { 287 { 288 G4Exception("G4DataInterpolation::Rati 288 G4Exception("G4DataInterpolation::RationalPolInterpolation()", "Error", 289 FatalException, "Coinciden 289 FatalException, "Coincident nodes !"); 290 } 290 } 291 mult = cd / mult; 291 mult = cd / mult; 292 d[i] = c[i + 1] * mult; 292 d[i] = c[i + 1] * mult; 293 c[i] = cof * mult; 293 c[i] = cof * mult; 294 } 294 } 295 y += (deltaY = (2 * k < (fNumber - j - 1) 295 y += (deltaY = (2 * k < (fNumber - j - 1) ? c[k + 1] : d[k--])); 296 } 296 } 297 delete[] c; 297 delete[] c; 298 delete[] d; 298 delete[] d; 299 299 300 return y; 300 return y; 301 } 301 } 302 302 303 ////////////////////////////////////////////// 303 ///////////////////////////////////////////////////////////////////////////// 304 // 304 // 305 // Cubic spline interpolation in point pX for 305 // Cubic spline interpolation in point pX for function given by the table: 306 // fArgument, fFunction. The constructor, whic 306 // fArgument, fFunction. The constructor, which creates fSecondDerivative, 307 // must be called before. The function works o 307 // must be called before. The function works optimal, if sequential calls 308 // are in random values of pX. 308 // are in random values of pX. 309 309 310 G4double G4DataInterpolation::CubicSplineInter 310 G4double G4DataInterpolation::CubicSplineInterpolation(G4double pX) const 311 { 311 { 312 G4int kLow = 0, kHigh = fNumber - 1, k = 0; 312 G4int kLow = 0, kHigh = fNumber - 1, k = 0; 313 313 314 // Searching in the table by means of bisect 314 // Searching in the table by means of bisection method. 315 // fArgument must be monotonic, either incre 315 // fArgument must be monotonic, either increasing or decreasing 316 316 317 while((kHigh - kLow) > 1) 317 while((kHigh - kLow) > 1) 318 { 318 { 319 k = (kHigh + kLow) >> 1; // compute midpo 319 k = (kHigh + kLow) >> 1; // compute midpoint 'bisection' 320 if(fArgument[k] > pX) 320 if(fArgument[k] > pX) 321 { 321 { 322 kHigh = k; 322 kHigh = k; 323 } 323 } 324 else 324 else 325 { 325 { 326 kLow = k; 326 kLow = k; 327 } 327 } 328 } // kLow and kHigh now bracket the input v 328 } // kLow and kHigh now bracket the input value of pX 329 G4double deltaHL = fArgument[kHigh] - fArgum 329 G4double deltaHL = fArgument[kHigh] - fArgument[kLow]; 330 if(!(deltaHL != 0.0)) 330 if(!(deltaHL != 0.0)) 331 { 331 { 332 G4Exception("G4DataInterpolation::CubicSpl 332 G4Exception("G4DataInterpolation::CubicSplineInterpolation()", "Error", 333 FatalException, "Bad fArgument 333 FatalException, "Bad fArgument input !"); 334 } 334 } 335 G4double a = (fArgument[kHigh] - pX) / delta 335 G4double a = (fArgument[kHigh] - pX) / deltaHL; 336 G4double b = (pX - fArgument[kLow]) / deltaH 336 G4double b = (pX - fArgument[kLow]) / deltaHL; 337 337 338 // Final evaluation of cubic spline polynomi 338 // Final evaluation of cubic spline polynomial for return 339 339 340 return a * fFunction[kLow] + b * fFunction[k 340 return a * fFunction[kLow] + b * fFunction[kHigh] + 341 ((a * a * a - a) * fSecondDerivative[ 341 ((a * a * a - a) * fSecondDerivative[kLow] + 342 (b * b * b - b) * fSecondDerivative[ 342 (b * b * b - b) * fSecondDerivative[kHigh]) * 343 deltaHL * deltaHL / 6.0; 343 deltaHL * deltaHL / 6.0; 344 } 344 } 345 345 346 ////////////////////////////////////////////// 346 /////////////////////////////////////////////////////////////////////////// 347 // 347 // 348 // Return cubic spline interpolation in the po 348 // Return cubic spline interpolation in the point pX which is located between 349 // fArgument[index] and fArgument[index+1]. It 349 // fArgument[index] and fArgument[index+1]. It is usually called in sequence 350 // of known from external analysis values of i 350 // of known from external analysis values of index. 351 351 352 G4double G4DataInterpolation::FastCubicSpline( 352 G4double G4DataInterpolation::FastCubicSpline(G4double pX, G4int index) const 353 { 353 { 354 G4double delta = fArgument[index + 1] - fArg 354 G4double delta = fArgument[index + 1] - fArgument[index]; 355 if(!(delta != 0.0)) 355 if(!(delta != 0.0)) 356 { 356 { 357 G4Exception("G4DataInterpolation::FastCubi 357 G4Exception("G4DataInterpolation::FastCubicSpline()", "Error", 358 FatalException, "Bad fArgument 358 FatalException, "Bad fArgument input !"); 359 } 359 } 360 G4double a = (fArgument[index + 1] - pX) / d 360 G4double a = (fArgument[index + 1] - pX) / delta; 361 G4double b = (pX - fArgument[index]) / delta 361 G4double b = (pX - fArgument[index]) / delta; 362 362 363 // Final evaluation of cubic spline polynomi 363 // Final evaluation of cubic spline polynomial for return 364 364 365 return a * fFunction[index] + b * fFunction[ 365 return a * fFunction[index] + b * fFunction[index + 1] + 366 ((a * a * a - a) * fSecondDerivative[ 366 ((a * a * a - a) * fSecondDerivative[index] + 367 (b * b * b - b) * fSecondDerivative[ 367 (b * b * b - b) * fSecondDerivative[index + 1]) * 368 delta * delta / 6.0; 368 delta * delta / 6.0; 369 } 369 } 370 370 371 ////////////////////////////////////////////// 371 //////////////////////////////////////////////////////////////////////////// 372 // 372 // 373 // Given argument pX, returns index k, so that 373 // Given argument pX, returns index k, so that pX bracketed by fArgument[k] 374 // and fArgument[k+1] 374 // and fArgument[k+1] 375 375 376 G4int G4DataInterpolation::LocateArgument(G4do 376 G4int G4DataInterpolation::LocateArgument(G4double pX) const 377 { 377 { 378 G4int kLow = -1, kHigh = fNumber, k = 0; 378 G4int kLow = -1, kHigh = fNumber, k = 0; 379 G4bool ascend = (fArgument[fNumber - 1] >= f 379 G4bool ascend = (fArgument[fNumber - 1] >= fArgument[0]); 380 while((kHigh - kLow) > 1) 380 while((kHigh - kLow) > 1) 381 { 381 { 382 k = (kHigh + kLow) >> 1; // compute midpo 382 k = (kHigh + kLow) >> 1; // compute midpoint 'bisection' 383 if((pX >= fArgument[k]) == ascend) 383 if((pX >= fArgument[k]) == ascend) 384 { 384 { 385 kLow = k; 385 kLow = k; 386 } 386 } 387 else 387 else 388 { 388 { 389 kHigh = k; 389 kHigh = k; 390 } 390 } 391 } 391 } 392 if(!(pX != fArgument[0])) 392 if(!(pX != fArgument[0])) 393 { 393 { 394 return 1; 394 return 1; 395 } 395 } 396 if(!(pX != fArgument[fNumber - 1])) << 396 else if(!(pX != fArgument[fNumber - 1])) 397 { 397 { 398 return fNumber - 2; 398 return fNumber - 2; 399 } 399 } 400 << 400 else 401 return kLow; << 401 return kLow; 402 << 403 } 402 } 404 403 405 ////////////////////////////////////////////// 404 ///////////////////////////////////////////////////////////////////////////// 406 // 405 // 407 // Given a value pX, returns a value 'index' s 406 // Given a value pX, returns a value 'index' such that pX is between 408 // fArgument[index] and fArgument[index+1]. fA 407 // fArgument[index] and fArgument[index+1]. fArgument MUST BE MONOTONIC, 409 // either increasing or decreasing. If index = 408 // either increasing or decreasing. If index = -1 or fNumber, this indicates 410 // that pX is out of range. The value index on 409 // that pX is out of range. The value index on input is taken as the initial 411 // approximation for index on output. 410 // approximation for index on output. 412 411 413 void G4DataInterpolation::CorrelatedSearch(G4d 412 void G4DataInterpolation::CorrelatedSearch(G4double pX, G4int& index) const 414 { 413 { 415 G4int kHigh = 0, k = 0, Increment = 0; 414 G4int kHigh = 0, k = 0, Increment = 0; 416 // ascend = true for ascending order of tabl 415 // ascend = true for ascending order of table, false otherwise 417 G4bool ascend = (fArgument[fNumber - 1] >= f 416 G4bool ascend = (fArgument[fNumber - 1] >= fArgument[0]); 418 if(index < 0 || index > fNumber - 1) 417 if(index < 0 || index > fNumber - 1) 419 { 418 { 420 index = -1; 419 index = -1; 421 kHigh = fNumber; 420 kHigh = fNumber; 422 } 421 } 423 else 422 else 424 { 423 { 425 Increment = 1; // What value would be t 424 Increment = 1; // What value would be the best ? 426 if((pX >= fArgument[index]) == ascend) 425 if((pX >= fArgument[index]) == ascend) 427 { 426 { 428 if(index == fNumber - 1) 427 if(index == fNumber - 1) 429 { 428 { 430 index = fNumber; 429 index = fNumber; 431 return; 430 return; 432 } 431 } 433 kHigh = index + 1; 432 kHigh = index + 1; 434 while((pX >= fArgument[kHigh]) == ascend 433 while((pX >= fArgument[kHigh]) == ascend) 435 { 434 { 436 index = kHigh; 435 index = kHigh; 437 Increment += Increment; // double the 436 Increment += Increment; // double the Increment 438 kHigh = index + Increment; 437 kHigh = index + Increment; 439 if(kHigh > (fNumber - 1)) 438 if(kHigh > (fNumber - 1)) 440 { 439 { 441 kHigh = fNumber; 440 kHigh = fNumber; 442 break; 441 break; 443 } 442 } 444 } 443 } 445 } 444 } 446 else 445 else 447 { 446 { 448 if(index == 0) 447 if(index == 0) 449 { 448 { 450 index = -1; 449 index = -1; 451 return; 450 return; 452 } 451 } 453 kHigh = --index; 452 kHigh = --index; 454 while((pX < fArgument[index]) == ascend) 453 while((pX < fArgument[index]) == ascend) 455 { 454 { 456 kHigh = index; 455 kHigh = index; 457 Increment <<= 1; // double the Increm 456 Increment <<= 1; // double the Increment 458 if(Increment >= kHigh) 457 if(Increment >= kHigh) 459 { 458 { 460 index = -1; 459 index = -1; 461 break; 460 break; 462 } 461 } 463 << 462 else 464 index = kHigh - Increment; << 463 { 465 << 464 index = kHigh - Increment; >> 465 } 466 } 466 } 467 } // Value bracketed 467 } // Value bracketed 468 } 468 } 469 // final bisection searching 469 // final bisection searching 470 470 471 while((kHigh - index) != 1) 471 while((kHigh - index) != 1) 472 { 472 { 473 k = (kHigh + index) >> 1; 473 k = (kHigh + index) >> 1; 474 if((pX >= fArgument[k]) == ascend) 474 if((pX >= fArgument[k]) == ascend) 475 { 475 { 476 index = k; 476 index = k; 477 } 477 } 478 else 478 else 479 { 479 { 480 kHigh = k; 480 kHigh = k; 481 } 481 } 482 } 482 } 483 if(!(pX != fArgument[fNumber - 1])) 483 if(!(pX != fArgument[fNumber - 1])) 484 { 484 { 485 index = fNumber - 2; 485 index = fNumber - 2; 486 } 486 } 487 if(!(pX != fArgument[0])) 487 if(!(pX != fArgument[0])) 488 { 488 { 489 index = 0; 489 index = 0; 490 } 490 } 491 return; 491 return; 492 } 492 } 493 493 494 // 494 // 495 // 495 // 496 ////////////////////////////////////////////// 496 //////////////////////////////////////////////////////////////////////////// 497 497