Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4TsitourasRK45 Class Reference

#include <G4TsitourasRK45.hh>

+ Inheritance diagram for G4TsitourasRK45:

Public Member Functions

 G4TsitourasRK45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 
 ~G4TsitourasRK45 ()
 
 G4TsitourasRK45 (const G4TsitourasRK45 &)=delete
 
G4TsitourasRK45 & operator= (const G4TsitourasRK45 &)=delete
 
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
 
void SetupInterpolation ()
 
void Interpolate (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
 
void interpolate (const G4double yInput[], const G4double dydx[], G4double yOut[], G4double Step, G4double tau)
 
G4double DistChord () const
 
G4int IntegratorOrder () const
 
- Public Member Functions inherited from G4MagIntegratorStepper
 G4MagIntegratorStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, G4bool isFSAL=false)
 
virtual ~G4MagIntegratorStepper ()=default
 
 G4MagIntegratorStepper (const G4MagIntegratorStepper &)=delete
 
G4MagIntegratorStepper & operator= (const G4MagIntegratorStepper &)=delete
 
virtual void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])=0
 
virtual G4double DistChord () const =0
 
void NormaliseTangentVector (G4double vec[6])
 
void NormalisePolarizationVector (G4double vec[12])
 
void RightHandSide (const G4double y[], G4double dydx[]) const
 
void RightHandSide (const G4double y[], G4double dydx[], G4double field[]) const
 
G4int GetNumberOfVariables () const
 
G4int GetNumberOfStateVariables () const
 
virtual G4int IntegratorOrder () const =0
 
G4int IntegrationOrder ()
 
G4EquationOfMotion * GetEquationOfMotion ()
 
const G4EquationOfMotion * GetEquationOfMotion () const
 
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
 
unsigned long GetfNoRHSCalls ()
 
void ResetfNORHSCalls ()
 
G4bool IsFSAL () const
 

Additional Inherited Members

- Protected Member Functions inherited from G4MagIntegratorStepper
void SetIntegrationOrder (G4int order)
 
void SetFSAL (G4bool flag=true)
 

Detailed Description

Definition at line 45 of file G4TsitourasRK45.hh.

Constructor & Destructor Documentation

◆ G4TsitourasRK45() [1/2]

G4TsitourasRK45::G4TsitourasRK45 ( G4EquationOfMotion *  EqRhs,
G4int  numberOfVariables = 6,
G4bool  primary = true 
)

Definition at line 39 of file G4TsitourasRK45.cc.

42 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
43{
44 const G4int numberOfVariables = noIntegrationVariables;
45
46 ak2 = new G4double[numberOfVariables] ;
47 ak3 = new G4double[numberOfVariables] ;
48 ak4 = new G4double[numberOfVariables] ;
49 ak5 = new G4double[numberOfVariables] ;
50 ak6 = new G4double[numberOfVariables] ;
51 ak7 = new G4double[numberOfVariables] ;
52 ak8 = new G4double[numberOfVariables] ;
53
54 // Must ensure space extra 'state' variables exists - i.e. yIn[7]
55 //
56 const G4int numStateMax = std::max(GetNumberOfStateVariables(), 8);
57 const G4int numStateVars = std::max(noIntegrationVariables,
58 numStateMax );
59
60 yTemp = new G4double[numStateVars] ;
61 yIn = new G4double[numStateVars] ;
62
63 fLastInitialVector = new G4double[numberOfVariables] ;
64 fLastFinalVector = new G4double[numberOfVariables] ;
65
66 fLastDyDx = new G4double[numberOfVariables];
67
68 fMidVector = new G4double[numberOfVariables];
69 fMidError = new G4double[numberOfVariables];
70
71 if( primary )
72 {
73 fAuxStepper = new G4TsitourasRK45(EqRhs, numberOfVariables, !primary);
74 }
75}
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4MagIntegratorStepper::GetNumberOfStateVariables
G4int GetNumberOfStateVariables() const

◆ ~G4TsitourasRK45()

G4TsitourasRK45::~G4TsitourasRK45 ( )

Definition at line 81 of file G4TsitourasRK45.cc.

82{
83 delete [] ak2;
84 delete [] ak3;
85 delete [] ak4;
86 delete [] ak5;
87 delete [] ak6;
88 delete [] ak7;
89 delete [] ak8;
90
91 delete [] yTemp;
92 delete [] yIn;
93
94 delete [] fLastInitialVector;
95 delete [] fLastFinalVector;
96 delete [] fLastDyDx;
97 delete [] fMidVector;
98 delete [] fMidError;
99
100 delete fAuxStepper;
101}

◆ G4TsitourasRK45() [2/2]

G4TsitourasRK45::G4TsitourasRK45 ( const G4TsitourasRK45 &  )
delete

Member Function Documentation

◆ DistChord()

G4double G4TsitourasRK45::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 296 of file G4TsitourasRK45.cc.

297{
298 G4double distLine, distChord;
299 G4ThreeVector initialPoint, finalPoint, midPoint;
300
301 // Store last initial and final points (they will be
302 // overwritten in self-Stepper call!)
303 //
304 initialPoint = G4ThreeVector( fLastInitialVector[0],
305 fLastInitialVector[1], fLastInitialVector[2]);
306 finalPoint = G4ThreeVector( fLastFinalVector[0],
307 fLastFinalVector[1], fLastFinalVector[2]);
308
309 // Do half a step using StepNoErr
310 //
311 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
312 fMidVector, fMidError );
313
314 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
315
316 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
317 // distance of Chord
318 //
319 if (initialPoint != finalPoint)
320 {
321 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
322 distChord = distLine;
323 }
324 else
325 {
326 distChord = (midPoint-initialPoint).mag();
327 }
328 return distChord;
329}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition: G4LineSection.hh:72
G4TsitourasRK45::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
Definition: G4TsitourasRK45.cc:116

◆ IntegratorOrder()

G4int G4TsitourasRK45::IntegratorOrder ( ) const
inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 81 of file G4TsitourasRK45.hh.

81{ return 4; }

◆ Interpolate()

void G4TsitourasRK45::Interpolate ( const G4double  yInput[],
const G4double  dydx[],
const G4double  Step,
G4double  yOut[],
G4double  tau 
)

Definition at line 249 of file G4TsitourasRK45.cc.

254{
255 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
256 // Coefficients for all the seven stages.
257
258 const G4int numberOfVariables = GetNumberOfVariables();
259
260 G4double tau0 = tau;
261
262 for(G4int i=0; i<numberOfVariables; ++i)
263 {
264 yIn[i] = yInput[i];
265 }
266
267 G4double tau_2 = tau0*tau0 ;
268 // tau_3 = tau0*tau_2,
269 // tau_4 = tau_2*tau_2;
270
271 bf1 = -1.0530884977290216*tau*(tau - 1.3299890189751412)*(tau_2 -
272 1.4364028541716351*tau + 0.7139816917074209);
273 bf2 = 0.1017*tau_2*(tau_2 - 2.1966568338249754*tau +
274 1.2949852507374631);
275 bf3 = 2.490627285651252793*tau_2*(tau_2 - 2.38535645472061657*tau
276 + 1.57803468208092486);
277 bf4 = -16.54810288924490272*(tau - 1.21712927295533244)*
278 (tau - 0.61620406037800089)*tau_2;
279 bf5 = 47.37952196281928122*(tau - 1.203071208372362603)*
280 (tau - 0.658047292653547382)*tau_2;
281 bf6 = -34.87065786149660974*(tau - 1.2)*(tau -
282 0.666666666666666667)*tau_2;
283 bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2;
284
285 // Putting together the coefficients calculated as the respective
286 // stage coefficients
287 //
288 for(G4int i=0; i<numberOfVariables; ++i)
289 {
290 yOut[i] = yIn[i] + Step*( bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i]
291 + bf4*ak4[i] + bf5*ak5[i] + bf6*ak6[i]
292 + bf7*ak7[i] ) ;
293 }
294}
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const

◆ interpolate()

void G4TsitourasRK45::interpolate ( const G4double  yInput[],
const G4double  dydx[],
G4double  yOut[],
G4double  Step,
G4double  tau 
)

◆ operator=()

G4TsitourasRK45 & G4TsitourasRK45::operator= ( const G4TsitourasRK45 &  )
delete

◆ SetupInterpolation()

void G4TsitourasRK45::SetupInterpolation ( )

Definition at line 243 of file G4TsitourasRK45.cc.

245{
246 // Nothing to be done
247}

◆ Stepper()

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

Implements G4MagIntegratorStepper.

Definition at line 116 of file G4TsitourasRK45.cc.

121{
122 const G4double b21 = 0.161 ,
123 b31 = -0.00848065549235698854 ,
124 b32 = 0.335480655492356989 ,
125
126 b41 = 2.89715305710549343 ,
127 b42 = -6.35944848997507484 ,
128 b43 = 4.36229543286958141 ,
129
130 b51 = 5.325864828439257,
131 b52 = -11.748883564062828,
132 b53 = 7.49553934288983621 ,
133 b54 = -0.09249506636175525,
134
135 b61 = 5.8614554429464200,
136 b62 = -12.9209693178471093 ,
137 b63 = 8.1593678985761586 ,
138 b64 = -0.071584973281400997,
139 b65 = -0.0282690503940683829,
140
141 b71 = 0.0964607668180652295 ,
142 b72 = 0.01,
143 b73 = 0.479889650414499575,
144 b74 = 1.37900857410374189,
145 b75 = -3.2900695154360807,
146 b76 = 2.32471052409977398,
147
148 // c1 = 0.001780011052226 ,
149 // c2 = 0.000816434459657 ,
150 // c3 = -0.007880878010262 ,
151 // c4 = 0.144711007173263 ,
152 // c5 = -0.582357165452555 ,
153 // c6 = 0.458082105929187 ,
154 // c7 = 1.0/66.0 ;
155
156 dc1 = 0.0935237485818927066 - b71 , // - 0.001780011052226,
157 dc2 = 0.00865288314156636761 - b72, // - 0.000816434459657,
158 dc3 = 0.492893099131431868 - b73 , // + 0.007880878010262,
159 dc4 = 1.14023541226785810 - b74 , // 0.144711007173263,
160 dc5 = - 2.3291801924393646 - b75, // + 0.582357165452555,
161 dc6 = 1.56887504931661552 - b76 , // - 0.458082105929187,
162 dc7 = 0.025; //- 1.0/66.0 ;
163
164 // dc1 = -3.0/1280.0,
165 // dc2 = 0.0,
166 // dc3 = 6561.0/632320.0,
167 // dc4 = -343.0/20800.0,
168 // dc5 = 243.0/12800.0,
169 // dc6 = -1.0/95.0,
170 // dc7 = 0.0 ;
171
172 const G4int numberOfVariables = GetNumberOfVariables();
173
174 // The number of variables to be integrated over
175 //
176 yOut[7] = yTemp[7] = yIn[7] = yInput[7];
177
178 // Saving yInput because yInput and yOut can be aliases for same array
179 //
180 for(G4int i=0; i<numberOfVariables; ++i)
181 {
182 yIn[i]=yInput[i];
183 }
184 // RightHandSide(yIn, dydx) ; // 1st Step - Not doing, getting passed
185
186 for(G4int i=0; i<numberOfVariables; ++i)
187 {
188 yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
189 }
190 RightHandSide(yTemp, ak2) ; // 2nd Stage
191
192 for(G4int i=0; i<numberOfVariables; ++i)
193 {
194 yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
195 }
196 RightHandSide(yTemp, ak3) ; // 3rd Stage
197
198 for(G4int i=0; i<numberOfVariables; ++i)
199 {
200 yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
201 }
202 RightHandSide(yTemp, ak4) ; // 4th Stage
203
204 for(G4int i=0; i<numberOfVariables; ++i)
205 {
206 yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
207 b54*ak4[i]) ;
208 }
209 RightHandSide(yTemp, ak5) ; // 5th Stage
210
211 for(G4int i=0; i<numberOfVariables; ++i)
212 {
213 yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
214 b64*ak4[i] + b65*ak5[i]) ;
215 }
216 RightHandSide(yTemp, ak6) ; // 6th Stage
217
218 for(G4int i=0; i<numberOfVariables; ++i)
219 {
220 yOut[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
221 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
222 }
223 RightHandSide(yOut, ak7); // 7th Stage
224
225 //Calculate the error in the step:
226 for(G4int i=0; i<numberOfVariables; ++i)
227 {
228 yErr[i] = Step*(dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
229 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ;
230
231 // Store Input and Final values, for possible use in calculating chord
232 //
233 fLastInitialVector[i] = yIn[i] ;
234 fLastFinalVector[i] = yOut[i];
235 fLastDyDx[i] = dydx[i];
236 }
237
238 fLastStepLength = Step;
239
240 return ;
241}
G4MagIntegratorStepper::RightHandSide
void RightHandSide(const G4double y[], G4double dydx[]) const

Referenced by DistChord().

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