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
00031
00032 #include "globals.hh"
00033 #include "G4SimpleIntegration.hh"
00034
00035
00036 G4int G4SimpleIntegration::fMaxDepth = 100 ;
00037
00038
00039 G4SimpleIntegration::G4SimpleIntegration( function pFunction )
00040 : fFunction(pFunction),
00041 fTolerance(.0001)
00042 {
00043 }
00044
00045 G4SimpleIntegration::G4SimpleIntegration( function pFunction,
00046 G4double pTolerance)
00047 : fFunction(pFunction),
00048 fTolerance(pTolerance)
00049 {
00050 }
00051
00052
00053 G4SimpleIntegration::~G4SimpleIntegration()
00054 {
00055 }
00056
00057
00058
00059 G4double
00060 G4SimpleIntegration::Trapezoidal(G4double xInitial,
00061 G4double xFinal,
00062 G4int iterationNumber )
00063 {
00064 G4double Step = (xFinal - xInitial)/iterationNumber ;
00065 G4double mean = (fFunction(xInitial) + fFunction(xFinal))*0.5 ;
00066 G4double x = xInitial ;
00067 for(G4int i=1;i<iterationNumber;i++)
00068 {
00069 x += Step ;
00070 mean += fFunction(x) ;
00071 }
00072 return mean*Step ;
00073 }
00074
00075 G4double
00076 G4SimpleIntegration::MidPoint(G4double xInitial,
00077 G4double xFinal,
00078 G4int iterationNumber )
00079 {
00080 G4double Step = (xFinal - xInitial)/iterationNumber ;
00081 G4double x = xInitial + 0.5*Step;
00082 G4double mean = fFunction(x) ;
00083 for(G4int i=1;i<iterationNumber;i++)
00084 {
00085 x += Step ;
00086 mean += fFunction(x) ;
00087 }
00088 return mean*Step ;
00089 }
00090
00091 G4double
00092 G4SimpleIntegration::Gauss(G4double xInitial,
00093 G4double xFinal,
00094 G4int iterationNumber )
00095 {
00096 G4double x=0.;
00097 static G4double root = 1.0/std::sqrt(3.0) ;
00098 G4double Step = (xFinal - xInitial)/(2.0*iterationNumber) ;
00099 G4double delta = Step*root ;
00100 G4double mean = 0.0 ;
00101 for(G4int i=0;i<iterationNumber;i++)
00102 {
00103 x = (2*i + 1)*Step ;
00104 mean += (fFunction(x+delta) + fFunction(x-delta)) ;
00105 }
00106 return mean*Step ;
00107 }
00108
00109 G4double
00110 G4SimpleIntegration::Simpson(G4double xInitial,
00111 G4double xFinal,
00112 G4int iterationNumber )
00113 {
00114 G4double Step = (xFinal - xInitial)/iterationNumber ;
00115 G4double x = xInitial ;
00116 G4double xPlus = xInitial + 0.5*Step ;
00117 G4double mean = (fFunction(xInitial) + fFunction(xFinal))*0.5 ;
00118 G4double sum = fFunction(xPlus) ;
00119 for(G4int i=1;i<iterationNumber;i++)
00120 {
00121 x += Step ;
00122 xPlus += Step ;
00123 mean += fFunction(x) ;
00124 sum += fFunction(xPlus) ;
00125 }
00126 mean += 2.0*sum ;
00127 return mean*Step/3.0 ;
00128 }
00129
00130
00131
00132
00133
00134 G4double
00135 G4SimpleIntegration::AdaptGaussIntegration( G4double xInitial,
00136 G4double xFinal )
00137 {
00138 G4int depth = 0 ;
00139 G4double sum = 0.0 ;
00140 AdaptGauss(xInitial,xFinal,sum,depth) ;
00141 return sum ;
00142 }
00143
00144
00145 G4double
00146 G4SimpleIntegration::Gauss( G4double xInitial,
00147 G4double xFinal )
00148 {
00149 static G4double root = 1.0/std::sqrt(3.0) ;
00150
00151 G4double xMean = (xInitial + xFinal)/2.0 ;
00152 G4double Step = (xFinal - xInitial)/2.0 ;
00153 G4double delta = Step*root ;
00154 G4double sum = (fFunction(xMean + delta) + fFunction(xMean - delta)) ;
00155
00156 return sum*Step ;
00157 }
00158
00159
00160 void
00161 G4SimpleIntegration::AdaptGauss( G4double xInitial,
00162 G4double xFinal,
00163 G4double& sum,
00164 G4int& depth )
00165 {
00166 if(depth >fMaxDepth)
00167 {
00168 G4Exception("G4SimpleIntegration::AdaptGauss()", "Error",
00169 FatalException, "Function varies too rapidly !") ;
00170 }
00171 G4double xMean = (xInitial + xFinal)/2.0 ;
00172 G4double leftHalf = Gauss(xInitial,xMean) ;
00173 G4double rightHalf = Gauss(xMean,xFinal) ;
00174 G4double full = Gauss(xInitial,xFinal) ;
00175 if(std::fabs(leftHalf+rightHalf-full) < fTolerance)
00176 {
00177 sum += full ;
00178 }
00179 else
00180 {
00181 depth++ ;
00182 AdaptGauss(xInitial,xMean,sum,depth) ;
00183 AdaptGauss(xMean,xFinal,sum,depth) ;
00184 }
00185 }