Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/HEPNumerics/src/G4ChebyshevApproximation.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/G4ChebyshevApproximation.cc (Version 11.3.0) and /global/HEPNumerics/src/G4ChebyshevApproximation.cc (Version 10.7.p4)


  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 // G4ChebyshevApproximation class implementati     26 // G4ChebyshevApproximation class implementation
 27 //                                                 27 //
 28 // Author: V.Grichine, 24.04.97                    28 // Author: V.Grichine, 24.04.97
 29 // -------------------------------------------     29 // --------------------------------------------------------------------
 30                                                    30 
 31 #include "G4ChebyshevApproximation.hh"             31 #include "G4ChebyshevApproximation.hh"
 32 #include "G4PhysicalConstants.hh"                  32 #include "G4PhysicalConstants.hh"
 33                                                    33 
 34 // Constructor for initialisation of the class     34 // Constructor for initialisation of the class data members.
 35 // It creates the array fChebyshevCof[0,...,fN     35 // It creates the array fChebyshevCof[0,...,fNumber-1], fNumber = n ;
 36 // which consists of Chebyshev coefficients de     36 // which consists of Chebyshev coefficients describing the function
 37 // pointed by pFunction. The values a and b fi     37 // pointed by pFunction. The values a and b fix the interval of validity
 38 // of the Chebyshev approximation.                 38 // of the Chebyshev approximation.
 39                                                    39 
 40 G4ChebyshevApproximation::G4ChebyshevApproxima     40 G4ChebyshevApproximation::G4ChebyshevApproximation(function pFunction, G4int n,
 41                                                    41                                                    G4double a, G4double b)
 42   : fFunction(pFunction)                           42   : fFunction(pFunction)
 43   , fNumber(n)                                     43   , fNumber(n)
 44   , fChebyshevCof(new G4double[fNumber])           44   , fChebyshevCof(new G4double[fNumber])
 45   , fMean(0.5 * (b + a))                           45   , fMean(0.5 * (b + a))
 46   , fDiff(0.5 * (b - a))                           46   , fDiff(0.5 * (b - a))
 47 {                                                  47 {
 48   G4int i = 0, j = 0;                              48   G4int i = 0, j = 0;
 49   G4double rootSum = 0.0, cofj = 0.0;              49   G4double rootSum = 0.0, cofj = 0.0;
 50   auto tempFunction = new G4double[fNumber];   <<  50   G4double* tempFunction = new G4double[fNumber];
 51   G4double weight        = 2.0 / fNumber;          51   G4double weight        = 2.0 / fNumber;
 52   G4double cof           = 0.5 * weight * pi;      52   G4double cof           = 0.5 * weight * pi;  // pi/n
 53                                                    53 
 54   for(i = 0; i < fNumber; ++i)                     54   for(i = 0; i < fNumber; ++i)
 55   {                                                55   {
 56     rootSum         = std::cos(cof * (i + 0.5)     56     rootSum         = std::cos(cof * (i + 0.5));
 57     tempFunction[i] = fFunction(rootSum * fDif     57     tempFunction[i] = fFunction(rootSum * fDiff + fMean);
 58   }                                                58   }
 59   for(j = 0; j < fNumber; ++j)                     59   for(j = 0; j < fNumber; ++j)
 60   {                                                60   {
 61     cofj    = cof * j;                             61     cofj    = cof * j;
 62     rootSum = 0.0;                                 62     rootSum = 0.0;
 63                                                    63 
 64     for(i = 0; i < fNumber; ++i)                   64     for(i = 0; i < fNumber; ++i)
 65     {                                              65     {
 66       rootSum += tempFunction[i] * std::cos(co     66       rootSum += tempFunction[i] * std::cos(cofj * (i + 0.5));
 67     }                                              67     }
 68     fChebyshevCof[j] = weight * rootSum;           68     fChebyshevCof[j] = weight * rootSum;
 69   }                                                69   }
 70   delete[] tempFunction;                           70   delete[] tempFunction;
 71 }                                                  71 }
 72                                                    72 
 73 // -------------------------------------------     73 // --------------------------------------------------------------------
 74 //                                                 74 //
 75 // Constructor for creation of Chebyshev coeff     75 // Constructor for creation of Chebyshev coefficients for mx-derivative
 76 // from pFunction. The value of mx ! MUST BE !     76 // from pFunction. The value of mx ! MUST BE ! < nx , because the result
 77 // array of fChebyshevCof will be of (nx-mx) s     77 // array of fChebyshevCof will be of (nx-mx) size.  The values a and b
 78 // fix the interval of validity of the Chebysh     78 // fix the interval of validity of the Chebyshev approximation.
 79                                                    79 
 80 G4ChebyshevApproximation::G4ChebyshevApproxima     80 G4ChebyshevApproximation::G4ChebyshevApproximation(function pFunction, G4int nx,
 81                                                    81                                                    G4int mx, G4double a,
 82                                                    82                                                    G4double b)
 83   : fFunction(pFunction)                           83   : fFunction(pFunction)
 84   , fNumber(nx)                                    84   , fNumber(nx)
 85   , fChebyshevCof(new G4double[fNumber])           85   , fChebyshevCof(new G4double[fNumber])
 86   , fMean(0.5 * (b + a))                           86   , fMean(0.5 * (b + a))
 87   , fDiff(0.5 * (b - a))                           87   , fDiff(0.5 * (b - a))
 88 {                                                  88 {
 89   if(nx <= mx)                                     89   if(nx <= mx)
 90   {                                                90   {
 91     G4Exception("G4ChebyshevApproximation::G4C     91     G4Exception("G4ChebyshevApproximation::G4ChebyshevApproximation()",
 92                 "InvalidCall", FatalException,     92                 "InvalidCall", FatalException, "Invalid arguments !");
 93   }                                                93   }
 94   G4int i = 0, j = 0;                              94   G4int i = 0, j = 0;
 95   G4double rootSum = 0.0, cofj = 0.0;              95   G4double rootSum = 0.0, cofj = 0.0;
 96   auto tempFunction = new G4double[fNumber];   <<  96   G4double* tempFunction = new G4double[fNumber];
 97   G4double weight        = 2.0 / fNumber;          97   G4double weight        = 2.0 / fNumber;
 98   G4double cof           = 0.5 * weight * pi;      98   G4double cof           = 0.5 * weight * pi;  // pi/nx
 99                                                    99 
100   for(i = 0; i < fNumber; ++i)                    100   for(i = 0; i < fNumber; ++i)
101   {                                               101   {
102     rootSum         = std::cos(cof * (i + 0.5)    102     rootSum         = std::cos(cof * (i + 0.5));
103     tempFunction[i] = fFunction(rootSum * fDif    103     tempFunction[i] = fFunction(rootSum * fDiff + fMean);
104   }                                               104   }
105   for(j = 0; j < fNumber; ++j)                    105   for(j = 0; j < fNumber; ++j)
106   {                                               106   {
107     cofj    = cof * j;                            107     cofj    = cof * j;
108     rootSum = 0.0;                                108     rootSum = 0.0;
109                                                   109 
110     for(i = 0; i < fNumber; ++i)                  110     for(i = 0; i < fNumber; ++i)
111     {                                             111     {
112       rootSum += tempFunction[i] * std::cos(co    112       rootSum += tempFunction[i] * std::cos(cofj * (i + 0.5));
113     }                                             113     }
114     fChebyshevCof[j] = weight * rootSum;  // c    114     fChebyshevCof[j] = weight * rootSum;  // corresponds to pFunction
115   }                                               115   }
116   // Chebyshev coefficients for (mx)-derivativ    116   // Chebyshev coefficients for (mx)-derivative of pFunction
117                                                   117 
118   for(i = 1; i <= mx; ++i)                        118   for(i = 1; i <= mx; ++i)
119   {                                               119   {
120     DerivativeChebyshevCof(tempFunction);         120     DerivativeChebyshevCof(tempFunction);
121     fNumber--;                                    121     fNumber--;
122     for(j = 0; j < fNumber; ++j)                  122     for(j = 0; j < fNumber; ++j)
123     {                                             123     {
124       fChebyshevCof[j] = tempFunction[j];  //     124       fChebyshevCof[j] = tempFunction[j];  // corresponds to (i)-derivative
125     }                                             125     }
126   }                                               126   }
127   delete[] tempFunction;  // delete of dynamic    127   delete[] tempFunction;  // delete of dynamically allocated tempFunction
128 }                                                 128 }
129                                                   129 
130 // -------------------------------------------    130 // ------------------------------------------------------
131 //                                                131 //
132 // Constructor for creation of Chebyshev coeff    132 // Constructor for creation of Chebyshev coefficients for integral
133 // from pFunction.                                133 // from pFunction.
134                                                   134 
135 G4ChebyshevApproximation::G4ChebyshevApproxima    135 G4ChebyshevApproximation::G4ChebyshevApproximation(function pFunction,
136                                                   136                                                    G4double a, G4double b,
137                                                   137                                                    G4int n)
138   : fFunction(pFunction)                          138   : fFunction(pFunction)
139   , fNumber(n)                                    139   , fNumber(n)
140   , fChebyshevCof(new G4double[fNumber])          140   , fChebyshevCof(new G4double[fNumber])
141   , fMean(0.5 * (b + a))                          141   , fMean(0.5 * (b + a))
142   , fDiff(0.5 * (b - a))                          142   , fDiff(0.5 * (b - a))
143 {                                                 143 {
144   G4int i = 0, j = 0;                             144   G4int i = 0, j = 0;
145   G4double rootSum = 0.0, cofj = 0.0;             145   G4double rootSum = 0.0, cofj = 0.0;
146   auto tempFunction = new G4double[fNumber];   << 146   G4double* tempFunction = new G4double[fNumber];
147   G4double weight        = 2.0 / fNumber;         147   G4double weight        = 2.0 / fNumber;
148   G4double cof           = 0.5 * weight * pi;     148   G4double cof           = 0.5 * weight * pi;  // pi/n
149                                                   149 
150   for(i = 0; i < fNumber; ++i)                    150   for(i = 0; i < fNumber; ++i)
151   {                                               151   {
152     rootSum         = std::cos(cof * (i + 0.5)    152     rootSum         = std::cos(cof * (i + 0.5));
153     tempFunction[i] = fFunction(rootSum * fDif    153     tempFunction[i] = fFunction(rootSum * fDiff + fMean);
154   }                                               154   }
155   for(j = 0; j < fNumber; ++j)                    155   for(j = 0; j < fNumber; ++j)
156   {                                               156   {
157     cofj    = cof * j;                            157     cofj    = cof * j;
158     rootSum = 0.0;                                158     rootSum = 0.0;
159                                                   159 
160     for(i = 0; i < fNumber; ++i)                  160     for(i = 0; i < fNumber; ++i)
161     {                                             161     {
162       rootSum += tempFunction[i] * std::cos(co    162       rootSum += tempFunction[i] * std::cos(cofj * (i + 0.5));
163     }                                             163     }
164     fChebyshevCof[j] = weight * rootSum;  // c    164     fChebyshevCof[j] = weight * rootSum;  // corresponds to pFunction
165   }                                               165   }
166   // Chebyshev coefficients for integral of pF    166   // Chebyshev coefficients for integral of pFunction
167                                                   167 
168   IntegralChebyshevCof(tempFunction);             168   IntegralChebyshevCof(tempFunction);
169   for(j = 0; j < fNumber; ++j)                    169   for(j = 0; j < fNumber; ++j)
170   {                                               170   {
171     fChebyshevCof[j] = tempFunction[j];  // co    171     fChebyshevCof[j] = tempFunction[j];  // corresponds to integral
172   }                                               172   }
173   delete[] tempFunction;  // delete of dynamic    173   delete[] tempFunction;  // delete of dynamically allocated tempFunction
174 }                                                 174 }
175                                                   175 
176 // -------------------------------------------    176 // ---------------------------------------------------------------
177 //                                                177 //
178 // Destructor deletes the array of Chebyshev c    178 // Destructor deletes the array of Chebyshev coefficients
179                                                   179 
180 G4ChebyshevApproximation::~G4ChebyshevApproxim    180 G4ChebyshevApproximation::~G4ChebyshevApproximation()
181 {                                                 181 {
182   delete[] fChebyshevCof;                         182   delete[] fChebyshevCof;
183 }                                                 183 }
184                                                   184 
185 // -------------------------------------------    185 // ---------------------------------------------------------------
186 //                                                186 //
187 // Access function for Chebyshev coefficients     187 // Access function for Chebyshev coefficients
188 //                                                188 //
189                                                   189 
190 G4double G4ChebyshevApproximation::GetChebyshe    190 G4double G4ChebyshevApproximation::GetChebyshevCof(G4int number) const
191 {                                                 191 {
192   if(number < 0 && number >= fNumber)             192   if(number < 0 && number >= fNumber)
193   {                                               193   {
194     G4Exception("G4ChebyshevApproximation::Get    194     G4Exception("G4ChebyshevApproximation::GetChebyshevCof()", "InvalidCall",
195                 FatalException, "Argument out     195                 FatalException, "Argument out of range !");
196   }                                               196   }
197   return fChebyshevCof[number];                   197   return fChebyshevCof[number];
198 }                                                 198 }
199                                                   199 
200 // -------------------------------------------    200 // --------------------------------------------------------------
201 //                                                201 //
202 // Evaluate the value of fFunction at the poin    202 // Evaluate the value of fFunction at the point x via the Chebyshev coefficients
203 // fChebyshevCof[0,...,fNumber-1]                 203 // fChebyshevCof[0,...,fNumber-1]
204                                                   204 
205 G4double G4ChebyshevApproximation::ChebyshevEv    205 G4double G4ChebyshevApproximation::ChebyshevEvaluation(G4double x) const
206 {                                                 206 {
207   G4double evaluate = 0.0, evaluate2 = 0.0, te    207   G4double evaluate = 0.0, evaluate2 = 0.0, temp = 0.0, xReduced = 0.0,
208            xReduced2 = 0.0;                       208            xReduced2 = 0.0;
209                                                   209 
210   if((x - fMean + fDiff) * (x - fMean - fDiff)    210   if((x - fMean + fDiff) * (x - fMean - fDiff) > 0.0)
211   {                                               211   {
212     G4Exception("G4ChebyshevApproximation::Che    212     G4Exception("G4ChebyshevApproximation::ChebyshevEvaluation()",
213                 "InvalidCall", FatalException,    213                 "InvalidCall", FatalException, "Invalid argument !");
214   }                                               214   }
215   xReduced  = (x - fMean) / fDiff;                215   xReduced  = (x - fMean) / fDiff;
216   xReduced2 = 2.0 * xReduced;                     216   xReduced2 = 2.0 * xReduced;
217   for(G4int i = fNumber - 1; i >= 1; --i)         217   for(G4int i = fNumber - 1; i >= 1; --i)
218   {                                               218   {
219     temp      = evaluate;                         219     temp      = evaluate;
220     evaluate  = xReduced2 * evaluate - evaluat    220     evaluate  = xReduced2 * evaluate - evaluate2 + fChebyshevCof[i];
221     evaluate2 = temp;                             221     evaluate2 = temp;
222   }                                               222   }
223   return xReduced * evaluate - evaluate2 + 0.5    223   return xReduced * evaluate - evaluate2 + 0.5 * fChebyshevCof[0];
224 }                                                 224 }
225                                                   225 
226 // -------------------------------------------    226 // ------------------------------------------------------------------
227 //                                                227 //
228 // Returns the array derCof[0,...,fNumber-2],     228 // Returns the array derCof[0,...,fNumber-2], the Chebyshev coefficients of the
229 // derivative of the function whose coefficien    229 // derivative of the function whose coefficients are fChebyshevCof
230                                                   230 
231 void G4ChebyshevApproximation::DerivativeCheby    231 void G4ChebyshevApproximation::DerivativeChebyshevCof(G4double derCof[]) const
232 {                                                 232 {
233   G4double cof        = 1.0 / fDiff;              233   G4double cof        = 1.0 / fDiff;
234   derCof[fNumber - 1] = 0.0;                      234   derCof[fNumber - 1] = 0.0;
235   derCof[fNumber - 2] = 2 * (fNumber - 1) * fC    235   derCof[fNumber - 2] = 2 * (fNumber - 1) * fChebyshevCof[fNumber - 1];
236   for(G4int i = fNumber - 3; i >= 0; --i)         236   for(G4int i = fNumber - 3; i >= 0; --i)
237   {                                               237   {
238     derCof[i] = derCof[i + 2] + 2 * (i + 1) *     238     derCof[i] = derCof[i + 2] + 2 * (i + 1) * fChebyshevCof[i + 1];
239   }                                               239   }
240   for(G4int j = 0; j < fNumber; ++j)              240   for(G4int j = 0; j < fNumber; ++j)
241   {                                               241   {
242     derCof[j] *= cof;                             242     derCof[j] *= cof;
243   }                                               243   }
244 }                                                 244 }
245                                                   245 
246 // -------------------------------------------    246 // ------------------------------------------------------------------------
247 //                                                247 //
248 // This function produces the array integralCo    248 // This function produces the array integralCof[0,...,fNumber-1] , the Chebyshev
249 // coefficients of the integral of the functio    249 // coefficients of the integral of the function whose coefficients are
250 // fChebyshevCof[]. The constant of integratio    250 // fChebyshevCof[]. The constant of integration is set so that the integral
251 // vanishes at the point (fMean - fDiff), i.e.    251 // vanishes at the point (fMean - fDiff), i.e. at the begining of the interval
252 // of validity (we start the integration from     252 // of validity (we start the integration from this point).
253 //                                                253 //
254                                                   254 
255 void G4ChebyshevApproximation::IntegralChebysh    255 void G4ChebyshevApproximation::IntegralChebyshevCof(
256   G4double integralCof[]) const                   256   G4double integralCof[]) const
257 {                                                 257 {
258   G4double cof = 0.5 * fDiff, sum = 0.0, facto    258   G4double cof = 0.5 * fDiff, sum = 0.0, factor = 1.0;
259   for(G4int i = 1; i < fNumber - 1; ++i)          259   for(G4int i = 1; i < fNumber - 1; ++i)
260   {                                               260   {
261     integralCof[i] = cof * (fChebyshevCof[i -     261     integralCof[i] = cof * (fChebyshevCof[i - 1] - fChebyshevCof[i + 1]) / i;
262     sum += factor * integralCof[i];               262     sum += factor * integralCof[i];
263     factor = -factor;                             263     factor = -factor;
264   }                                               264   }
265   integralCof[fNumber - 1] = cof * fChebyshevC    265   integralCof[fNumber - 1] = cof * fChebyshevCof[fNumber - 2] / (fNumber - 1);
266   sum += factor * integralCof[fNumber - 1];       266   sum += factor * integralCof[fNumber - 1];
267   integralCof[0] = 2.0 * sum;  // set the cons    267   integralCof[0] = 2.0 * sum;  // set the constant of integration
268 }                                                 268 }
269                                                   269