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 11.0.p3,)


** Warning: Cannot open xref database.

  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 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,     
 36 // data members                                   
 37                                                   
 38 G4DataInterpolation::G4DataInterpolation(G4dou    
 39                                          G4int    
 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.    
 54 // fSecondDerivative[0,...fNumber-1] which is     
 55 // the function                                   
 56                                                   
 57 G4DataInterpolation::G4DataInterpolation(G4dou    
 58                                          G4int    
 59                                          G4dou    
 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 =     
 67   const G4double maxDerivative = 0.99e30;         
 68   auto u                  = new G4double[fNumb    
 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    
 84            ((fFunction[1] - fFunction[0]) / (f    
 85             pFirstDerStart);                      
 86   }                                               
 87                                                   
 88   // Decomposition loop for tridiagonal algori    
 89   // and u[i] are used for temporary storage o    
 90                                                   
 91   for(i = 1; i < fNumber - 1; ++i)                
 92   {                                               
 93     sig =                                         
 94       (fArgument[i] - fArgument[i - 1]) / (fAr    
 95     p                    = sig * fSecondDeriva    
 96     fSecondDerivative[i] = (sig - 1.0) / p;       
 97     u[i] =                                        
 98       (fFunction[i + 1] - fFunction[i]) / (fAr    
 99       (fFunction[i] - fFunction[i - 1]) / (fAr    
100     u[i] =                                        
101       (6.0 * u[i] / (fArgument[i + 1] - fArgum    
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] - fArgume    
113       (pFirstDerFinish - (fFunction[fNumber -     
114                            (fArgument[fNumber     
115   }                                               
116   fSecondDerivative[fNumber - 1] =                
117     (un - qn * u[fNumber - 2]) / (qn * fSecond    
118                                                   
119   // The backsubstitution loop for the triagon    
120   // a linear system of equations.                
121                                                   
122   for(G4int k = fNumber - 2; k >= 0; --k)         
123   {                                               
124     fSecondDerivative[k] =                        
125       fSecondDerivative[k] * fSecondDerivative    
126   }                                               
127   delete[] u;                                     
128 }                                                 
129                                                   
130 //////////////////////////////////////////////    
131 //                                                
132 // Destructor deletes dynamically created arra    
133 // fFunction and fSecondDerivative, all have d    
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), wher    
145 // degree such that P(fArgument[i]) = fFunctio    
146 // This is Lagrange's form of interpolation an    
147 // algorithm                                      
148                                                   
149 G4double G4DataInterpolation::PolynomInterpola    
150                                                   
151 {                                                 
152   G4int i = 0, j = 1, k = 0;                      
153   G4double mult = 0.0, difi = 0.0, deltaLow =     
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::Poly    
181                     FatalException, "Coinciden    
182       }                                           
183       mult = cd / mult;                           
184       d[i] = deltaUp * mult;                      
185       c[i] = deltaLow * mult;                     
186     }                                             
187     y += (deltaY = (2 * k < (fNumber - j - 1)     
188   }                                               
189   delete[] c;                                     
190   delete[] d;                                     
191                                                   
192   return y;                                       
193 }                                                 
194                                                   
195 //////////////////////////////////////////////    
196 //                                                
197 // Given arrays fArgument[0,..,fNumber-1] and     
198 // function calculates an array of coefficient    
199 // usually (fNumber>10) better accuracy for po    
200 // with PolynomInterpolation function. They co    
201 // calculations and some other applications.      
202                                                   
203 void G4DataInterpolation::PolIntCoefficient(G4    
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;     
219     {                                             
220       tempArgument[j] -= fArgument[i] * tempAr    
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 *     
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 * fArgumen    
237     }                                             
238   }                                               
239   delete[] tempArgument;                          
240 }                                                 
241                                                   
242 //////////////////////////////////////////////    
243 //                                                
244 // The function returns diagonal rational func    
245 // algorithm of Neville type) Pn(x)/Qm(x) wher    
246 // Tests showed the method is not stable and h    
247 // with polynomial interpolation ?!               
248                                                   
249 G4double G4DataInterpolation::RationalPolInter    
250                                                   
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    
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 pr    
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 in    
287       {                                           
288         G4Exception("G4DataInterpolation::Rati    
289                     FatalException, "Coinciden    
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)     
296   }                                               
297   delete[] c;                                     
298   delete[] d;                                     
299                                                   
300   return y;                                       
301 }                                                 
302                                                   
303 //////////////////////////////////////////////    
304 //                                                
305 // Cubic spline interpolation in point pX for     
306 // fArgument, fFunction. The constructor, whic    
307 // must be called before. The function works o    
308 // are in random values of pX.                    
309                                                   
310 G4double G4DataInterpolation::CubicSplineInter    
311 {                                                 
312   G4int kLow = 0, kHigh = fNumber - 1, k = 0;     
313                                                   
314   // Searching in the table by means of bisect    
315   // fArgument must be monotonic, either incre    
316                                                   
317   while((kHigh - kLow) > 1)                       
318   {                                               
319     k = (kHigh + kLow) >> 1;  // compute midpo    
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 v    
329   G4double deltaHL = fArgument[kHigh] - fArgum    
330   if(!(deltaHL != 0.0))                           
331   {                                               
332     G4Exception("G4DataInterpolation::CubicSpl    
333                 FatalException, "Bad fArgument    
334   }                                               
335   G4double a = (fArgument[kHigh] - pX) / delta    
336   G4double b = (pX - fArgument[kLow]) / deltaH    
337                                                   
338   // Final evaluation of cubic spline polynomi    
339                                                   
340   return a * fFunction[kLow] + b * fFunction[k    
341          ((a * a * a - a) * fSecondDerivative[    
342           (b * b * b - b) * fSecondDerivative[    
343            deltaHL * deltaHL / 6.0;               
344 }                                                 
345                                                   
346 //////////////////////////////////////////////    
347 //                                                
348 // Return cubic spline interpolation in the po    
349 // fArgument[index] and fArgument[index+1]. It    
350 // of known from external analysis values of i    
351                                                   
352 G4double G4DataInterpolation::FastCubicSpline(    
353 {                                                 
354   G4double delta = fArgument[index + 1] - fArg    
355   if(!(delta != 0.0))                             
356   {                                               
357     G4Exception("G4DataInterpolation::FastCubi    
358                 FatalException, "Bad fArgument    
359   }                                               
360   G4double a = (fArgument[index + 1] - pX) / d    
361   G4double b = (pX - fArgument[index]) / delta    
362                                                   
363   // Final evaluation of cubic spline polynomi    
364                                                   
365   return a * fFunction[index] + b * fFunction[    
366          ((a * a * a - a) * fSecondDerivative[    
367           (b * b * b - b) * fSecondDerivative[    
368            delta * delta / 6.0;                   
369 }                                                 
370                                                   
371 //////////////////////////////////////////////    
372 //                                                
373 // Given argument pX, returns index k, so that    
374 // and fArgument[k+1]                             
375                                                   
376 G4int G4DataInterpolation::LocateArgument(G4do    
377 {                                                 
378   G4int kLow = -1, kHigh = fNumber, k = 0;        
379   G4bool ascend = (fArgument[fNumber - 1] >= f    
380   while((kHigh - kLow) > 1)                       
381   {                                               
382     k = (kHigh + kLow) >> 1;  // compute midpo    
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' s    
408 // fArgument[index] and fArgument[index+1]. fA    
409 // either increasing or decreasing. If index =    
410 // that pX is out of range. The value index on    
411 // approximation for index on output.             
412                                                   
413 void G4DataInterpolation::CorrelatedSearch(G4d    
414 {                                                 
415   G4int kHigh = 0, k = 0, Increment = 0;          
416   // ascend = true for ascending order of tabl    
417   G4bool ascend = (fArgument[fNumber - 1] >= f    
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 t    
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    
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 Increm    
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