Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4ConstRK4 Class Reference

#include <G4ConstRK4.hh>

+ Inheritance diagram for G4ConstRK4:

Public Member Functions

 G4ConstRK4 (G4Mag_EqRhs *EquationMotion, G4int numberOfStateVariables=8)
 
 ~G4ConstRK4 ()
 
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
 
void DumbStepper (const G4double yIn[], const G4double dydx[], G4double h, G4double yOut[])
 
G4double DistChord () const
 
void RightHandSideConst (const G4double y[], G4double dydx[]) const
 
void GetConstField (const G4double y[], G4double Field[])
 
G4int IntegratorOrder () const
 
- Public Member Functions inherited from G4MagErrorStepper
 G4MagErrorStepper (G4EquationOfMotion *EqRhs, G4int numberOfVariables, G4int numStateVariables=12)
 
virtual ~G4MagErrorStepper ()
 
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
 
virtual void DumbStepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[])=0
 
G4double DistChord () const
 
- Public Member Functions inherited from G4MagIntegratorStepper
 G4MagIntegratorStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12)
 
virtual ~G4MagIntegratorStepper ()
 
virtual void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])=0
 
virtual G4double DistChord () const =0
 
virtual void ComputeRightHandSide (const G4double y[], G4double dydx[])
 
void NormaliseTangentVector (G4double vec[6])
 
void NormalisePolarizationVector (G4double vec[12])
 
void RightHandSide (const double y[], double dydx[])
 
G4int GetNumberOfVariables () const
 
G4int GetNumberOfStateVariables () const
 
virtual G4int IntegratorOrder () const =0
 
G4EquationOfMotionGetEquationOfMotion ()
 
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
 

Detailed Description

Definition at line 51 of file G4ConstRK4.hh.

Constructor & Destructor Documentation

◆ G4ConstRK4()

G4ConstRK4::G4ConstRK4 ( G4Mag_EqRhs EquationMotion,
G4int  numberOfStateVariables = 8 
)

Definition at line 42 of file G4ConstRK4.cc.

43 : G4MagErrorStepper(EqRhs, 6, numStateVariables)
44{
45 // const G4int numberOfVariables= 6;
46 if( numStateVariables < 8 )
47 {
48 std::ostringstream message;
49 message << "The number of State variables at least 8 " << G4endl
50 << "Instead it is - numStateVariables= " << numStateVariables;
51 G4Exception("G4ConstRK4::G4ConstRK4()", "GeomField0002",
52 FatalException, message, "Use another Stepper!");
53 }
54
55 fEq = EqRhs;
56 yMiddle = new G4double[8];
57 dydxMid = new G4double[8];
58 yInitial = new G4double[8];
59 yOneStep = new G4double[8];
60
61 dydxm = new G4double[8];
62 dydxt = new G4double[8];
63 yt = new G4double[8];
64 Field[0]=0.; Field[1]=0.; Field[2]=0.;
65}
@ FatalException
double G4double
Definition: G4Types.hh:64
#define G4endl
Definition: G4ios.hh:52
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

◆ ~G4ConstRK4()

G4ConstRK4::~G4ConstRK4 ( )

Definition at line 71 of file G4ConstRK4.cc.

72{
73 delete [] yMiddle;
74 delete [] dydxMid;
75 delete [] yInitial;
76 delete [] yOneStep;
77 delete [] dydxm;
78 delete [] dydxt;
79 delete [] yt;
80}

Member Function Documentation

◆ DistChord()

G4double G4ConstRK4::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 213 of file G4ConstRK4.cc.

214{
215 G4double distLine, distChord;
216
217 if (fInitialPoint != fFinalPoint)
218 {
219 distLine= G4LineSection::Distline( fMidPoint, fInitialPoint, fFinalPoint );
220 // This is a class method that gives distance of Mid
221 // from the Chord between the Initial and Final points
222 distChord = distLine;
223 }
224 else
225 {
226 distChord = (fMidPoint-fInitialPoint).mag();
227 }
228 return distChord;
229}
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)

◆ DumbStepper()

void G4ConstRK4::DumbStepper ( const G4double  yIn[],
const G4double  dydx[],
G4double  h,
G4double  yOut[] 
)
virtual

Implements G4MagErrorStepper.

Definition at line 92 of file G4ConstRK4.cc.

96{
97 G4double hh = h*0.5 , h6 = h/6.0 ;
98
99 // 1st Step K1=h*dydx
100 yt[5] = yIn[5] + hh*dydx[5] ;
101 yt[4] = yIn[4] + hh*dydx[4] ;
102 yt[3] = yIn[3] + hh*dydx[3] ;
103 yt[2] = yIn[2] + hh*dydx[2] ;
104 yt[1] = yIn[1] + hh*dydx[1] ;
105 yt[0] = yIn[0] + hh*dydx[0] ;
106 RightHandSideConst(yt,dydxt) ;
107
108 // 2nd Step K2=h*dydxt
109 yt[5] = yIn[5] + hh*dydxt[5] ;
110 yt[4] = yIn[4] + hh*dydxt[4] ;
111 yt[3] = yIn[3] + hh*dydxt[3] ;
112 yt[2] = yIn[2] + hh*dydxt[2] ;
113 yt[1] = yIn[1] + hh*dydxt[1] ;
114 yt[0] = yIn[0] + hh*dydxt[0] ;
115 RightHandSideConst(yt,dydxm) ;
116
117 // 3rd Step K3=h*dydxm
118 // now dydxm=(K2+K3)/h
119 yt[5] = yIn[5] + h*dydxm[5] ;
120 dydxm[5] += dydxt[5] ;
121 yt[4] = yIn[4] + h*dydxm[4] ;
122 dydxm[4] += dydxt[4] ;
123 yt[3] = yIn[3] + h*dydxm[3] ;
124 dydxm[3] += dydxt[3] ;
125 yt[2] = yIn[2] + h*dydxm[2] ;
126 dydxm[2] += dydxt[2] ;
127 yt[1] = yIn[1] + h*dydxm[1] ;
128 dydxm[1] += dydxt[1] ;
129 yt[0] = yIn[0] + h*dydxm[0] ;
130 dydxm[0] += dydxt[0] ;
131 RightHandSideConst(yt,dydxt) ;
132
133 // 4th Step K4=h*dydxt
134 yOut[5] = yIn[5]+h6*(dydx[5]+dydxt[5]+2.0*dydxm[5]);
135 yOut[4] = yIn[4]+h6*(dydx[4]+dydxt[4]+2.0*dydxm[4]);
136 yOut[3] = yIn[3]+h6*(dydx[3]+dydxt[3]+2.0*dydxm[3]);
137 yOut[2] = yIn[2]+h6*(dydx[2]+dydxt[2]+2.0*dydxm[2]);
138 yOut[1] = yIn[1]+h6*(dydx[1]+dydxt[1]+2.0*dydxm[1]);
139 yOut[0] = yIn[0]+h6*(dydx[0]+dydxt[0]+2.0*dydxm[0]);
140
141} // end of DumbStepper ....................................................
void RightHandSideConst(const G4double y[], G4double dydx[]) const
Definition: G4ConstRK4.hh:96

Referenced by Stepper().

◆ GetConstField()

void G4ConstRK4::GetConstField ( const G4double  y[],
G4double  Field[] 
)
inline

Definition at line 114 of file G4ConstRK4.hh.

115{
116 G4double PositionAndTime[4];
117
118 PositionAndTime[0] = y[0];
119 PositionAndTime[1] = y[1];
120 PositionAndTime[2] = y[2];
121 // Global Time
122 PositionAndTime[3] = y[7];
123 fEq -> GetFieldValue(PositionAndTime, B) ;
124}

Referenced by Stepper().

◆ IntegratorOrder()

G4int G4ConstRK4::IntegratorOrder ( ) const
inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 76 of file G4ConstRK4.hh.

76{ return 4; }

Referenced by Stepper().

◆ RightHandSideConst()

void G4ConstRK4::RightHandSideConst ( const G4double  y[],
G4double  dydx[] 
) const
inline

Definition at line 96 of file G4ConstRK4.hh.

98{
99
100 G4double momentum_mag_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5];
101 G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
102
103 G4double cof =fEq->FCof()*inv_momentum_magnitude;
104
105 dydx[0] = y[3]*inv_momentum_magnitude; // (d/ds)x = Vx/V
106 dydx[1] = y[4]*inv_momentum_magnitude; // (d/ds)y = Vy/V
107 dydx[2] = y[5]*inv_momentum_magnitude; // (d/ds)z = Vz/V
108
109 dydx[3] = cof*(y[4]*Field[2] - y[5]*Field[1]) ; // Ax = a*(Vy*Bz - Vz*By)
110 dydx[4] = cof*(y[5]*Field[0] - y[3]*Field[2]) ; // Ay = a*(Vz*Bx - Vx*Bz)
111 dydx[5] = cof*(y[3]*Field[1] - y[4]*Field[0]) ; // Az = a*(Vx*By - Vy*Bx)
112}
G4double FCof() const
Definition: G4Mag_EqRhs.hh:83

Referenced by DumbStepper(), and Stepper().

◆ Stepper()

void G4ConstRK4::Stepper ( const G4double  y[],
const G4double  dydx[],
G4double  h,
G4double  yout[],
G4double  yerr[] 
)
virtual

Implements G4MagIntegratorStepper.

Definition at line 148 of file G4ConstRK4.cc.

153{
154 const G4int nvar = 6; // number of variables integrated
155 const G4int maxvar= GetNumberOfStateVariables();
156
157 // Correction for Richardson extrapolation
158 G4double correction = 1. / ( (1 << IntegratorOrder()) -1 );
159
160 G4int i;
161
162 // Saving yInput because yInput and yOutput can be aliases for same array
163 for (i=0; i<maxvar; i++) { yInitial[i]= yInput[i]; }
164
165 // Must copy the part of the state *not* integrated to the output
166 for (i=nvar; i<maxvar; i++) { yOutput[i]= yInput[i]; }
167
168 // yInitial[7]= yInput[7]; // The time is typically needed
169 yMiddle[7] = yInput[7]; // Copy the time from initial value
170 yOneStep[7] = yInput[7]; // As it contributes to final value of yOutput ?
171 // yOutput[7] = yInput[7]; // -> dumb stepper does it too for RK4
172 yError[7] = 0.0;
173
174 G4double halfStep = hstep * 0.5;
175
176 // Do two half steps
177 //
178 GetConstField(yInitial,Field);
179 DumbStepper (yInitial, dydx, halfStep, yMiddle);
180 RightHandSideConst(yMiddle, dydxMid);
181 DumbStepper (yMiddle, dydxMid, halfStep, yOutput);
182
183 // Store midpoint, chord calculation
184 //
185 fMidPoint = G4ThreeVector( yMiddle[0], yMiddle[1], yMiddle[2]);
186
187 // Do a full Step
188 //
189 DumbStepper(yInitial, dydx, hstep, yOneStep);
190 for(i=0;i<nvar;i++)
191 {
192 yError [i] = yOutput[i] - yOneStep[i] ;
193 yOutput[i] += yError[i]*correction ;
194 // Provides accuracy increased by 1 order via the
195 // Richardson extrapolation
196 }
197
198 fInitialPoint = G4ThreeVector( yInitial[0], yInitial[1], yInitial[2]);
199 fFinalPoint = G4ThreeVector( yOutput[0], yOutput[1], yOutput[2]);
200
201 return;
202}
CLHEP::Hep3Vector G4ThreeVector
int G4int
Definition: G4Types.hh:66
void DumbStepper(const G4double yIn[], const G4double dydx[], G4double h, G4double yOut[])
Definition: G4ConstRK4.cc:92
G4int IntegratorOrder() const
Definition: G4ConstRK4.hh:76
void GetConstField(const G4double y[], G4double Field[])
Definition: G4ConstRK4.hh:114
G4int GetNumberOfStateVariables() const

The documentation for this class was generated from the following files: