00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "G4DataInterpolation.hh"
00030
00032
00033
00034
00035
00036 G4DataInterpolation::G4DataInterpolation( G4double pX[],
00037 G4double pY[],
00038 G4int number )
00039 : fArgument(new G4double[number]),
00040 fFunction(new G4double[number]),
00041 fSecondDerivative(0),
00042 fNumber(number)
00043 {
00044 for(G4int i=0;i<fNumber;i++)
00045 {
00046 fArgument[i] = pX[i] ;
00047 fFunction[i] = pY[i] ;
00048 }
00049 }
00050
00052
00053
00054
00055
00056
00057
00058 G4DataInterpolation::G4DataInterpolation( G4double pX[],
00059 G4double pY[],
00060 G4int number,
00061 G4double pFirstDerStart,
00062 G4double pFirstDerFinish )
00063 : fArgument(new G4double[number]),
00064 fFunction(new G4double[number]),
00065 fSecondDerivative(new G4double[number]),
00066 fNumber(number)
00067 {
00068 G4int i=0 ;
00069 G4double p=0.0, qn=0.0, sig=0.0, un=0.0 ;
00070 const G4double maxDerivative = 0.99e30 ;
00071 G4double* u = new G4double[fNumber - 1] ;
00072
00073 for(i=0;i<fNumber;i++)
00074 {
00075 fArgument[i] = pX[i] ;
00076 fFunction[i] = pY[i] ;
00077 }
00078 if(pFirstDerStart > maxDerivative)
00079 {
00080 fSecondDerivative[0] = 0.0 ;
00081 u[0] = 0.0 ;
00082 }
00083 else
00084 {
00085 fSecondDerivative[0] = -0.5 ;
00086 u[0] = (3.0/(fArgument[1]-fArgument[0]))
00087 * ((fFunction[1]-fFunction[0])/(fArgument[1]-fArgument[0])
00088 - pFirstDerStart) ;
00089 }
00090
00091
00092
00093
00094 for(i=1;i<fNumber-1;i++)
00095 {
00096 sig = (fArgument[i]-fArgument[i-1])/(fArgument[i+1]-fArgument[i-1]) ;
00097 p = sig*fSecondDerivative[i-1] + 2.0 ;
00098 fSecondDerivative[i] = (sig - 1.0)/p ;
00099 u[i] = (fFunction[i+1]-fFunction[i])/(fArgument[i+1]-fArgument[i]) -
00100 (fFunction[i]-fFunction[i-1])/(fArgument[i]-fArgument[i-1]) ;
00101 u[i] =(6.0*u[i]/(fArgument[i+1]-fArgument[i-1]) - sig*u[i-1])/p ;
00102 }
00103 if(pFirstDerFinish > maxDerivative)
00104 {
00105 qn = 0.0 ;
00106 un = 0.0 ;
00107 }
00108 else
00109 {
00110 qn = 0.5 ;
00111 un = (3.0/(fArgument[fNumber-1]-fArgument[fNumber-2]))
00112 * (pFirstDerFinish - (fFunction[fNumber-1]-fFunction[fNumber-2])
00113 / (fArgument[fNumber-1]-fArgument[fNumber-2])) ;
00114 }
00115 fSecondDerivative[fNumber-1] = (un - qn*u[fNumber-2])/
00116 (qn*fSecondDerivative[fNumber-2] + 1.0) ;
00117
00118
00119
00120
00121 for(G4int k=fNumber-2;k>=0;k--)
00122 {
00123 fSecondDerivative[k] = fSecondDerivative[k]*fSecondDerivative[k+1] + u[k];
00124 }
00125 delete[] u ;
00126 }
00127
00129
00130
00131
00132
00133 G4DataInterpolation::~G4DataInterpolation()
00134 {
00135 delete [] fArgument ;
00136 delete [] fFunction ;
00137 if(fSecondDerivative) { delete [] fSecondDerivative; }
00138 }
00139
00141
00142
00143
00144
00145
00146
00147 G4double
00148 G4DataInterpolation::PolynomInterpolation(G4double pX,
00149 G4double& deltaY ) const
00150 {
00151 G4int i=0, j=1, k=0 ;
00152 G4double mult=0.0, difi=0.0, deltaLow=0.0, deltaUp=0.0, cd=0.0, y=0.0 ;
00153 G4double* c = new G4double[fNumber] ;
00154 G4double* d = new G4double[fNumber] ;
00155 G4double diff = std::fabs(pX-fArgument[0]) ;
00156 for(i=0;i<fNumber;i++)
00157 {
00158 difi = std::fabs(pX-fArgument[i]) ;
00159 if(difi <diff)
00160 {
00161 k = i ;
00162 diff = difi ;
00163 }
00164 c[i] = fFunction[i] ;
00165 d[i] = fFunction[i] ;
00166 }
00167 y = fFunction[k--] ;
00168 for(j=1;j<fNumber;j++)
00169 {
00170 for(i=0;i<fNumber-j;i++)
00171 {
00172 deltaLow = fArgument[i] - pX ;
00173 deltaUp = fArgument[i+j] - pX ;
00174 cd = c[i+1] - d[i] ;
00175 mult = deltaLow - deltaUp ;
00176 if (!(mult != 0.0))
00177 {
00178 G4Exception("G4DataInterpolation::PolynomInterpolation()",
00179 "Error", FatalException, "Coincident nodes !") ;
00180 }
00181 mult = cd/mult ;
00182 d[i] = deltaUp*mult ;
00183 c[i] = deltaLow*mult ;
00184 }
00185 y += (deltaY = (2*k < (fNumber - j -1) ? c[k+1] : d[k--] )) ;
00186 }
00187 delete[] c ;
00188 delete[] d ;
00189
00190 return y ;
00191 }
00192
00194
00195
00196
00197
00198
00199
00200
00201 void
00202 G4DataInterpolation::PolIntCoefficient( G4double cof[]) const
00203 {
00204 G4int i=0, j=0 ;
00205 G4double factor;
00206 G4double reducedY=0.0, mult=1.0 ;
00207 G4double* tempArgument = new G4double[fNumber] ;
00208
00209 for(i=0;i<fNumber;i++)
00210 {
00211 tempArgument[i] = cof[i] = 0.0 ;
00212 }
00213 tempArgument[fNumber-1] = -fArgument[0] ;
00214
00215 for(i=1;i<fNumber;i++)
00216 {
00217 for(j=fNumber-1-i;j<fNumber-1;j++)
00218 {
00219 tempArgument[j] -= fArgument[i]*tempArgument[j+1] ;
00220 }
00221 tempArgument[fNumber-1] -= fArgument[i] ;
00222 }
00223 for(i=0;i<fNumber;i++)
00224 {
00225 factor = fNumber ;
00226 for(j=fNumber-1;j>=1;j--)
00227 {
00228 factor = j*tempArgument[j] + factor*fArgument[i] ;
00229 }
00230 reducedY = fFunction[i]/factor ;
00231 mult = 1.0 ;
00232 for(j=fNumber-1;j>=0;j--)
00233 {
00234 cof[j] += mult*reducedY ;
00235 mult = tempArgument[j] + mult*fArgument[i] ;
00236 }
00237 }
00238 delete[] tempArgument ;
00239 }
00240
00242
00243
00244
00245
00246
00247
00248 G4double
00249 G4DataInterpolation::RationalPolInterpolation(G4double pX,
00250 G4double& deltaY ) const
00251 {
00252 G4int i=0, j=1, k=0 ;
00253 const G4double tolerance = 1.6e-24 ;
00254 G4double mult=0.0, difi=0.0, cd=0.0, y=0.0, cof=0.0 ;
00255 G4double* c = new G4double[fNumber] ;
00256 G4double* d = new G4double[fNumber] ;
00257 G4double diff = std::fabs(pX-fArgument[0]) ;
00258 for(i=0;i<fNumber;i++)
00259 {
00260 difi = std::fabs(pX-fArgument[i]) ;
00261 if (!(difi != 0.0))
00262 {
00263 y = fFunction[i] ;
00264 deltaY = 0.0 ;
00265 delete[] c ;
00266 delete[] d ;
00267 return y ;
00268 }
00269 else if(difi < diff)
00270 {
00271 k = i ;
00272 diff = difi ;
00273 }
00274 c[i] = fFunction[i] ;
00275 d[i] = fFunction[i] + tolerance ;
00276 }
00277 y = fFunction[k--] ;
00278 for(j=1;j<fNumber;j++)
00279 {
00280 for(i=0;i<fNumber-j;i++)
00281 {
00282 cd = c[i+1] - d[i] ;
00283 difi = fArgument[i+j] - pX ;
00284 cof = (fArgument[i] - pX)*d[i]/difi ;
00285 mult = cof - c[i+1] ;
00286 if (!(mult != 0.0))
00287 {
00288 G4Exception("G4DataInterpolation::RationalPolInterpolation()",
00289 "Error", FatalException, "Coincident nodes !") ;
00290 }
00291 mult = cd/mult ;
00292 d[i] = c[i+1]*mult ;
00293 c[i] = cof*mult ;
00294 }
00295 y += (deltaY = (2*k < (fNumber - j - 1) ? c[k+1] : d[k--] )) ;
00296 }
00297 delete[] c ;
00298 delete[] d ;
00299
00300 return y ;
00301 }
00302
00304
00305
00306
00307
00308
00309
00310 G4double
00311 G4DataInterpolation::CubicSplineInterpolation(G4double pX) const
00312 {
00313 G4int kLow=0, kHigh=fNumber-1, k=0 ;
00314
00315
00316
00317
00318 while((kHigh - kLow) > 1)
00319 {
00320 k = (kHigh + kLow) >> 1 ;
00321 if(fArgument[k] > pX)
00322 {
00323 kHigh = k ;
00324 }
00325 else
00326 {
00327 kLow = k ;
00328 }
00329 }
00330 G4double deltaHL = fArgument[kHigh] - fArgument[kLow] ;
00331 if (!(deltaHL != 0.0))
00332 {
00333 G4Exception("G4DataInterpolation::CubicSplineInterpolation()",
00334 "Error", FatalException, "Bad fArgument input !") ;
00335 }
00336 G4double a = (fArgument[kHigh] - pX)/deltaHL ;
00337 G4double b = (pX - fArgument[kLow])/deltaHL ;
00338
00339
00340
00341 return a*fFunction[kLow] + b*fFunction[kHigh] +
00342 ((a*a*a - a)*fSecondDerivative[kLow] +
00343 (b*b*b - b)*fSecondDerivative[kHigh])*deltaHL*deltaHL/6.0 ;
00344 }
00345
00347
00348
00349
00350
00351
00352 G4double
00353 G4DataInterpolation::FastCubicSpline(G4double pX,
00354 G4int index) const
00355 {
00356 G4double delta = fArgument[index+1] - fArgument[index] ;
00357 if (!(delta != 0.0))
00358 {
00359 G4Exception("G4DataInterpolation::FastCubicSpline()",
00360 "Error", FatalException, "Bad fArgument input !") ;
00361 }
00362 G4double a = (fArgument[index+1] - pX)/delta ;
00363 G4double b = (pX - fArgument[index])/delta ;
00364
00365
00366
00367 return a*fFunction[index] + b*fFunction[index+1] +
00368 ((a*a*a - a)*fSecondDerivative[index] +
00369 (b*b*b - b)*fSecondDerivative[index+1])*delta*delta/6.0 ;
00370 }
00371
00373
00374
00375
00376
00377 G4int
00378 G4DataInterpolation::LocateArgument(G4double pX) const
00379 {
00380 G4int kLow=-1, kHigh=fNumber, k=0 ;
00381 G4bool ascend=(fArgument[fNumber-1] >= fArgument[0]) ;
00382 while((kHigh - kLow) > 1)
00383 {
00384 k = (kHigh + kLow) >> 1 ;
00385 if( (pX >= fArgument[k]) == ascend)
00386 {
00387 kLow = k ;
00388 }
00389 else
00390 {
00391 kHigh = k ;
00392 }
00393 }
00394 if (!(pX != fArgument[0]))
00395 {
00396 return 1 ;
00397 }
00398 else if (!(pX != fArgument[fNumber-1]))
00399 {
00400 return fNumber - 2 ;
00401 }
00402 else return kLow ;
00403 }
00404
00406
00407
00408
00409
00410
00411
00412
00413 void
00414 G4DataInterpolation::CorrelatedSearch( G4double pX,
00415 G4int& index ) const
00416 {
00417 G4int kHigh=0, k=0, Increment=0 ;
00418
00419 G4bool ascend = (fArgument[fNumber-1] >= fArgument[0]) ;
00420 if(index < 0 || index > fNumber-1)
00421 {
00422 index = -1 ;
00423 kHigh = fNumber ;
00424 }
00425 else
00426 {
00427 Increment = 1 ;
00428 if((pX >= fArgument[index]) == ascend)
00429 {
00430 if(index == fNumber -1)
00431 {
00432 index = fNumber ;
00433 return ;
00434 }
00435 kHigh = index + 1 ;
00436 while((pX >= fArgument[kHigh]) == ascend)
00437 {
00438 index = kHigh ;
00439 Increment += Increment ;
00440 kHigh = index + Increment ;
00441 if(kHigh > (fNumber - 1))
00442 {
00443 kHigh = fNumber ;
00444 break ;
00445 }
00446 }
00447 }
00448 else
00449 {
00450 if(index == 0)
00451 {
00452 index = -1 ;
00453 return ;
00454 }
00455 kHigh = index-- ;
00456 while((pX < fArgument[index]) == ascend)
00457 {
00458 kHigh = index ;
00459 Increment <<= 1 ;
00460 if(Increment >= kHigh)
00461 {
00462 index = -1 ;
00463 break ;
00464 }
00465 else
00466 {
00467 index = kHigh - Increment ;
00468 }
00469 }
00470 }
00471 }
00472
00473
00474 while((kHigh - index) != 1)
00475 {
00476 k = (kHigh + index) >> 1 ;
00477 if((pX >= fArgument[k]) == ascend)
00478 {
00479 index = k ;
00480 }
00481 else
00482 {
00483 kHigh = k ;
00484 }
00485 }
00486 if (!(pX != fArgument[fNumber-1]))
00487 {
00488 index = fNumber - 2 ;
00489 }
00490 if (!(pX != fArgument[0]))
00491 {
00492 index = 0 ;
00493 }
00494 return ;
00495 }
00496
00497
00498