Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4DormandPrince745 Class Reference

#include <G4DormandPrince745.hh>

Inheritance diagram for G4DormandPrince745:
Inheritance graph
[legend]
Collaboration diagram for G4DormandPrince745:
Collaboration graph
[legend]

Public Member Functions

 G4DormandPrince745 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 
 ~G4DormandPrince745 ()
 
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
 
void SetupInterpolation_low ()
 
void Interpolate_low (G4double yOut[], G4double tau)
 
void SetupInterpolation ()
 
void Interpolate (G4double tau, G4double yOut[])
 
void SetupInterpolation_high ()
 
void Interpolate_high (G4double yOut[], G4double tau)
 
G4double DistChord () const
 
G4double DistChord2 () const
 
G4double DistChord3 () const
 
G4double DistLine (G4double yStart[], G4double yMid[], G4double yEnd[]) 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 ()
 
G4EquationOfMotion * GetEquationOfMotion ()
 
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 52 of file G4DormandPrince745.hh.

Constructor & Destructor Documentation

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

Definition at line 71 of file G4DormandPrince745.cc.

74 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
75  fAuxStepper(0)
76 {
77  const G4int numberOfVariables = // noIntegrationVariables;
78  std::max( noIntegrationVariables,
79  ( ( (noIntegrationVariables-1)/4 + 1 ) * 4 ) );
80  // For better alignment with cache-line
81 
82  //New Chunk of memory being created for use by the stepper
83 
84  //ak_i - for storing intermediate RHS
85  ak2 = new G4double[numberOfVariables];
86  ak3 = new G4double[numberOfVariables];
87  ak4 = new G4double[numberOfVariables];
88  ak5 = new G4double[numberOfVariables];
89  ak6 = new G4double[numberOfVariables];
90  ak7 = new G4double[numberOfVariables];
91  // Also always allocate arrays for interpolation stages
92  ak8 = new G4double[numberOfVariables];
93  ak9 = new G4double[numberOfVariables];
94 
95  // Must ensure space for extra 'state' variables exists - i.e. yIn[7]
96  const G4int numStateVars =
97  std::max(noIntegrationVariables,
98  std::max( GetNumberOfStateVariables(), 8)
99  );
100  yTemp = new G4double[numStateVars];
101  yIn = new G4double[numStateVars];
102 
103  fLastInitialVector = new G4double[numStateVars] ;
104  fLastFinalVector = new G4double[numStateVars] ;
105 
106  // fLastDyDx = new G4double[numberOfVariables];
107 
108  fMidVector = new G4double[numStateVars];
109  fMidError = new G4double[numStateVars];
110 
111  yTemp = new G4double[numberOfVariables] ;
112  yIn = new G4double[numberOfVariables] ;
113 
114  fLastInitialVector = new G4double[numberOfVariables] ;
115  fLastFinalVector = new G4double[numberOfVariables] ;
116  fInitialDyDx = new G4double[numberOfVariables];
117 
118  fMidVector = new G4double[numberOfVariables];
119  fMidError = new G4double[numberOfVariables];
120  if( primary )
121  {
122  fAuxStepper = new G4DormandPrince745(EqRhs, numberOfVariables,
123  !primary);
124  }
125 }
G4DormandPrince745::G4DormandPrince745
G4DormandPrince745(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition: G4DormandPrince745.cc:71
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::G4MagIntegratorStepper
G4MagIntegratorStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, bool isFSAL=false)
Definition: G4MagIntegratorStepper.cc:36
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4MagIntegratorStepper::GetNumberOfStateVariables
G4int GetNumberOfStateVariables() const
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4DormandPrince745::~G4DormandPrince745 ( )

Definition at line 128 of file G4DormandPrince745.cc.

129 {
130  //clear all previously allocated memory for stepper and DistChord
131  delete[] ak2;
132  delete[] ak3;
133  delete[] ak4;
134  delete[] ak5;
135  delete[] ak6;
136  delete[] ak7;
137  // Used only for interpolation
138  delete[] ak8;
139  delete[] ak9;
140 
141  delete[] yTemp;
142  delete[] yIn;
143 
144  delete[] fLastInitialVector;
145  delete[] fLastFinalVector;
146  delete[] fInitialDyDx;
147  delete[] fMidVector;
148  delete[] fMidError;
149 
150  delete fAuxStepper;
151 }

Member Function Documentation

G4double G4DormandPrince745::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 335 of file G4DormandPrince745.cc.

336 {
337  //Coefficients for halfway interpolation
338  const G4double
339  hf1 = 5783653.0/57600000.0 ,
340  hf2 = 0. ,
341  hf3 = 466123.0/1192500.0 ,
342  hf4 = -41347.0/1920000.0 ,
343  hf5 = 16122321.0/339200000.0 ,
344  hf6 = -7117.0/20000.0,
345  hf7 = 183.0/10000.0 ;
346 
347  for(int i=0; i<3; i++){
348  fMidVector[i] = fLastInitialVector[i] + fLastStepLength*(
349  hf1*fInitialDyDx[i] + hf2*ak2[i] + hf3*ak3[i] + hf4*ak4[i] +
350  hf5*ak5[i] + hf6*ak6[i] + hf7*ak7[i] );
351  }
352 
353  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
354  // distance of Chord
355 
356  return DistLine( fLastInitialVector, fMidVector, fLastFinalVector);
357 }
G4DormandPrince745::DistLine
G4double DistLine(G4double yStart[], G4double yMid[], G4double yEnd[]) const
Definition: G4DormandPrince745.cc:298
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4DormandPrince745::DistChord2 ( ) const

Definition at line 322 of file G4DormandPrince745.cc.

323 {
324  // Copy the values of stages from this (original) into the Aux Stepper
325  *fAuxStepper = *this;
326 
327  //Preparing for the interpolation
328  fAuxStepper->SetupInterpolation(); // (fLastInitialVector, fInitialDyDx, fLastStepLength);
329  //Interpolate to half step
330  fAuxStepper->Interpolate( /*fLastInitialVector, fInitialDyDx, fLastStepLength,*/ 0.5, fAuxStepper->fMidVector);
331 
332  return DistLine( fLastInitialVector, fAuxStepper->fMidVector, fLastFinalVector);
333 }
G4DormandPrince745::Interpolate
void Interpolate(G4double tau, G4double yOut[])
Definition: G4DormandPrince745.hh:88
G4DormandPrince745::DistLine
G4double DistLine(G4double yStart[], G4double yMid[], G4double yEnd[]) const
Definition: G4DormandPrince745.cc:298

Here is the call graph for this function:

G4double G4DormandPrince745::DistChord3 ( ) const

Definition at line 360 of file G4DormandPrince745.cc.

361 {
362  // Do half a step using StepNoErr
363  fAuxStepper->Stepper( fLastInitialVector, fInitialDyDx, 0.5 * fLastStepLength,
364  fAuxStepper->fMidVector, fAuxStepper->fMidError) ;
365  return DistLine( fLastInitialVector, fAuxStepper->fMidVector, fLastFinalVector);
366 }
G4DormandPrince745::DistLine
G4double DistLine(G4double yStart[], G4double yMid[], G4double yEnd[]) const
Definition: G4DormandPrince745.cc:298
G4DormandPrince745::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
Definition: G4DormandPrince745.cc:180

Here is the call graph for this function:

G4double G4DormandPrince745::DistLine ( G4double  yStart[],
G4double  yMid[],
G4double  yEnd[] 
) const

Definition at line 298 of file G4DormandPrince745.cc.

299 {
300  G4double distLine, distChord;
301  G4ThreeVector initialPoint, finalPoint, midPoint;
302 
303  initialPoint = G4ThreeVector( yStart[0], yStart[1], yStart[2]);
304  finalPoint = G4ThreeVector( yEnd[0], yEnd[1], yEnd[2]);
305  midPoint = G4ThreeVector( yMid[0], yMid[1], yMid[2]);
306 
307  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
308  // distance of Chord
309  if (initialPoint != finalPoint)
310  {
311  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
312  distChord = distLine;
313  }
314  else
315  {
316  distChord = (midPoint-initialPoint).mag();
317  }
318  return distChord;
319 }
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition: G4LineSection.hh:84
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4int G4DormandPrince745::IntegratorOrder ( ) const
inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 117 of file G4DormandPrince745.hh.

117 {return 4; }
void G4DormandPrince745::Interpolate ( G4double  tau,
G4double  yOut[] 
)
inline

Definition at line 88 of file G4DormandPrince745.hh.

95  {
96  Interpolate_low( /* yInput, dydx, Step, */ yOut, tau);
97  // Interpolate_high( /* yInput, dydx, Step, */ yOut, tau);
98  }
G4DormandPrince745::Interpolate_low
void Interpolate_low(G4double yOut[], G4double tau)
Definition: G4DormandPrince745.cc:381

Here is the call graph for this function:

Here is the caller graph for this function:

void G4DormandPrince745::Interpolate_high ( G4double  yOut[],
G4double  tau 
)

Definition at line 478 of file G4DormandPrince745.cc.

482  {
483  //Define the coefficients for the polynomials
484  G4double bi[10][5], b[10];
485  const G4int numberOfVariables = this->GetNumberOfVariables();
486  const G4double *dydx = fInitialDyDx;
487  // const G4double fullStep = fLastStepLength;
488 
489  // If given requestedStep in argument:
490  // G4double tau = requestedStep / fLastStepLength;
491 
492  // COEFFICIENTS OF bi[1]
493  bi[1][0] = 1.0 ,
494  bi[1][1] = -38039.0/7040.0 ,
495  bi[1][2] = 125923.0/10560.0 ,
496  bi[1][3] = -19683.0/1760.0 ,
497  bi[1][4] = 3303.0/880.0 ,
498  // --------------------------------------------------------
499  //
500  // COEFFICIENTS OF bi[2]
501  bi[2][0] = 0.0 ,
502  bi[2][1] = 0.0 ,
503  bi[2][2] = 0.0 ,
504  bi[2][3] = 0.0 ,
505  bi[2][4] = 0.0 ,
506  // --------------------------------------------------------
507  //
508  // COEFFICIENTS OF bi[3]
509  bi[3][0] = 0.0 ,
510  bi[3][1] = -12500.0/4081.0 ,
511  bi[3][2] = 205000.0/12243.0 ,
512  bi[3][3] = -90000.0/4081.0 ,
513  bi[3][4] = 36000.0/4081.0 ,
514  // --------------------------------------------------------
515  //
516  // COEFFICIENTS OF bi[4]
517  bi[4][0] = 0.0 ,
518  bi[4][1] = -3125.0/704.0 ,
519  bi[4][2] = 25625.0/1056.0 ,
520  bi[4][3] = -5625.0/176.0 ,
521  bi[4][4] = 1125.0/88.0 ,
522  // --------------------------------------------------------
523  //
524  // COEFFICIENTS OF bi[5]
525  bi[5][0] = 0.0 ,
526  bi[5][1] = 164025.0/74624.0 ,
527  bi[5][2] = -448335.0/37312.0 ,
528  bi[5][3] = 295245.0/18656.0 ,
529  bi[5][4] = -59049.0/9328.0 ,
530  // --------------------------------------------------------
531  //
532  // COEFFICIENTS OF bi[6]
533  bi[6][0] = 0.0 ,
534  bi[6][1] = -25.0/28.0 ,
535  bi[6][2] = 205.0/42.0 ,
536  bi[6][3] = -45.0/7.0 ,
537  bi[6][4] = 18.0/7.0 ,
538  // --------------------------------------------------------
539  //
540  // COEFFICIENTS OF bi[7]
541  bi[7][0] = 0.0 ,
542  bi[7][1] = -2.0/11.0 ,
543  bi[7][2] = 73.0/55.0 ,
544  bi[7][3] = -171.0/55.0 ,
545  bi[7][4] = 108.0/55.0 ,
546  // --------------------------------------------------------
547  //
548  // COEFFICIENTS OF bi[8]
549  bi[8][0] = 0.0 ,
550  bi[8][1] = 189.0/22.0 ,
551  bi[8][2] = -1593.0/55.0 ,
552  bi[8][3] = 3537.0/110.0 ,
553  bi[8][4] = -648.0/55.0 ,
554  // --------------------------------------------------------
555  //
556  // COEFFICIENTS OF bi[9]
557  bi[9][0] = 0.0 ,
558  bi[9][1] = 351.0/110.0 ,
559  bi[9][2] = -999.0/55.0 ,
560  bi[9][3] = 2943.0/110.0 ,
561  bi[9][4] = -648.0/55.0 ;
562  // --------------------------------------------------------
563 
564  // for(G4int i = 0; i< numberOfVariables; i++) { yIn[i] = yInput[i]; }
565 
566  // Calculating the polynomials :
567 #if 1
568  for(int iStage=1; iStage<=9; iStage++){
569  b[iStage] = 0;
570  }
571 
572  for(int j=0; j<=4; j++){
573  G4double tauPower = 1.0;
574  for(int iStage=1; iStage<=9; iStage++){
575  b[iStage] += bi[iStage][j]*tauPower;
576  }
577  tauPower *= tau;
578  }
579 #else
580  G4double tau0 = tau;
581 
582  for(int i=1; i<=9; i++){ //Here i is NOT the coordinate no. , it's stage no.
583  b[i] = 0;
584  tau = 1.0;
585  for(int j=0; j<=4; j++){
586  b[i] += bi[i][j]*tau;
587  tau*=tau0;
588  }
589  }
590 #endif
591 
592  G4double stepLen = fLastStepLength * tau;
593  for(int i=0; i<numberOfVariables; i++){ //Here i IS the cooridnate no.
594  yOut[i] = yIn[i] + stepLen *(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
595  b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
596  b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] );
597  }
598 
599 }
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4DormandPrince745::Interpolate_low ( G4double  yOut[],
G4double  tau 
)

Definition at line 381 of file G4DormandPrince745.cc.

386 {
387  G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
388  // Coefficients for all the seven stages.
389  G4double Step = fLastStepLength;
390  const G4double *dydx= fInitialDyDx;
391 
392  const G4int numberOfVariables= this->GetNumberOfVariables();
393 
394  // for(int i=0;i<numberOfVariables;i++) { yIn[i]=yInput[i]; }
395 
396  const G4double
397  tau_2 = tau * tau,
398  tau_3 = tau * tau_2,
399  tau_4 = tau_2 * tau_2;
400 
401  bf1 = (157015080.0*tau_4 - 13107642775.0*tau_3+ 34969693132.0*tau_2- 32272833064.0*tau
402  + 11282082432.0)/11282082432.0,
403  bf2 = 0.0 ,
404  bf3 = - 100.0*tau*(15701508.0*tau_3 - 914128567.0*tau_2 + 2074956840.0*tau
405  - 1323431896.0)/32700410799.0,
406  bf4 = 25.0*tau*(94209048.0*tau_3- 1518414297.0*tau_2+ 2460397220.0*tau - 889289856.0)/5641041216.0 ,
407  bf5 = -2187.0*tau*(52338360.0*tau_3 - 451824525.0*tau_2 + 687873124.0*tau - 259006536.0)/199316789632.0 ,
408  bf6 = 11.0*tau*(106151040.0*tau_3- 661884105.0*tau_2 + 946554244.0*tau - 361440756.0)/2467955532.0 ,
409  bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 82437520.0*tau + 44764047.0)/ 29380423.0 ;
410 
411  //Putting together the coefficients calculated as the respective stage coefficients
412  for( int i=0; i<numberOfVariables; i++){
413  yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] + bf4*ak4[i]
414  + bf5*ak5[i] + bf6*ak6[i] + bf7*ak7[i] ) ;
415  }
416 }
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4DormandPrince745::SetupInterpolation ( )
inline

Definition at line 78 of file G4DormandPrince745.hh.

82  {
83  SetupInterpolation_low( /* yInput, dydx, Step */ );
84  // SetupInterpolation_high( /* yInput, dydx, Step */ );
85  }

Here is the call graph for this function:

Here is the caller graph for this function:

void G4DormandPrince745::SetupInterpolation_high ( )

Definition at line 427 of file G4DormandPrince745.cc.

429  {
430 
431  //Coefficients for the additional stages :
432  const G4double
433  b81 = 6245.0/62208.0 ,
434  b82 = 0.0 ,
435  b83 = 8875.0/103032.0 ,
436  b84 = -125.0/1728.0 ,
437  b85 = 801.0/13568.0 ,
438  b86 = -13519.0/368064.0 ,
439  b87 = 11105.0/368064.0 ,
440 
441  b91 = 632855.0/4478976.0 ,
442  b92 = 0.0 ,
443  b93 = 4146875.0/6491016.0 ,
444  b94 = 5490625.0/14183424.0 ,
445  b95 = -15975.0/108544.0 ,
446  b96 = 8295925.0/220286304.0 ,
447  b97 = -1779595.0/62938944.0 ,
448  b98 = -805.0/4104.0 ;
449 
450  const G4int numberOfVariables= this->GetNumberOfVariables();
451  const G4double *dydx = fInitialDyDx;
452  const G4double Step = fLastStepLength;
453 
454  // Saving yInput because yInput and yOut can be aliases for same array
455  // for(int i=0;i<numberOfVariables;i++) { yIn[i]=yInput[i]; }
456  // yTemp[7] = yIn[7];
457 
458  //Evaluate the extra stages :
459  for(int i=0;i<numberOfVariables;i++)
460  {
461  yTemp[i] = yIn[i] + Step*(b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
462  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
463  b87*ak7[i]);
464  }
465  RightHandSide(yTemp, ak8); //8th Stage
466 
467  for(int i=0;i<numberOfVariables;i++)
468  {
469  yTemp[i] = yIn[i] + Step*(b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
470  b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
471  b97*ak7[i] + b98*ak8[i] );
472  }
473  RightHandSide(yTemp, ak9); //9th Stage
474 }
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
G4MagIntegratorStepper::RightHandSide
void RightHandSide(const double y[], double dydx[])
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4DormandPrince745::SetupInterpolation_low ( )

Definition at line 376 of file G4DormandPrince745.cc.

377 {
378  //Nothing to be done
379 }

Here is the caller graph for this function:

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

Implements G4MagIntegratorStepper.

Definition at line 180 of file G4DormandPrince745.cc.

185 {
186  G4int i;
187 
188  //The various constants defined on the basis of butcher tableu
189  const G4double //G4double - only once
190  b21 = 0.2 ,
191 
192  b31 = 3.0/40.0, b32 = 9.0/40.0 ,
193 
194  b41 = 44.0/45.0, b42 = -56.0/15.0, b43 = 32.0/9.0,
195 
196  b51 = 19372.0/6561.0, b52 = -25360.0/2187.0, b53 = 64448.0/6561.0,
197  b54 = -212.0/729.0 ,
198 
199  b61 = 9017.0/3168.0 , b62 = -355.0/33.0,
200  b63 = 46732.0/5247.0 , b64 = 49.0/176.0 ,
201  b65 = -5103.0/18656.0 ,
202 
203  b71 = 35.0/384.0, b72 = 0.,
204  b73 = 500.0/1113.0, b74 = 125.0/192.0,
205  b75 = -2187.0/6784.0, b76 = 11.0/84.0,
206 
207  //Sum of columns, sum(bij) = ei
208 // e1 = 0. ,
209 // e2 = 1.0/5.0 ,
210 // e3 = 3.0/10.0 ,
211 // e4 = 4.0/5.0 ,
212 // e5 = 8.0/9.0 ,
213 // e6 = 1.0 ,
214 // e7 = 1.0 ,
215 
216 // Difference between the higher and the lower order method coeff. :
217  // b7j are the coefficients of higher order
218 
219  dc1 = -( b71 - 5179.0/57600.0),
220  dc2 = -( b72 - .0),
221  dc3 = -( b73 - 7571.0/16695.0),
222  dc4 = -( b74 - 393.0/640.0),
223  dc5 = -( b75 + 92097.0/339200.0),
224  dc6 = -( b76 - 187.0/2100.0),
225  dc7 = -( - 1.0/40.0 ); //end of declaration
226 
227  const G4int numberOfVariables= this->GetNumberOfVariables();
228 
229  // The number of variables to be integrated over
230  yOut[7] = yTemp[7] = yInput[7];
231  // Saving yInput because yInput and yOut can be aliases for same array
232 
233  for(i=0;i<numberOfVariables;i++)
234  {
235  yIn[i]=yInput[i];
236  }
237 
238  // RightHandSide(yIn, DyDx) ;
239  // 1st Step - Not doing, getting passed
240 
241  for(i=0;i<numberOfVariables;i++)
242  {
243  yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
244  }
245  RightHandSide(yTemp, ak2) ; // 2nd stage
246 
247  for(i=0;i<numberOfVariables;i++)
248  {
249  yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
250  }
251  RightHandSide(yTemp, ak3) ; // 3rd stage
252 
253  for(i=0;i<numberOfVariables;i++)
254  {
255  yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
256  }
257  RightHandSide(yTemp, ak4) ; // 4th stage
258 
259  for(i=0;i<numberOfVariables;i++)
260  {
261  yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
262  b54*ak4[i]) ;
263  }
264  RightHandSide(yTemp, ak5) ; // 5th stage
265 
266  for(i=0;i<numberOfVariables;i++)
267  {
268  yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
269  b64*ak4[i] + b65*ak5[i]) ;
270  }
271  RightHandSide(yTemp, ak6) ; // 6th stage
272 
273  for(i=0;i<numberOfVariables;i++)
274  {
275  yOut[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
276  b74*ak4[i] + b75*ak5[i] + b76*ak6[i] );
277  }
278  RightHandSide(yOut, ak7); //7th and Final stage
279 
280  for(i=0;i<numberOfVariables;i++)
281  {
282  yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
283  dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) + 1.5e-18 ;
284 
285  // Store Input and Final values, for possible use in calculating chord
286  fLastInitialVector[i] = yIn[i] ;
287  fLastFinalVector[i] = yOut[i];
288  fInitialDyDx[i] = DyDx[i];
289  }
290 
291  fLastStepLength = Step;
292 
293  return ;
294 }
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
G4MagIntegratorStepper::RightHandSide
void RightHandSide(const double y[], double dydx[])
G4double
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: