Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/HEPNumerics/src/G4DataInterpolation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /global/HEPNumerics/src/G4DataInterpolation.cc (Version 11.3.0) and /global/HEPNumerics/src/G4DataInterpolation.cc (Version 8.3.p2)


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