00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include "G4ChebyshevApproximation.hh"
00031 #include "G4PhysicalConstants.hh"
00032
00033
00034
00035
00036
00037
00038
00039 G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
00040 G4int n,
00041 G4double a,
00042 G4double b )
00043 : fFunction(pFunction), fNumber(n),
00044 fChebyshevCof(new G4double[fNumber]),
00045 fMean(0.5*(b+a)), fDiff(0.5*(b-a))
00046 {
00047 G4int i=0, j=0 ;
00048 G4double rootSum=0.0, cofj=0.0 ;
00049 G4double* tempFunction = new G4double[fNumber] ;
00050 G4double weight = 2.0/fNumber ;
00051 G4double cof = 0.5*weight*pi ;
00052
00053 for (i=0;i<fNumber;i++)
00054 {
00055 rootSum = std::cos(cof*(i+0.5)) ;
00056 tempFunction[i]= fFunction(rootSum*fDiff+fMean) ;
00057 }
00058 for (j=0;j<fNumber;j++)
00059 {
00060 cofj = cof*j ;
00061 rootSum = 0.0 ;
00062
00063 for (i=0;i<fNumber;i++)
00064 {
00065 rootSum += tempFunction[i]*std::cos(cofj*(i+0.5)) ;
00066 }
00067 fChebyshevCof[j] = weight*rootSum ;
00068 }
00069 delete[] tempFunction ;
00070 }
00071
00072
00073
00074
00075
00076
00077
00078
00079 G4ChebyshevApproximation::
00080 G4ChebyshevApproximation( function pFunction,
00081 G4int nx, G4int mx,
00082 G4double a, G4double b )
00083 : fFunction(pFunction), fNumber(nx),
00084 fChebyshevCof(new G4double[fNumber]),
00085 fMean(0.5*(b+a)), fDiff(0.5*(b-a))
00086 {
00087 if(nx <= mx)
00088 {
00089 G4Exception("G4ChebyshevApproximation::G4ChebyshevApproximation()",
00090 "InvalidCall", FatalException, "Invalid arguments !") ;
00091 }
00092 G4int i=0, j=0 ;
00093 G4double rootSum = 0.0, cofj=0.0;
00094 G4double* tempFunction = new G4double[fNumber] ;
00095 G4double weight = 2.0/fNumber ;
00096 G4double cof = 0.5*weight*pi ;
00097
00098 for (i=0;i<fNumber;i++)
00099 {
00100 rootSum = std::cos(cof*(i+0.5)) ;
00101 tempFunction[i] = fFunction(rootSum*fDiff+fMean) ;
00102 }
00103 for (j=0;j<fNumber;j++)
00104 {
00105 cofj = cof*j ;
00106 rootSum = 0.0 ;
00107
00108 for (i=0;i<fNumber;i++)
00109 {
00110 rootSum += tempFunction[i]*std::cos(cofj*(i+0.5)) ;
00111 }
00112 fChebyshevCof[j] = weight*rootSum ;
00113 }
00114
00115
00116 for(i=1;i<=mx;i++)
00117 {
00118 DerivativeChebyshevCof(tempFunction) ;
00119 fNumber-- ;
00120 for(j=0;j<fNumber;j++)
00121 {
00122 fChebyshevCof[j] = tempFunction[j] ;
00123 }
00124 }
00125 delete[] tempFunction ;
00126 }
00127
00128
00129
00130
00131
00132
00133 G4ChebyshevApproximation::G4ChebyshevApproximation( function pFunction,
00134 G4double a,
00135 G4double b,
00136 G4int n )
00137 : fFunction(pFunction), fNumber(n),
00138 fChebyshevCof(new G4double[fNumber]),
00139 fMean(0.5*(b+a)), fDiff(0.5*(b-a))
00140 {
00141 G4int i=0, j=0;
00142 G4double rootSum=0.0, cofj=0.0;
00143 G4double* tempFunction = new G4double[fNumber] ;
00144 G4double weight = 2.0/fNumber;
00145 G4double cof = 0.5*weight*pi ;
00146
00147 for (i=0;i<fNumber;i++)
00148 {
00149 rootSum = std::cos(cof*(i+0.5)) ;
00150 tempFunction[i]= fFunction(rootSum*fDiff+fMean) ;
00151 }
00152 for (j=0;j<fNumber;j++)
00153 {
00154 cofj = cof*j ;
00155 rootSum = 0.0 ;
00156
00157 for (i=0;i<fNumber;i++)
00158 {
00159 rootSum += tempFunction[i]*std::cos(cofj*(i+0.5)) ;
00160 }
00161 fChebyshevCof[j] = weight*rootSum ;
00162 }
00163
00164
00165 IntegralChebyshevCof(tempFunction) ;
00166 for(j=0;j<fNumber;j++)
00167 {
00168 fChebyshevCof[j] = tempFunction[j] ;
00169 }
00170 delete[] tempFunction ;
00171 }
00172
00173
00174
00175
00176
00177
00178
00179 G4ChebyshevApproximation::~G4ChebyshevApproximation()
00180 {
00181 delete[] fChebyshevCof ;
00182 }
00183
00184
00185
00186
00187
00188
00189
00190 G4double
00191 G4ChebyshevApproximation::GetChebyshevCof(G4int number) const
00192 {
00193 if(number < 0 && number >= fNumber)
00194 {
00195 G4Exception("G4ChebyshevApproximation::GetChebyshevCof()",
00196 "InvalidCall", FatalException, "Argument out of range !") ;
00197 }
00198 return fChebyshevCof[number] ;
00199 }
00200
00201
00202
00203
00204
00205
00206 G4double
00207 G4ChebyshevApproximation::ChebyshevEvaluation(G4double x) const
00208 {
00209 G4double evaluate = 0.0, evaluate2 = 0.0, temp = 0.0,
00210 xReduced = 0.0, xReduced2 = 0.0 ;
00211
00212 if ((x-fMean+fDiff)*(x-fMean-fDiff) > 0.0)
00213 {
00214 G4Exception("G4ChebyshevApproximation::ChebyshevEvaluation()",
00215 "InvalidCall", FatalException, "Invalid argument !") ;
00216 }
00217 xReduced = (x-fMean)/fDiff ;
00218 xReduced2 = 2.0*xReduced ;
00219 for (G4int i=fNumber-1;i>=1;i--)
00220 {
00221 temp = evaluate ;
00222 evaluate = xReduced2*evaluate - evaluate2 + fChebyshevCof[i] ;
00223 evaluate2 = temp ;
00224 }
00225 return xReduced*evaluate - evaluate2 + 0.5*fChebyshevCof[0] ;
00226 }
00227
00228
00229
00230
00231
00232
00233 void
00234 G4ChebyshevApproximation::DerivativeChebyshevCof(G4double derCof[]) const
00235 {
00236 G4double cof = 1.0/fDiff ;
00237 derCof[fNumber-1] = 0.0 ;
00238 derCof[fNumber-2] = 2*(fNumber-1)*fChebyshevCof[fNumber-1] ;
00239 for(G4int i=fNumber-3;i>=0;i--)
00240 {
00241 derCof[i] = derCof[i+2] + 2*(i+1)*fChebyshevCof[i+1] ;
00242 }
00243 for(G4int j=0;j<fNumber;j++)
00244 {
00245 derCof[j] *= cof ;
00246 }
00247 }
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 void
00259 G4ChebyshevApproximation::IntegralChebyshevCof(G4double integralCof[]) const
00260 {
00261 G4double cof = 0.5*fDiff, sum = 0.0, factor = 1.0 ;
00262 for(G4int i=1;i<fNumber-1;i++)
00263 {
00264 integralCof[i] = cof*(fChebyshevCof[i-1] - fChebyshevCof[i+1])/i ;
00265 sum += factor*integralCof[i] ;
00266 factor = -factor ;
00267 }
00268 integralCof[fNumber-1] = cof*fChebyshevCof[fNumber-2]/(fNumber-1) ;
00269 sum += factor*integralCof[fNumber-1] ;
00270 integralCof[0] = 2.0*sum ;
00271 }