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 #include "G4RKG3_Stepper.hh"
00032 #include "G4LineSection.hh"
00033 #include "G4Mag_EqRhs.hh"
00034
00035 G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs *EqRhs)
00036 : G4MagIntegratorStepper(EqRhs,6), hStep(0.)
00037 {
00038 }
00039
00040 G4RKG3_Stepper::~G4RKG3_Stepper()
00041 {
00042 }
00043
00044 void G4RKG3_Stepper::Stepper( const G4double yInput[8],
00045 const G4double dydx[6],
00046 G4double Step,
00047 G4double yOut[8],
00048 G4double yErr[])
00049 {
00050 G4double B[3];
00051 G4int nvar = 6 ;
00052 G4int i;
00053 G4double by15 = 1. / 15. ;
00054
00055 G4double yTemp[8], dydxTemp[6], yIn[8] ;
00056
00057 for(i=0;i<nvar;i++) yIn[i]=yInput[i];
00058 yIn[6] = yInput[6];
00059 yIn[7] = yInput[7];
00060 G4double h = Step * 0.5;
00061 hStep=Step;
00062
00063
00064 StepNoErr(yIn, dydx,h, yTemp,B) ;
00065
00066
00067 for(i=0;i<3;i++)BfldIn[i]=B[i];
00068
00069
00070
00071 GetEquationOfMotion()->EvaluateRhsGivenB(yTemp,B,dydxTemp) ;
00072 StepNoErr(yTemp,dydxTemp,h,yOut,B);
00073
00074
00075
00076 fyMidPoint = G4ThreeVector( yTemp[0], yTemp[1], yTemp[2]);
00077
00078
00079
00080 h *= 2 ;
00081 StepNoErr(yIn,dydx,h,yTemp,B);
00082 for(i=0;i<nvar;i++)
00083 {
00084 yErr[i] = yOut[i] - yTemp[i] ;
00085 yOut[i] += yErr[i]*by15 ;
00086 }
00087
00088
00089
00090 fyInitial = G4ThreeVector( yIn[0], yIn[1], yIn[2]);
00091 fpInitial = G4ThreeVector( yIn[3], yIn[4], yIn[5]);
00092 fyFinal = G4ThreeVector( yOut[0], yOut[1], yOut[2]);
00093
00094
00095 }
00096
00097
00098
00099
00100
00101
00102
00103
00104 void G4RKG3_Stepper::StepWithEst( const G4double*,
00105 const G4double*,
00106 G4double,
00107 G4double*,
00108 G4double&,
00109 G4double&,
00110 const G4double*,
00111 G4double* )
00112
00113 {
00114 G4Exception("G4RKG3_Stepper::StepWithEst()", "GeomField0001",
00115 FatalException, "Method no longer used.");
00116 }
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126 void G4RKG3_Stepper::StepNoErr(const G4double tIn[8],
00127 const G4double dydx[6],
00128 G4double Step,
00129 G4double tOut[8],
00130 G4double B[3] )
00131
00132 {
00133
00134
00135 G4double K1[7],K2[7],K3[7],K4[7] ;
00136 G4double tTemp[8], yderiv[6] ;
00137
00138
00139
00140 G4double mom,inverse_mom;
00141 G4int i ;
00142 const G4double c1=0.5,c2=0.125,c3=1./6.;
00143
00144
00145
00146
00147 mom=std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4]+tIn[5]*tIn[5]);
00148 inverse_mom=1./mom;
00149 for(i=0;i<3;i++)
00150 {
00151 K1[i] = Step * dydx[i+3]*inverse_mom;
00152 tTemp[i] = tIn[i] + Step*(c1*tIn[i+3]*inverse_mom + c2*K1[i]) ;
00153 tTemp[i+3] = tIn[i+3] + c1*K1[i]*mom ;
00154
00155 }
00156
00157 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
00158
00159
00160 for(i=0;i<3;i++)
00161 {
00162 K2[i] = Step * yderiv[i+3]*inverse_mom;
00163 tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ;
00164 }
00165
00166
00167 GetEquationOfMotion()->EvaluateRhsGivenB(tTemp,B,yderiv) ;
00168
00169 for(i=0;i<3;i++)
00170 {
00171 K3[i] = Step * yderiv[i+3]*inverse_mom;
00172 tTemp[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom + c1*K3[i]) ;
00173 tTemp[i+3] = tIn[i+3] + K3[i]*mom ;
00174 }
00175
00176
00177
00178 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
00179
00180
00181 for(i=0;i<3;i++)
00182 {
00183 K4[i] = Step * yderiv[i+3]*inverse_mom;
00184 tOut[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom+ (K1[i] + K2[i] + K3[i])*c3) ;
00185 tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2[i] + 2*K3[i] +K4[i])*c3 ;
00186 }
00187 tOut[6] = tIn[6];
00188 tOut[7] = tIn[7];
00189
00190
00191
00192 }
00193
00194
00195
00196
00197 G4double G4RKG3_Stepper::DistChord() const
00198 {
00199
00200
00201 G4double distChord,distLine;
00202
00203 if (fyInitial != fyFinal) {
00204 distLine= G4LineSection::Distline(fyMidPoint,fyInitial,fyFinal );
00205
00206 distChord = distLine;
00207 }else{
00208 distChord = (fyMidPoint-fyInitial).mag();
00209 }
00210
00211
00212 return distChord;
00213
00214 }
00215