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