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 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // G4ChebyshevApproximation class implementation
 27 //
 28 // Author: V.Grichine, 24.04.97
 29 // --------------------------------------------------------------------
 30 
 31 #include "G4ChebyshevApproximation.hh"
 32 #include "G4PhysicalConstants.hh"
 33 
 34 // Constructor for initialisation of the class data members.
 35 // It creates the array fChebyshevCof[0,...,fNumber-1], fNumber = n ;
 36 // which consists of Chebyshev coefficients describing the function
 37 // pointed by pFunction. The values a and b fix the interval of validity
 38 // of the Chebyshev approximation.
 39 
 40 G4ChebyshevApproximation::G4ChebyshevApproximation(function pFunction, G4int n,
 41                                                    G4double a, G4double b)
 42   : fFunction(pFunction)
 43   , fNumber(n)
 44   , fChebyshevCof(new G4double[fNumber])
 45   , fMean(0.5 * (b + a))
 46   , fDiff(0.5 * (b - a))
 47 {
 48   G4int i = 0, j = 0;
 49   G4double rootSum = 0.0, cofj = 0.0;
 50   auto tempFunction = new G4double[fNumber];
 51   G4double weight        = 2.0 / fNumber;
 52   G4double cof           = 0.5 * weight * pi;  // pi/n
 53 
 54   for(i = 0; i < fNumber; ++i)
 55   {
 56     rootSum         = std::cos(cof * (i + 0.5));
 57     tempFunction[i] = fFunction(rootSum * fDiff + fMean);
 58   }
 59   for(j = 0; j < fNumber; ++j)
 60   {
 61     cofj    = cof * j;
 62     rootSum = 0.0;
 63 
 64     for(i = 0; i < fNumber; ++i)
 65     {
 66       rootSum += tempFunction[i] * std::cos(cofj * (i + 0.5));
 67     }
 68     fChebyshevCof[j] = weight * rootSum;
 69   }
 70   delete[] tempFunction;
 71 }
 72 
 73 // --------------------------------------------------------------------
 74 //
 75 // Constructor for creation of Chebyshev coefficients for mx-derivative
 76 // from pFunction. The value of mx ! MUST BE ! < nx , because the result
 77 // array of fChebyshevCof will be of (nx-mx) size.  The values a and b
 78 // fix the interval of validity of the Chebyshev approximation.
 79 
 80 G4ChebyshevApproximation::G4ChebyshevApproximation(function pFunction, G4int nx,
 81                                                    G4int mx, G4double a,
 82                                                    G4double b)
 83   : fFunction(pFunction)
 84   , fNumber(nx)
 85   , fChebyshevCof(new G4double[fNumber])
 86   , fMean(0.5 * (b + a))
 87   , fDiff(0.5 * (b - a))
 88 {
 89   if(nx <= mx)
 90   {
 91     G4Exception("G4ChebyshevApproximation::G4ChebyshevApproximation()",
 92                 "InvalidCall", FatalException, "Invalid arguments !");
 93   }
 94   G4int i = 0, j = 0;
 95   G4double rootSum = 0.0, cofj = 0.0;
 96   auto tempFunction = new G4double[fNumber];
 97   G4double weight        = 2.0 / fNumber;
 98   G4double cof           = 0.5 * weight * pi;  // pi/nx
 99 
100   for(i = 0; i < fNumber; ++i)
101   {
102     rootSum         = std::cos(cof * (i + 0.5));
103     tempFunction[i] = fFunction(rootSum * fDiff + fMean);
104   }
105   for(j = 0; j < fNumber; ++j)
106   {
107     cofj    = cof * j;
108     rootSum = 0.0;
109 
110     for(i = 0; i < fNumber; ++i)
111     {
112       rootSum += tempFunction[i] * std::cos(cofj * (i + 0.5));
113     }
114     fChebyshevCof[j] = weight * rootSum;  // corresponds to pFunction
115   }
116   // Chebyshev coefficients for (mx)-derivative of pFunction
117 
118   for(i = 1; i <= mx; ++i)
119   {
120     DerivativeChebyshevCof(tempFunction);
121     fNumber--;
122     for(j = 0; j < fNumber; ++j)
123     {
124       fChebyshevCof[j] = tempFunction[j];  // corresponds to (i)-derivative
125     }
126   }
127   delete[] tempFunction;  // delete of dynamically allocated tempFunction
128 }
129 
130 // ------------------------------------------------------
131 //
132 // Constructor for creation of Chebyshev coefficients for integral
133 // from pFunction.
134 
135 G4ChebyshevApproximation::G4ChebyshevApproximation(function pFunction,
136                                                    G4double a, G4double b,
137                                                    G4int n)
138   : fFunction(pFunction)
139   , fNumber(n)
140   , fChebyshevCof(new G4double[fNumber])
141   , fMean(0.5 * (b + a))
142   , fDiff(0.5 * (b - a))
143 {
144   G4int i = 0, j = 0;
145   G4double rootSum = 0.0, cofj = 0.0;
146   auto tempFunction = new G4double[fNumber];
147   G4double weight        = 2.0 / fNumber;
148   G4double cof           = 0.5 * weight * pi;  // pi/n
149 
150   for(i = 0; i < fNumber; ++i)
151   {
152     rootSum         = std::cos(cof * (i + 0.5));
153     tempFunction[i] = fFunction(rootSum * fDiff + fMean);
154   }
155   for(j = 0; j < fNumber; ++j)
156   {
157     cofj    = cof * j;
158     rootSum = 0.0;
159 
160     for(i = 0; i < fNumber; ++i)
161     {
162       rootSum += tempFunction[i] * std::cos(cofj * (i + 0.5));
163     }
164     fChebyshevCof[j] = weight * rootSum;  // corresponds to pFunction
165   }
166   // Chebyshev coefficients for integral of pFunction
167 
168   IntegralChebyshevCof(tempFunction);
169   for(j = 0; j < fNumber; ++j)
170   {
171     fChebyshevCof[j] = tempFunction[j];  // corresponds to integral
172   }
173   delete[] tempFunction;  // delete of dynamically allocated tempFunction
174 }
175 
176 // ---------------------------------------------------------------
177 //
178 // Destructor deletes the array of Chebyshev coefficients
179 
180 G4ChebyshevApproximation::~G4ChebyshevApproximation()
181 {
182   delete[] fChebyshevCof;
183 }
184 
185 // ---------------------------------------------------------------
186 //
187 // Access function for Chebyshev coefficients
188 //
189 
190 G4double G4ChebyshevApproximation::GetChebyshevCof(G4int number) const
191 {
192   if(number < 0 && number >= fNumber)
193   {
194     G4Exception("G4ChebyshevApproximation::GetChebyshevCof()", "InvalidCall",
195                 FatalException, "Argument out of range !");
196   }
197   return fChebyshevCof[number];
198 }
199 
200 // --------------------------------------------------------------
201 //
202 // Evaluate the value of fFunction at the point x via the Chebyshev coefficients
203 // fChebyshevCof[0,...,fNumber-1]
204 
205 G4double G4ChebyshevApproximation::ChebyshevEvaluation(G4double x) const
206 {
207   G4double evaluate = 0.0, evaluate2 = 0.0, temp = 0.0, xReduced = 0.0,
208            xReduced2 = 0.0;
209 
210   if((x - fMean + fDiff) * (x - fMean - fDiff) > 0.0)
211   {
212     G4Exception("G4ChebyshevApproximation::ChebyshevEvaluation()",
213                 "InvalidCall", FatalException, "Invalid argument !");
214   }
215   xReduced  = (x - fMean) / fDiff;
216   xReduced2 = 2.0 * xReduced;
217   for(G4int i = fNumber - 1; i >= 1; --i)
218   {
219     temp      = evaluate;
220     evaluate  = xReduced2 * evaluate - evaluate2 + fChebyshevCof[i];
221     evaluate2 = temp;
222   }
223   return xReduced * evaluate - evaluate2 + 0.5 * fChebyshevCof[0];
224 }
225 
226 // ------------------------------------------------------------------
227 //
228 // Returns the array derCof[0,...,fNumber-2], the Chebyshev coefficients of the
229 // derivative of the function whose coefficients are fChebyshevCof
230 
231 void G4ChebyshevApproximation::DerivativeChebyshevCof(G4double derCof[]) const
232 {
233   G4double cof        = 1.0 / fDiff;
234   derCof[fNumber - 1] = 0.0;
235   derCof[fNumber - 2] = 2 * (fNumber - 1) * fChebyshevCof[fNumber - 1];
236   for(G4int i = fNumber - 3; i >= 0; --i)
237   {
238     derCof[i] = derCof[i + 2] + 2 * (i + 1) * fChebyshevCof[i + 1];
239   }
240   for(G4int j = 0; j < fNumber; ++j)
241   {
242     derCof[j] *= cof;
243   }
244 }
245 
246 // ------------------------------------------------------------------------
247 //
248 // This function produces the array integralCof[0,...,fNumber-1] , the Chebyshev
249 // coefficients of the integral of the function whose coefficients are
250 // fChebyshevCof[]. The constant of integration is set so that the integral
251 // vanishes at the point (fMean - fDiff), i.e. at the begining of the interval
252 // of validity (we start the integration from this point).
253 //
254 
255 void G4ChebyshevApproximation::IntegralChebyshevCof(
256   G4double integralCof[]) const
257 {
258   G4double cof = 0.5 * fDiff, sum = 0.0, factor = 1.0;
259   for(G4int i = 1; i < fNumber - 1; ++i)
260   {
261     integralCof[i] = cof * (fChebyshevCof[i - 1] - fChebyshevCof[i + 1]) / i;
262     sum += factor * integralCof[i];
263     factor = -factor;
264   }
265   integralCof[fNumber - 1] = cof * fChebyshevCof[fNumber - 2] / (fNumber - 1);
266   sum += factor * integralCof[fNumber - 1];
267   integralCof[0] = 2.0 * sum;  // set the constant of integration
268 }
269