Geant4: G4DataInterpolation Class Reference

G4DataInterpolation::G4DataInterpolation	(	G4double	pX[],
		G4double	pY[],
		G4int	number
	)

Definition at line 36 of file G4DataInterpolation.cc.

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 }

G4DataInterpolation::G4DataInterpolation	(	G4double	pX[],
		G4double	pY[],
		G4int	number,
		G4double	pFirstDerStart,
		G4double	pFirstDerFinish
	)

Definition at line 58 of file G4DataInterpolation.cc.

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    // Decomposition loop for tridiagonal algorithm. fSecondDerivative[i]
00092    // and u[i] are used for temporary storage of the decomposed factors.
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    // The backsubstitution loop for the triagonal algorithm of solving
00119    // a linear system of equations.
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 }

G4DataInterpolation::~G4DataInterpolation ( )

Definition at line 133 of file G4DataInterpolation.cc.

00134 {
00135    delete [] fArgument ;
00136    delete [] fFunction ;
00137    if(fSecondDerivative)  { delete [] fSecondDerivative; }
00138 }

void G4DataInterpolation::CorrelatedSearch	(	G4double	pX,
		G4int &	index
	)			const

Definition at line 414 of file G4DataInterpolation.cc.

00416 {
00417    G4int kHigh=0, k=0, Increment=0 ;
00418    // ascend = true for ascending order of table, false otherwise
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 ;   //   What value would be the best ?
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 ;      // double the 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 ;      // double the Increment
00460             if(Increment >= kHigh)
00461             {
00462                index = -1 ;
00463                break ;
00464             }
00465             else
00466             {
00467                index = kHigh - Increment ;
00468             }
00469          }
00470       }        // Value bracketed
00471    }
00472    // final bisection searching
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 }

G4double G4DataInterpolation::CubicSplineInterpolation ( G4double pX ) const

Definition at line 311 of file G4DataInterpolation.cc.

References FatalException, and G4Exception().

00312 {
00313    G4int kLow=0, kHigh=fNumber-1, k=0 ;
00314    
00315    // Searching in the table by means of bisection method. 
00316    // fArgument must be monotonic, either increasing or decreasing
00317    
00318    while((kHigh - kLow) > 1)
00319    {
00320       k = (kHigh + kLow) >> 1 ; // compute midpoint 'bisection'
00321       if(fArgument[k] > pX)
00322       {
00323          kHigh = k ;
00324       }
00325       else
00326       {
00327          kLow = k ;
00328       }
00329    }        // kLow and kHigh now bracket the input value of pX
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    // Final evaluation of cubic spline polynomial for return
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 }

G4double G4DataInterpolation::FastCubicSpline	(	G4double	pX,
		G4int	index
	)			const

Definition at line 353 of file G4DataInterpolation.cc.

References FatalException, and G4Exception().

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    // Final evaluation of cubic spline polynomial for return
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 }

G4int G4DataInterpolation::LocateArgument ( G4double pX ) const

Definition at line 378 of file G4DataInterpolation.cc.

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 ; // compute midpoint 'bisection'
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 }

void G4DataInterpolation::PolIntCoefficient ( G4double cof[] ) const

Definition at line 202 of file G4DataInterpolation.cc.

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 }

G4double G4DataInterpolation::PolynomInterpolation	(	G4double	pX,
		G4double &	deltaY
	)			const

Definition at line 148 of file G4DataInterpolation.cc.

References FatalException, and G4Exception().

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 }

G4double G4DataInterpolation::RationalPolInterpolation	(	G4double	pX,
		G4double &	deltaY
	)			const

Definition at line 249 of file G4DataInterpolation.cc.

References FatalException, and G4Exception().

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 ;   // to prevent rare zero/zero cases
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))    // function to be interpolated has pole at pX
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 }


Public Member Functions
	G4DataInterpolation (G4double pX[], G4double pY[], G4int number)
	G4DataInterpolation (G4double pX[], G4double pY[], G4int number, G4double pFirstDerStart, G4double pFirstDerFinish)
	~G4DataInterpolation ()
G4double	PolynomInterpolation (G4double pX, G4double &deltaY) const
void	PolIntCoefficient (G4double cof[]) const
G4double	RationalPolInterpolation (G4double pX, G4double &deltaY) const
G4double	CubicSplineInterpolation (G4double pX) const
G4double	FastCubicSpline (G4double pX, G4int index) const
G4int	LocateArgument (G4double pX) const
void	CorrelatedSearch (G4double pX, G4int &index) const

G4DataInterpolation Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation