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
00033
00034
00035
00036
00037
00038
00039
00040
00041 #include "G4HelixMixedStepper.hh"
00042 #include "G4PhysicalConstants.hh"
00043 #include "G4ClassicalRK4.hh"
00044 #include "G4CashKarpRKF45.hh"
00045 #include "G4SimpleRunge.hh"
00046 #include "G4HelixImplicitEuler.hh"
00047 #include "G4HelixExplicitEuler.hh"
00048 #include "G4HelixSimpleRunge.hh"
00049 #include "G4ExactHelixStepper.hh"
00050 #include "G4ExplicitEuler.hh"
00051 #include "G4ImplicitEuler.hh"
00052 #include "G4SimpleHeum.hh"
00053 #include "G4RKG3_Stepper.hh"
00054
00055 #include "G4ThreeVector.hh"
00056 #include "G4LineSection.hh"
00057 G4HelixMixedStepper::G4HelixMixedStepper(G4Mag_EqRhs *EqRhs,G4int fStepperNumber)
00058 : G4MagHelicalStepper(EqRhs)
00059
00060 {
00061 SetVerbose(1); fNumCallsRK4=0; fNumCallsHelix=0;
00062 if(!fStepperNumber) fStepperNumber=4;
00063 fRK4Stepper = SetupStepper(EqRhs, fStepperNumber);
00064 }
00065
00066
00067 G4HelixMixedStepper::~G4HelixMixedStepper() {
00068
00069 delete(fRK4Stepper);
00070 if (fVerbose>0){ PrintCalls();};
00071 }
00072 void G4HelixMixedStepper::Stepper( const G4double yInput[7],
00073 const G4double dydx[7],
00074 G4double Step,
00075 G4double yOut[7],
00076 G4double yErr[])
00077
00078 {
00079
00080
00081
00082 G4ThreeVector Bfld;
00083 MagFieldEvaluate(yInput, Bfld);
00084
00085 G4double Bmag = Bfld.mag();
00086 const G4double *pIn = yInput+3;
00087 G4ThreeVector initVelocity= G4ThreeVector( pIn[0], pIn[1], pIn[2]);
00088 G4double velocityVal = initVelocity.mag();
00089 G4double R_1;
00090 G4double Ang_curve;
00091
00092 R_1=std::abs(GetInverseCurve(velocityVal,Bmag));
00093 Ang_curve=R_1*Step;
00094 SetAngCurve(Ang_curve);
00095 SetCurve(std::abs(1/R_1));
00096
00097
00098 if(Ang_curve<0.33*pi){
00099 fNumCallsRK4++;
00100 fRK4Stepper->Stepper(yInput,dydx,Step,yOut,yErr);
00101
00102
00103 }
00104 else{
00105 fNumCallsHelix++;
00106 const G4int nvar = 6 ;
00107 G4int i;
00108 G4double yTemp[7], yIn[7] ;
00109 G4double yTemp2[7];
00110 G4ThreeVector Bfld_midpoint;
00111
00112 for(i=0;i<nvar;i++) yIn[i]=yInput[i];
00113
00114 G4double h = Step * 0.5;
00115
00116 AdvanceHelix(yIn, Bfld, h, yTemp,yTemp2);
00117 MagFieldEvaluate(yTemp, Bfld_midpoint) ;
00118 AdvanceHelix(yTemp, Bfld_midpoint, h, yOut);
00119
00120 for(i=0;i<nvar;i++) {
00121 yErr[i] = yOut[i] - yTemp2[i] ;
00122
00123 }
00124 }
00125
00126
00127
00128
00129 }
00130
00131 void
00132 G4HelixMixedStepper::DumbStepper( const G4double yIn[],
00133 G4ThreeVector Bfld,
00134 G4double h,
00135 G4double yOut[])
00136 {
00137
00138
00139 AdvanceHelix(yIn, Bfld, h, yOut);
00140
00141
00142
00143 }
00144
00145 G4double G4HelixMixedStepper::DistChord() const
00146 {
00147
00148
00149
00150
00151 G4double distChord;
00152 G4double Ang_curve=GetAngCurve();
00153
00154
00155 if(Ang_curve<=pi){
00156 distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve));
00157 }
00158 else
00159 if(Ang_curve<twopi){
00160 distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve)));
00161 }
00162 else{
00163 distChord=2.*GetRadHelix();
00164 }
00165
00166
00167
00168 return distChord;
00169
00170 }
00171
00172 void G4HelixMixedStepper::PrintCalls()
00173 {
00174 G4cout<<"In HelixMixedStepper::Number of calls to smallStepStepper = "<<fNumCallsRK4
00175 <<" and Number of calls to Helix = "<<fNumCallsHelix<<G4endl;
00176 }
00177
00178
00179
00180 G4MagIntegratorStepper* G4HelixMixedStepper:: SetupStepper(G4Mag_EqRhs* pE, G4int StepperNumber)
00181 {
00182 G4MagIntegratorStepper* pStepper;
00183 if (fVerbose>0)G4cout<<"In G4HelixMixedStepper Stepper for small steps is ";
00184 switch ( StepperNumber )
00185 {
00186 case 0: pStepper = new G4ExplicitEuler( pE ); if (fVerbose>0)G4cout<<"G4ExplicitEuler"<<G4endl; break;
00187 case 1: pStepper = new G4ImplicitEuler( pE ); if (fVerbose>0)G4cout<<"G4ImplicitEuler"<<G4endl; break;
00188 case 2: pStepper = new G4SimpleRunge( pE ); if (fVerbose>0)G4cout<<"G4SimpleRunge"<<G4endl; break;
00189 case 3: pStepper = new G4SimpleHeum( pE ); if (fVerbose>0)G4cout<<"G4SimpleHeum"<<G4endl;break;
00190 case 4: pStepper = new G4ClassicalRK4( pE ); if (fVerbose>0)G4cout<<"G4ClassicalRK4"<<G4endl; break;
00191 case 5: pStepper = new G4HelixExplicitEuler( pE ); if (fVerbose>0)G4cout<<"G4HelixExplicitEuler"<<G4endl; break;
00192 case 6: pStepper = new G4HelixImplicitEuler( pE ); if (fVerbose>0)G4cout<<"G4HelixImplicitEuler"<<G4endl; break;
00193 case 7: pStepper = new G4HelixSimpleRunge( pE ); if (fVerbose>0)G4cout<<"G4HelixSimpleRunge"<<G4endl; break;
00194 case 8: pStepper = new G4CashKarpRKF45( pE ); if (fVerbose>0)G4cout<<"G4CashKarpRKF45"<<G4endl; break;
00195 case 9: pStepper = new G4ExactHelixStepper( pE ); if (fVerbose>0)G4cout<<"G4ExactHelixStepper"<<G4endl; break;
00196 case 10: pStepper = new G4RKG3_Stepper( pE ); if (fVerbose>0)G4cout<<"G4RKG3_Stepper"<<G4endl; break;
00197
00198 default: pStepper = new G4ClassicalRK4( pE );G4cout<<"Default G4ClassicalRK4"<<G4endl; break;
00199
00200 }
00201 return pStepper;
00202 }