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