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

#include <G4BogackiShampine45.hh>

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

Public Member Functions

 G4BogackiShampine45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 
 ~G4BogackiShampine45 ()
 
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
 
void SetupInterpolation ()
 
void Interpolate (G4double tau, G4double yOut[])
 
void SetupInterpolationHigh ()
 
void InterpolateHigh (G4double tau, G4double yOut[]) const
 
G4double DistChord () const
 
G4int IntegratorOrder () const
 
void GetLastDydx (G4double dyDxLast[])
 
void PrepareConstants ()
 
- 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 51 of file G4BogackiShampine45.hh.

Constructor & Destructor Documentation

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

Definition at line 76 of file G4BogackiShampine45.cc.

79  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
80  fAuxStepper(0),
81  fPreparedInterpolation(false)
82 {
83 
84  const G4int numberOfVariables = noIntegrationVariables;
85 
86  //New Chunk of memory being created for use by the stepper
87 
88  //aki - for storing intermediate RHS
89  ak2 = new G4double[numberOfVariables];
90  ak3 = new G4double[numberOfVariables];
91  ak4 = new G4double[numberOfVariables];
92  ak5 = new G4double[numberOfVariables];
93  ak6 = new G4double[numberOfVariables];
94  ak7 = new G4double[numberOfVariables];
95  ak8 = new G4double[numberOfVariables];
96  ak9 = new G4double[numberOfVariables];
97  ak10 = new G4double[numberOfVariables];
98  ak11 = new G4double[numberOfVariables];
99 
100  for (int i = 0; i < 6; i++) {
101  p[i]= new G4double[numberOfVariables];
102  }
103 
104  assert ( GetNumberOfStateVariables() >= 8 );
105  const G4int numStateVars = std::max(noIntegrationVariables,
106  GetNumberOfStateVariables() );
107 
108  // Must ensure space extra 'state' variables exists - i.e. yIn[7]
109  yTemp = new G4double[numStateVars];
110  yIn = new G4double[numStateVars] ;
111 
112  fLastInitialVector = new G4double[numStateVars] ;
113  fLastFinalVector = new G4double[numStateVars] ;
114  fLastDyDx = new G4double[numberOfVariables]; // Only derivatives
115 
116  fMidVector = new G4double[numberOfVariables]; // new G4double[numStateVars];
117  fMidError = new G4double[numberOfVariables]; // new G4double[numStateVars];
118 
119  if( ! fPreparedConstants )
120  PrepareConstants();
121 
122  if( primary )
123  {
124  fAuxStepper = new G4BogackiShampine45(EqRhs, numberOfVariables, false);
125  }
126 }
G4BogackiShampine45::G4BogackiShampine45
G4BogackiShampine45(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition: G4BogackiShampine45.cc:76
ak2
static const G4double ak2
Definition: G4MuPairProductionModel.cc:94
p
const char * p
Definition: xmltok.h:285
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:

G4BogackiShampine45::~G4BogackiShampine45 ( )

Definition at line 130 of file G4BogackiShampine45.cc.

130  {
131  //clear all previously allocated memory for stepper and DistChord
132  delete[] ak2;
133  delete[] ak3;
134  delete[] ak4;
135  delete[] ak5;
136  delete[] ak6;
137  delete[] ak7;
138  delete[] ak8;
139  delete[] ak9;
140  delete[] ak10;
141  delete[] ak11;
142 
143  for (int i = 0; i < 6; i++) {
144  delete[] p[i];
145  }
146 
147  delete[] yTemp;
148  delete[] yIn;
149 
150  delete[] fLastInitialVector;
151  delete[] fLastFinalVector;
152  delete[] fLastDyDx;
153  delete[] fMidVector;
154  delete[] fMidError;
155 
156  delete fAuxStepper;
157 }
ak2
static const G4double ak2
Definition: G4MuPairProductionModel.cc:94
p
const char * p
Definition: xmltok.h:285

Member Function Documentation

G4double G4BogackiShampine45::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 314 of file G4BogackiShampine45.cc.

315 {
316  G4double distLine, distChord;
317  G4ThreeVector initialPoint, finalPoint, midPoint;
318 
319  // Store last initial and final points (they will be overwritten in self-Stepper call!)
320  initialPoint = G4ThreeVector( fLastInitialVector[0],
321  fLastInitialVector[1], fLastInitialVector[2]);
322  finalPoint = G4ThreeVector( fLastFinalVector[0],
323  fLastFinalVector[1], fLastFinalVector[2]);
324 
325 #if 1
326  // Old method -- Do half a step using StepNoErr
327  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
328  fMidVector, fMidError);
329 #else
330  // New method -- Using interpolation, requires only 3 extra stages (ie 3 extra field evaluations )
331 
332  // Use Interpolation, instead of auxiliary stepper to evaluate midpoint
333  if( ! fPreparedInterpolation ) {
334  G4BogackiShampine45 *cThis= const_cast<G4BogackiShampine45 *>(this);
335  cThis-> SetupInterpolationHigh(); // ( fLastInitialVector, fLastDyDx, fLastStepLength );
336  }
337  //For calculating the output at the tau fraction of Step
338  G4double tau = 0.5;
339  // cThis->InterpolateHigh( /* fLastInitialVector, fLastDyDx, fLastStepLength, */, fMidVector, tau);
340  // Old arguments: ( /*yInput, dydx, step,*/ yOut, tau );
341  this->InterpolateHigh( tau, fMidVector );
342 #endif
343 
344  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
345 
346  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
347  // distance of Chord
348 
349  if (initialPoint != finalPoint)
350  {
351  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
352  distChord = distLine;
353  }
354  else
355  {
356  distChord = (midPoint-initialPoint).mag();
357  }
358  return distChord;
359 }
G4BogackiShampine45::InterpolateHigh
void InterpolateHigh(G4double tau, G4double yOut[]) const
Definition: G4BogackiShampine45.cc:551
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
G4BogackiShampine45::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
Definition: G4BogackiShampine45.cc:178
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4BogackiShampine45::GetLastDydx ( G4double  dyDxLast[])

Definition at line 164 of file G4BogackiShampine45.cc.

165 {
166  const G4int numberOfVariables= this->GetNumberOfVariables();
167 
168  for(G4int i=0; i < numberOfVariables; i++ ){
169  dyDxLast[i] = ak9[i];
170  }
171 }
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const

Here is the call graph for this function:

G4int G4BogackiShampine45::IntegratorOrder ( ) const
inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 99 of file G4BogackiShampine45.hh.

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

Definition at line 78 of file G4BogackiShampine45.hh.

83  {
84  InterpolateHigh( tau, yOut);
85  // InterpolateHigh( yInput, dydx, Step, yOut, tau);
86  }
G4BogackiShampine45::InterpolateHigh
void InterpolateHigh(G4double tau, G4double yOut[]) const
Definition: G4BogackiShampine45.cc:551

Here is the call graph for this function:

void G4BogackiShampine45::InterpolateHigh ( G4double  tau,
G4double  yOut[] 
) const

Definition at line 551 of file G4BogackiShampine45.cc.

553 {
554  G4int numberOfVariables = this->GetNumberOfVariables();
555  assert( fPreparedConstants);
556 
557  G4Exception("G4BogackiShampine45::InterpolateHigh()", "GeomField0001",
558  FatalException, "Method is not yet validated.");
559 
560  // const G4double *yIn= fLastInitialVector;
561  // const G4double *dydx= fLastDyDx;
562  const G4double Step = fLastStepLength;
563 
564  // for(G4int i = 0; i< numberOfVariables; i++) yIn[i] = yInput[i];
565 #if 1
566  G4int nwant = numberOfVariables;
567  const G4int norder= 6;
568  G4int l, k;
569 
570  for (l = 0; l < nwant; l++) {
571  yOut[l] = p[norder-1][l] * tau;
572  }
573  for (k = norder - 2; k >= 1; k--) {
574  for (l = 0; l < nwant; l++) {
575  yOut[l] = ( yOut[l] + p[k][l] ) * tau;
576  }
577  }
578  for (l = 0; l < nwant; l++) {
579  yOut[l] = ( yOut[l] + Step * ak8[l] ) * tau + yIn[l];
580  }
581  // The derivative at the end-point is nextDydx[i] = ak8[i];
582 #else
583  // The scheme tries to do the same as the DormandPrince745 routine, but fails
584  G4double b[12];
585  const G4double *dydx= fLastDyDx;
586 
587  G4double tau0 = tau;
588 
589  for(int iStage=1; iStage<=11; iStage++){ // iStage = stage number
590  b[iStage] = 0.0;
591  tau = tau0;
592  for(int j=6; j>=1; j--){ // j reversed
593  b[iStage] += bi[iStage][j] * tau;
594  tau *= tau0;
595  }
596  }
597 
598  for(int i=0; i<numberOfVariables; i++){
599  yOut[i] = yIn[i] + Step*(b[1] * dydx[i] + b[2] * ak2[i] + b[3] * ak3[i] +
600  b[4] * ak4[i] + b[5] * ak5[i] + b[6] * ak6[i] +
601  b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i] +
602  b[10] * ak10[i] + b[11] * ak11[i] );
603  }
604 #endif
605 }
ak2
static const G4double ak2
Definition: G4MuPairProductionModel.cc:94
p
const char * p
Definition: xmltok.h:285
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
test.b
tuple b
Definition: test.py:12
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
FatalException
Definition: G4ExceptionSeverity.hh:60
Step
Definition: Step.hh:41
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4BogackiShampine45::PrepareConstants ( )

Definition at line 480 of file G4BogackiShampine45.cc.

481 {
482  for(int i=1; i<= 11; i++)
483  bi[i][1] = 0.0 ;
484 
485  for(int i=1; i<=6; i++)
486  bi[2][i] = 0.0 ;
487 
488  bi[1][6] = -12134338393.0 / 1050809760.0 ,
489  bi[1][5] = -1620741229.0 / 50038560.0 ,
490  bi[1][4] = -2048058893.0 / 59875200.0 ,
491  bi[1][3] = -87098480009.0 / 5254048800.0 ,
492  bi[1][2] = -11513270273.0 / 3502699200.0 ,
493  //
494  bi[3][6] = -33197340367.0 / 1218433216.0 ,
495  bi[3][5] = -539868024987.0 / 6092166080.0 ,
496  bi[3][4] = -39991188681.0 / 374902528.0 ,
497  bi[3][3] = -69509738227.0 / 1218433216.0 ,
498  bi[3][2] = -29327744613.0 / 2436866432.0 ,
499  //
500  bi[4][6] = -284800997201.0 / 19905339168.0 ,
501  bi[4][5] = -7896875450471.0 / 165877826400.0 ,
502  bi[4][4] = -333945812879.0 / 5671036800.0 ,
503  bi[4][3] = -16209923456237.0 / 497633479200.0 ,
504  bi[4][2] = -2382590741699.0 / 331755652800.0 ,
505  //
506  bi[5][6] = -540919.0 / 741312.0 ,
507  bi[5][5] = -103626067.0 / 43243200.0 ,
508  bi[5][4] = -633779.0 / 211200.0 ,
509  bi[5][3] = -32406787.0 / 18532800.0 ,
510  bi[5][2] = -36591193.0 / 86486400.0 ,
511  //
512  bi[6][6] = 7157998304.0 / 374350977.0 ,
513  bi[6][5] = 30405842464.0 / 623918295.0 ,
514  bi[6][4] = 183022264.0 / 5332635.0 ,
515  bi[6][3] = -3357024032.0 / 1871754885.0 ,
516  bi[6][2] = -611586736.0 / 89131185.0 ,
517  //
518  bi[7][6] = -138073.0 / 9408.0 ,
519  bi[7][5] = -719433.0 / 15680.0 ,
520  bi[7][4] = -1620541.0 / 31360.0 ,
521  bi[7][3] = -385151.0 / 15680.0 ,
522  bi[7][2] = -65403.0 / 15680.0 ,
523  //
524  bi[8][6] = 1245.0 / 64.0 ,
525  bi[8][5] = 3991.0 / 64.0 ,
526  bi[8][4] = 4715.0 / 64.0 ,
527  bi[8][3] = 2501.0 / 64.0 ,
528  bi[8][2] = 149.0 / 16.0 ,
529  bi[8][1] = 1.0 ,
530  //
531  bi[9][6] = 55.0 / 3.0 ,
532  bi[9][5] = 71.0 ,
533  bi[9][4] = 103.0 ,
534  bi[9][3] = 199.0 / 3.0 ,
535  bi[9][2] = 16.0 ,
536  //
537  bi[10][6] = -1774004627.0 / 75810735.0 ,
538  bi[10][5] = -1774004627.0 / 25270245.0 ,
539  bi[10][4] = -26477681.0 / 359975.0 ,
540  bi[10][3] = -11411880511.0 / 379053675.0 ,
541  bi[10][2] = -423642896.0 / 126351225.0 ,
542  //
543  bi[11][6] = 35.0 ,
544  bi[11][5] = 105.0 ,
545  bi[11][4] = 117.0 ,
546  bi[11][3] = 59.0 ,
547  bi[11][2] = 12.0 ;
548  fPreparedConstants= true;
549 }

Here is the caller graph for this function:

void G4BogackiShampine45::SetupInterpolation ( )
inline

Definition at line 71 of file G4BogackiShampine45.hh.

73  {
74  SetupInterpolationHigh(); // ( yInput, dydx, Step);
75  }

Here is the call graph for this function:

void G4BogackiShampine45::SetupInterpolationHigh ( )

Definition at line 361 of file G4BogackiShampine45.cc.

363 {
364  //Coefficients for the additional stages :
365  const G4double
366  a91 = 455.0/6144.0 ,
367  a92 = 0.0 ,
368  a93 = 10256301.0/35409920.0 ,
369  a94 = 2307361.0/17971200.0 ,
370  a95 = -387.0/102400.0 ,
371  a96 = 73.0/5130.0 ,
372  a97 = -7267.0/215040.0 ,
373  a98 = 1.0/32.0 ,
374 
375  a101 = -837888343715.0/13176988637184.0 ,
376  a102 = 30409415.0/52955362.0 ,
377  a103 = -48321525963.0/759168069632.0 ,
378  a104 = 8530738453321.0/197654829557760.0 ,
379  a105 = 1361640523001.0/1626788720640.0 ,
380  a106 = -13143060689.0/38604458898.0 ,
381  a107 = 18700221969.0/379584034816.0 ,
382  a108 = -5831595.0/847285792.0 ,
383  a109 = -5183640.0/26477681.0 ,
384 
385  a111 = 98719073263.0/1551965184000.0 ,
386  a112 = 1307.0/123552.0 ,
387  a113 = 4632066559387.0/70181753241600.0 ,
388  a114 = 7828594302389.0/382182512025600.0 ,
389  a115 = 40763687.0/11070259200.0 ,
390  a116 = 34872732407.0/224610586200.0 ,
391  a117 = -2561897.0/30105600.0 ,
392  a118 = 1.0/10.0 ,
393  a119 = -1.0/10.0 ,
394  a1110 = -1403317093.0/11371610250.0 ;
395 
396  const G4int numberOfVariables= this->GetNumberOfVariables();
397 
398  // const G4double *yIn= fLastInitialVector;
399  const G4double *dydx= fLastDyDx;
400  const G4double Step = fLastStepLength;
401 
402  yTemp[7] = yIn[7];
403 
404  //Evaluate the extra stages :
405  for(int i=0; i<numberOfVariables; i++){
406  yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
407  a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
408  a97*ak7[i] + a98*ak8[i] );
409  }
410 
411  RightHandSide(yTemp, ak9); //9th stage
412 
413  for(int i=0; i<numberOfVariables; i++){
414  yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
415  a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
416  a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
417  }
418 
419  RightHandSide(yTemp, ak10); //10th stage
420 
421  for(int i=0; i<numberOfVariables; i++){
422  yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
423  a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
424  a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
425  a1110*ak10[i] );
426  }
427  RightHandSide(yTemp, ak11); //11th stage
428 
429  // In future we can restrict the number of variables interpolated
430  int nwant = numberOfVariables;
431 
432 // Form the coefficients of the interpolating polynomial in its shifted
433 // and scaled form. The terms are grouped to minimize the errors
434 // of the transformation, to cope with ill-conditioning. ( From RKSUITE )
435 //
436  for (int l = 0; l < nwant; l++) {
437  // Coefficient of tau^6
438  p[5][l] = bi[5][6]*ak5[l] +
439  ((bi[10][6]*ak10[l] + bi[8][6]*ak8[l]) +
440  (bi[7][6]*ak7[l] + bi[6][6]*ak6[l])) +
441  ((bi[4][6]*ak4[l] + bi[9][6]*ak9[l]) +
442  (bi[3][6]*ak3[l] + bi[11][6]*ak11[l]) +
443  bi[1][6]*dydx[l]);
444  // Coefficient of tau^5
445  p[4][l] = (bi[10][5]*ak10[l] + bi[9][5]*ak9[l]) +
446  ((bi[7][5]*ak7[l] + bi[6][5]*ak6[l]) +
447  bi[5][5]*ak5[l]) + ((bi[4][5]*ak4[l] +
448  bi[8][5]*ak8[l]) + (bi[3][5]*ak3[l] +
449  bi[11][5]*ak11[l]) + bi[1][5]*dydx[l]);
450  // Coefficient of tau^4
451  p[3][l] = ((bi[4][4]*ak4[l] + bi[8][4]*ak8[l]) +
452  (bi[7][4]*ak7[l] + bi[6][4]*ak6[l]) +
453  bi[5][4]*ak5[l]) + ((bi[10][4]*ak10[l] +
454  bi[9][4]*ak9[l]) + (bi[3][4]*ak3[l] +
455  bi[11][4]*ak11[l]) + bi[1][4]*dydx[l]);
456  // Coefficient of tau^3
457  p[2][l] = bi[5][3]*ak5[l] + bi[6][3]*ak6[l] +
458  ((bi[3][3]*ak3[l] + bi[9][3]*ak9[l]) +
459  (bi[10][3]*ak10[l]+ bi[8][3]*ak8[l]) + bi[1][3]*dydx[l]) +
460  ((bi[4][3]*ak4[l] + bi[11][3]*ak11[l]) + bi[7][3]*ak7[l]);
461  // Coefficient of tau^2
462  p[1][l] = bi[5][2]*ak5[l] + ((bi[6][2]*ak6[l] +
463  bi[8][2]*ak8[l]) + bi[1][2]*dydx[l]) +
464  ((bi[3][2]*ak3[l] + bi[9][2]*ak9[l]) +
465  bi[10][2]*ak10[l])+ ((bi[4][2]*ak4[l] +
466  bi[11][2]*ak2[l]) + bi[7][2]*ak7[l]);
467  }
468 //
469 // Scale all the coefficients by the step size.
470 //
471  for (int i = 0; i < 6; i++) {
472  for (int l = 0; l < nwant; l++) {
473  p[i][l] *= Step;
474  }
475  }
476 
477  fPreparedInterpolation= true;
478 }
ak2
static const G4double ak2
Definition: G4MuPairProductionModel.cc:94
p
const char * p
Definition: xmltok.h:285
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
Step
Definition: Step.hh:41
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:

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

Implements G4MagIntegratorStepper.

Definition at line 178 of file G4BogackiShampine45.cc.

183 {
184  G4int i;
185 
186  // Constants from the Butcher tableu
187  const G4double
188  b21 = 1.0/6.0 ,
189  b31 = 2.0/27.0 , b32 = 4.0/27.0,
190 
191  b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
192 
193  b51 = 68.0/297.0, b52 = -4.0/11.0,
194  b53 = 42.0/143.0, b54 = 1960.0/3861.0,
195 
196  b61 = 597.0/22528.0, b62 = 81.0/352.0,
197  b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
198  b65 = 4617.0/20480.0,
199 
200  b71 = 174197.0/959244.0, b72 = -30942.0/79937.0,
201  b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
202  b75 = -29421.0/29068.0, b76 = 482048.0/414219.0,
203 
204  b81 = 587.0/8064.0, b82 = 0.0,
205  b83 = 4440339.0/15491840.0, b84 = 24353.0/124800.0,
206  b85 = 387.0/44800.0, b86 = 2152.0/5985.0,
207  b87 = 7267.0/94080.0;
208 
209 // c1 = 2479.0/34992.0,
210 // c2 = 0.0,
211 // c3 = 123.0/416.0,
212 // c4 = 612941.0/3411720.0,
213 // c5 = 43.0/1440.0,
214 // c6 = 2272.0/6561.0,
215 // c7 = 79937.0/1113912.0,
216 // c8 = 3293.0/556956.0,
217 
218  // For the embedded higher order method only the difference of values
219  // taken and is used directly later (instead of defining the last row
220  // of Butcher table in separate constants and taking the
221  // difference)
222  const G4double
223  dc1 = b81 - 2479.0 / 34992.0 ,
224  dc2 = 0.0,
225  dc3 = b83 - 123.0 / 416.0 ,
226  dc4 = b84 - 612941.0 / 3411720.0,
227  dc5 = b85 - 43.0 / 1440.0,
228  dc6 = b86 - 2272.0 / 6561.0,
229  dc7 = b87 - 79937.0 / 1113912.0,
230  dc8 = - 3293.0 / 556956.0;
231 
232  const G4int numberOfVariables= this->GetNumberOfVariables();
233 
234  // The number of variables to be integrated over
235  yOut[7] = yTemp[7] = yIn[7];
236  // Saving yInput because yInput and yOut can be aliases for same array
237 
238  for(i=0;i<numberOfVariables;i++)
239  {
240  yIn[i]=yInput[i];
241  }
242 
243  // RightHandSide(yIn, dydx) ;
244  // 1st Step - Not doing, getting passed
245 
246  for(i=0;i<numberOfVariables;i++)
247  {
248  yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
249  }
250  RightHandSide(yTemp, ak2) ; // 2nd Step
251 
252  for(i=0;i<numberOfVariables;i++)
253  {
254  yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
255  }
256  RightHandSide(yTemp, ak3) ; // 3rd Step
257 
258  for(i=0;i<numberOfVariables;i++)
259  {
260  yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
261  }
262  RightHandSide(yTemp, ak4) ; // 4th Step
263 
264  for(i=0;i<numberOfVariables;i++)
265  {
266  yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
267  b54*ak4[i]) ;
268  }
269  RightHandSide(yTemp, ak5) ; // 5th Step
270 
271  for(i=0;i<numberOfVariables;i++)
272  {
273  yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
274  b64*ak4[i] + b65*ak5[i]) ;
275  }
276  RightHandSide(yTemp, ak6) ; // 6th Step
277 
278  for(i=0;i<numberOfVariables;i++)
279  {
280  yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
281  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
282  }
283  RightHandSide(yTemp, ak7); //7th Step
284 
285  for(i=0;i<numberOfVariables;i++)
286  {
287  yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] +
288  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
289  b87*ak7[i]);
290  }
291  RightHandSide(yOut, ak8); //8th Step - Final one Using FSAL
292 
293  for(i=0;i<numberOfVariables;i++)
294  {
295  yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
296  dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
297 
298  // Store Input and Final values, for possible use in calculating chord
299  fLastInitialVector[i] = yIn[i] ;
300  fLastFinalVector[i] = yOut[i];
301  fLastDyDx[i] = DyDx[i];
302  }
303 
304  fLastStepLength = Step;
305  fPreparedInterpolation= false;
306 
307  return ;
308 }
ak2
static const G4double ak2
Definition: G4MuPairProductionModel.cc:94
G4int
int G4int
Definition: G4Types.hh:78
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
Step
Definition: Step.hh:41
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: