Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TsitourasRK45 Class Reference

#include <G4TsitourasRK45.hh>

Inheritance diagram for G4TsitourasRK45:
Collaboration diagram for G4TsitourasRK45:

Public Member Functions

 G4TsitourasRK45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 
 ~G4TsitourasRK45 ()
 
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, bool isFSAL=false)
 
virtual ~G4MagIntegratorStepper ()
 
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
 
G4int IntegrationOrder ()
 
G4EquationOfMotionGetEquationOfMotion ()
 
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
 
unsigned long GetfNoRHSCalls ()
 
void ResetfNORHSCalls ()
 
bool IsFSAL ()
 

Additional Inherited Members

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

Detailed Description

Definition at line 39 of file G4TsitourasRK45.hh.

Constructor & Destructor Documentation

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

Definition at line 51 of file G4TsitourasRK45.cc.

54  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
55  fLastStepLength(0.), fAuxStepper(0)
56 {
57  const G4int numberOfVariables = noIntegrationVariables;
58  G4cout << "G4TsitourasRK45 constructor called." << G4endl;
59 
60  ak2 = new G4double[numberOfVariables] ;
61  ak3 = new G4double[numberOfVariables] ;
62  ak4 = new G4double[numberOfVariables] ;
63  ak5 = new G4double[numberOfVariables] ;
64  ak6 = new G4double[numberOfVariables] ;
65  ak7 = new G4double[numberOfVariables] ;
66  ak8 = new G4double[numberOfVariables] ;
67 
68 
69  // Must ensure space extra 'state' variables exists - i.e. yIn[7]
70  const G4int numStateMax = std::max(GetNumberOfStateVariables(), 8);
71  const G4int numStateVars = std::max(noIntegrationVariables,
72  numStateMax );
73  // GetNumberOfStateVariables() );
74 
75  yTemp = new G4double[numStateVars] ;
76  yIn = new G4double[numStateVars] ;
77 
78  fLastInitialVector = new G4double[numberOfVariables] ;
79  fLastFinalVector = new G4double[numberOfVariables] ;
80 
81  fLastDyDx = new G4double[numberOfVariables];
82 
83  fMidVector = new G4double[numberOfVariables];
84  fMidError = new G4double[numberOfVariables];
85  if( primary )
86  {
87  fAuxStepper = new G4TsitourasRK45(EqRhs, numberOfVariables, !primary);
88  }
89 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4MagIntegratorStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, bool isFSAL=false)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4int GetNumberOfStateVariables() const
#define G4endl
Definition: G4ios.hh:61
G4TsitourasRK45(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4TsitourasRK45::~G4TsitourasRK45 ( )

Definition at line 95 of file G4TsitourasRK45.cc.

96 {
97  delete[] ak2;
98  delete[] ak3;
99  delete[] ak4;
100  delete[] ak5;
101  delete[] ak6;
102  delete[] ak7;
103  delete[] ak8;
104 
105  delete[] yTemp;
106  delete[] yIn;
107 
108  delete[] fLastInitialVector;
109  delete[] fLastFinalVector;
110  delete[] fLastDyDx;
111  delete[] fMidVector;
112  delete[] fMidError;
113 
114  delete fAuxStepper;
115 }

Member Function Documentation

G4double G4TsitourasRK45::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 309 of file G4TsitourasRK45.cc.

310 {
311  G4double distLine, distChord;
312  G4ThreeVector initialPoint, finalPoint, midPoint;
313 
314  // Store last initial and final points (they will be overwritten in self-Stepper call!)
315  initialPoint = G4ThreeVector( fLastInitialVector[0],
316  fLastInitialVector[1], fLastInitialVector[2]);
317  finalPoint = G4ThreeVector( fLastFinalVector[0],
318  fLastFinalVector[1], fLastFinalVector[2]);
319 
320  // Do half a step using StepNoErr
321 
322  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
323  fMidVector, fMidError );
324 
325  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
326 
327  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
328  // distance of Chord
329 
330 
331  if (initialPoint != finalPoint)
332  {
333  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
334  distChord = distLine;
335  }
336  else
337  {
338  distChord = (midPoint-initialPoint).mag();
339  }
340  return distChord;
341 }
CLHEP::Hep3Vector G4ThreeVector
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4int G4TsitourasRK45::IntegratorOrder ( ) const
inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 72 of file G4TsitourasRK45.hh.

72 {return 4; }
void G4TsitourasRK45::Interpolate ( const G4double  yInput[],
const G4double  dydx[],
const G4double  Step,
G4double  yOut[],
G4double  tau 
)

Definition at line 266 of file G4TsitourasRK45.cc.

266  {
267 
268 
269  G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
270  // Coefficients for all the seven stages.
271 
272  const G4int numberOfVariables= this->GetNumberOfVariables();
273 
274  G4double tau0 = tau;
275 
276  for(int i=0;i<numberOfVariables;i++)
277  {
278  yIn[i]=yInput[i];
279  }
280 
281  G4double
282  tau_2 = tau0*tau0 ;
283 // tau_3 = tau0*tau_2,
284 // tau_4 = tau_2*tau_2;
285 
286  bf1 = -1.0530884977290216*tau*(tau - 1.3299890189751412)*(tau_2 -
287  1.4364028541716351*tau + 0.7139816917074209),
288  bf2 = 0.1017*tau_2*(tau_2 - 2.1966568338249754*tau +
289  1.2949852507374631),
290  bf3 = 2.490627285651252793*tau_2*(tau_2 - 2.38535645472061657*tau
291  + 1.57803468208092486) ,
292  bf4 = -16.54810288924490272*(tau - 1.21712927295533244)*
293  (tau - 0.61620406037800089)*tau_2,
294  bf5 = 47.37952196281928122*(tau - 1.203071208372362603)*
295  (tau - 0.658047292653547382)*tau_2,
296  bf6 = -34.87065786149660974*(tau - 1.2)*(tau -
297  0.666666666666666667)*tau_2,
298  bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2;
299 
300  //Putting together the coefficients calculated as the respective stage coefficients
301  for( int i=0; i<numberOfVariables; i++){
302  yOut[i] = yIn[i] + Step*( bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] + bf4*ak4[i]
303  + bf5*ak5[i] + bf6*ak6[i] + bf7*ak7[i] ) ;
304  }
305 }
int G4int
Definition: G4Types.hh:78
G4int GetNumberOfVariables() const
Definition: Step.hh:41
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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

Definition at line 260 of file G4TsitourasRK45.cc.

261 {
262  //Nothing to be done
263 }
void G4TsitourasRK45::Stepper ( const G4double  y[],
const G4double  dydx[],
G4double  h,
G4double  yout[],
G4double  yerr[] 
)
virtual

Implements G4MagIntegratorStepper.

Definition at line 129 of file G4TsitourasRK45.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:


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