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

#include <G4DormandPrinceRK56.hh>

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

Public Member Functions

 G4DormandPrinceRK56 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 
 ~G4DormandPrinceRK56 ()
 
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
 
G4double DistChord () const
 
G4int IntegratorOrder () const
 
 G4DormandPrinceRK56 (const G4DormandPrinceRK56 &)
 
G4DormandPrinceRK56 & operator= (const G4DormandPrinceRK56 &)
 
void SetupInterpolate_low (const G4double yInput[], const G4double dydx[], const G4double Step)
 
void Interpolate_low (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
 
void SetupInterpolation ()
 
void SetupInterpolate (const G4double yInput[], const G4double dydx[], const G4double Step)
 
void Interpolate (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
 
void Interpolate (G4double tau, G4double yOut[])
 
void SetupInterpolate_high (const G4double yInput[], const G4double dydx[], const G4double Step)
 
void Interpolate_high (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
 
- 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 G4DormandPrinceRK56.hh.

Constructor & Destructor Documentation

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

Definition at line 47 of file G4DormandPrinceRK56.cc.

50 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
51  fAuxStepper(0)
52 {
53 
54  const G4int numberOfVariables = noIntegrationVariables;
55 
56  //New Chunk of memory being created for use by the stepper
57 
58  //aki - for storing intermediate RHS
59  ak2 = new G4double[numberOfVariables];
60  ak3 = new G4double[numberOfVariables];
61  ak4 = new G4double[numberOfVariables];
62  ak5 = new G4double[numberOfVariables];
63  ak6 = new G4double[numberOfVariables];
64  ak7 = new G4double[numberOfVariables];
65  ak8 = new G4double[numberOfVariables];
66  ak9 = new G4double[numberOfVariables];
67 
68  const G4int numStateVars = std::max(noIntegrationVariables, 8);
69  yTemp = new G4double[numStateVars];
70  yIn = new G4double[numStateVars] ;
71 
72  fLastInitialVector = new G4double[numStateVars] ;
73  fLastFinalVector = new G4double[numStateVars] ;
74 
75  fLastDyDx = new G4double[numStateVars];
76 
77  fMidVector = new G4double[numStateVars];
78  fMidError = new G4double[numStateVars];
79 
80  if( primary )
81  {
82  fAuxStepper = new G4DormandPrinceRK56(EqRhs, numberOfVariables,
83  !primary);
84  }
85 }
G4DormandPrinceRK56::G4DormandPrinceRK56
G4DormandPrinceRK56(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition: G4DormandPrinceRK56.cc:47
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
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4DormandPrinceRK56::~G4DormandPrinceRK56 ( )

Definition at line 89 of file G4DormandPrinceRK56.cc.

89  {
90  //clear all previously allocated memory for stepper and DistChord
91  delete[] ak2;
92  delete[] ak3;
93  delete[] ak4;
94  delete[] ak5;
95  delete[] ak6;
96  delete[] ak7;
97  delete[] ak8;
98  delete[] ak9;
99 
100  delete[] yTemp;
101  delete[] yIn;
102 
103  delete[] fLastInitialVector;
104  delete[] fLastFinalVector;
105  delete[] fLastDyDx;
106  delete[] fMidVector;
107  delete[] fMidError;
108 
109  delete fAuxStepper;
110 
111 }
G4DormandPrinceRK56::G4DormandPrinceRK56 ( const G4DormandPrinceRK56 &  )

Member Function Documentation

G4double G4DormandPrinceRK56::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 372 of file G4DormandPrinceRK56.cc.

373 {
374  G4double distLine, distChord;
375  G4ThreeVector initialPoint, finalPoint, midPoint;
376 
377  // Store last initial and final points (they will be overwritten in self-Stepper call!)
378  initialPoint = G4ThreeVector( fLastInitialVector[0],
379  fLastInitialVector[1], fLastInitialVector[2]);
380  finalPoint = G4ThreeVector( fLastFinalVector[0],
381  fLastFinalVector[1], fLastFinalVector[2]);
382 
383  // Do half a step using StepNoErr
384 
385  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
386  fMidVector, fMidError );
387 
388  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
389 
390  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
391  // distance of Chord
392 
393 
394  if (initialPoint != finalPoint)
395  {
396  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
397  distChord = distLine;
398  }
399  else
400  {
401  distChord = (midPoint-initialPoint).mag();
402  }
403  return distChord;
404 }
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
G4DormandPrinceRK56::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
Definition: G4DormandPrinceRK56.cc:119
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4int G4DormandPrinceRK56::IntegratorOrder ( ) const
inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 71 of file G4DormandPrinceRK56.hh.

71 { return 5; }
void G4DormandPrinceRK56::Interpolate ( const G4double  yInput[],
const G4double  dydx[],
const G4double  Step,
G4double  yOut[],
G4double  tau 
)
inline

Definition at line 98 of file G4DormandPrinceRK56.hh.

102  {
103  Interpolate_low( yInput, dydx, Step, yOut, tau);
104  }
G4DormandPrinceRK56::Interpolate_low
void Interpolate_low(const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
Definition: G4DormandPrinceRK56.cc:451

Here is the call graph for this function:

Here is the caller graph for this function:

void G4DormandPrinceRK56::Interpolate ( G4double  tau,
G4double  yOut[] 
)
inline

Definition at line 106 of file G4DormandPrinceRK56.hh.

107  { Interpolate( fLastInitialVector, fLastDyDx, fLastStepLength, yOut, tau ); }
G4DormandPrinceRK56::Interpolate
void Interpolate(const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
Definition: G4DormandPrinceRK56.hh:98

Here is the call graph for this function:

void G4DormandPrinceRK56::Interpolate_high ( const G4double  yInput[],
const G4double  dydx[],
const G4double  Step,
G4double  yOut[],
G4double  tau 
)

Definition at line 604 of file G4DormandPrinceRK56.cc.

609 {
610 
611  //Define the coefficients for the polynomials
612  G4double bi[13][6], b[13];
613  G4int numberOfVariables = this->GetNumberOfVariables();
614 
615 
616  // COEFFICIENTS OF bi[ 1]
617  bi[1][0] = 1.0 ,
618  bi[1][1] = -18487.0/2880.0 ,
619  bi[1][2] = 139189.0/7200.0 ,
620  bi[1][3] = -53923.0/1800.0 ,
621  bi[1][4] = 13811.0/600,
622  bi[1][5] = -2071.0/300,
623  // --------------------------------------------------------
624  //
625  // COEFFICIENTS OF bi[2]
626  bi[2][0] = 0.0 ,
627  bi[2][1] = 0.0 ,
628  bi[2][2] = 0.0 ,
629  bi[2][3] = 0.0 ,
630  bi[2][4] = 0.0 ,
631  bi[2][5] = 0.0 ,
632  // --------------------------------------------------------
633  //
634  // COEFFICIENTS OF bi[3]
635  bi[3][0] = 0.0 ,
636  bi[3][1] = 0.0 ,
637  bi[3][2] = 0.0 ,
638  bi[3][3] = 0.0 ,
639  bi[3][4] = 0.0 ,
640  bi[3][5] = 0.0 ,
641  // --------------------------------------------------------
642  //
643  // COEFFICIENTS OF bi[4]
644  bi[4][0] = 0.0 ,
645  bi[4][1] = -30208.0/14157.0 ,
646  bi[4][2] = 1147904.0/70785.0 ,
647  bi[4][3] = -241664.0/5445.0 ,
648  bi[4][4] = 241664.0/4719.0 ,
649  bi[4][5] = -483328.0/23595.0 ,
650  // --------------------------------------------------------
651  //
652  // COEFFICIENTS OF bi[5]
653  bi[5][0] = 0.0 ,
654  bi[5][1] = -177147.0/32912.0 ,
655  bi[5][2] = 3365793.0/82280.0 ,
656  bi[5][3] = -2302911.0/20570.0 ,
657  bi[5][4] = 531441.0/4114.0 ,
658  bi[5][5] = -531441.0/10285.0 ,
659  // --------------------------------------------------------
660  //
661  // COEFFICIENTS OF bi[6]
662  bi[6][0] = 0.0 ,
663  bi[6][1] = 536.0/141.0 ,
664  bi[6][2] = -20368.0/705.0 ,
665  bi[6][3] = 55744.0/705.0 ,
666  bi[6][4] = -4288.0/47.0 ,
667  bi[6][5] = 8576.0/235,
668  // --------------------------------------------------------
669  //
670  // COEFFICIENTS OF bi[7]
671  bi[7][0] = 0.0 ,
672  bi[7][1] = -1977326743.0/723932352.0 ,
673  bi[7][2] = 37569208117.0/1809830880.0 ,
674  bi[7][3] = -1977326743.0/34804440.0 ,
675  bi[7][4] = 1977326743.0/30163848.0 ,
676  bi[7][5] = -1977326743.0/75409620.0 ,
677  // --------------------------------------------------------
678  //
679  // COEFFICIENTS OF bi[8]
680  bi[8][0] = 0.0 ,
681  bi[8][1] = 259.0/144.0 ,
682  bi[8][2] = -4921.0/360.0 ,
683  bi[8][3] = 3367.0/90.0 ,
684  bi[8][4] = -259.0/6.0 ,
685  bi[8][5] = 259.0/15.0 ,
686  // --------------------------------------------------------
687  //
688  // COEFFICIENTS OF bi[9]
689  bi[9][0] = 0.0 ,
690  bi[9][1] = 62.0/105.0 ,
691  bi[9][2] = -2381.0/525.0 ,
692  bi[9][3] = 949.0/75.0 ,
693  bi[9][4] = -2636.0/175.0 ,
694  bi[9][5] = 1112.0/175.0 ,
695  // --------------------------------------------------------
696  //
697  // COEFFICIENTS OF bi[10]
698  bi[10][0] = 0.0 ,
699  bi[10][1] = 43.0/3.0 ,
700  bi[10][2] = -1534.0/15.0 ,
701  bi[10][3] = 3767.0/15.0 ,
702  bi[10][4] = -1264.0/5.0 ,
703  bi[10][5] = 448.0/5.0 ,
704  // --------------------------------------------------------
705  //
706  // COEFFICIENTS OF bi[11]
707  bi[11][0] = 0.0 ,
708  bi[11][1] = 63.0/5.0 ,
709  bi[11][2] = -1494.0/25.0 ,
710  bi[11][3] = 2907.0/25.0 ,
711  bi[11][4] = -2592.0/25.0 ,
712  bi[11][5] = 864.0/25.0 ,
713  // --------------------------------------------------------
714  //
715  // COEFFICIENTS OF bi[12]
716  bi[12][0] = 0.0 ,
717  bi[12][1] = -576.0/35.0 ,
718  bi[12][2] = 19584.0/175.0 ,
719  bi[12][3] = -6336.0/25.0 ,
720  bi[12][4] = 41472.0/175.0 ,
721  bi[12][5] = -13824.0/175.0 ;
722  // --------------------------------------------------------
723 
724 
725  for(G4int i = 0; i< numberOfVariables; i++)
726  yIn[i] = yInput[i];
727 
728  G4double tau0 = tau;
729  // Calculating the polynomials (coefficents for the respective stages) :
730 
731  for(int i=1; i<=12; i++){ //Here i is NOT the coordinate no. , it's stage no.
732  b[i] = 0;
733  tau = 1.0;
734  for(int j=0; j<=5; j++){
735  b[i] += bi[i][j]*tau;
736  tau*=tau0;
737  }
738  }
739 
740 // Calculating the interpolation at the fraction tau of the step using the polynomial
741 // coefficients and the respective stages
742 
743  for(int i=0; i<numberOfVariables; i++){ //Here i IS the cooridnate no.
744  yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
745  b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
746  b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
747  b[10]*ak10[i] + b[11]*ak11[i] + b[12]*ak12[i]);
748  }
749 }
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 G4DormandPrinceRK56::Interpolate_low ( const G4double  yInput[],
const G4double  dydx[],
const G4double  Step,
G4double  yOut[],
G4double  tau 
)

Definition at line 451 of file G4DormandPrinceRK56.cc.

455  {
456  {
457 
458  G4double
459  bf1, bf4, bf5, bf6, bf7, bf8, bf9, bf10;
460 
461  G4double tau0 = tau;
462  const G4int numberOfVariables= this->GetNumberOfVariables();
463 
464  for(int i=0;i<numberOfVariables;i++)
465  {
466  yIn[i]=yInput[i];
467  }
468 
469  G4double
470  tau_2 = tau0*tau0 ,
471  tau_3 = tau0*tau_2,
472  tau_4 = tau_2*tau_2;
473 
474  //bf2 = bf3 = 0
475  bf1 = (66480.0*tau_4 - 206243.0*tau_3 + 237786.0*tau_2 - 124793.0*tau + 28800.0)/28800.0 ,
476  bf4 = -16.0*tau*(45312.0*tau_3 - 125933.0*tau_2 + 119706.0*tau -40973.0)/70785.0 ,
477  bf5 = -2187.0*tau*(19440.0*tau_3 - 45743.0*tau_2 + 34786.0*tau - 9293.0)/1645600.0 ,
478  bf6 = tau*(12864.0*tau_3 - 30653.0*tau_2 + 23786.0*tau - 6533.0)/705.0 ,
479  bf7 = -5764801.0*tau*(16464.0*tau_3 - 32797.0*tau_2 + 17574.0*tau - 1927.0)/7239323520.0 ,
480  bf8 = 37.0*tau*(336.0*tau_3 - 661.0*tau_2 + 342.0*tau -31.0)/1440.0 ,
481  bf9 = tau*(tau-1.0)*(16.0*tau_2 - 15.0*tau +3.0)/4.0 ,
482  bf10 = 8.0*tau*(tau - 1.0)*(tau - 1.0)*(2.0*tau - 1.0) ;
483 
484  for( int i=0; i<numberOfVariables; i++){
485  yOut[i] = yIn[i] + Step*tau*( bf1*dydx[i] + bf4*ak4[i] + bf5*ak5[i] +
486  bf6*ak6[i] + bf7*ak7[i] + bf8*ak8[i] +
487  bf9*ak9[i] + bf10*ak10_low[i] ) ;
488  }
489 
490 
491 
492  }
493 
494 }
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:

G4DormandPrinceRK56& G4DormandPrinceRK56::operator= ( const G4DormandPrinceRK56 &  )
void G4DormandPrinceRK56::SetupInterpolate ( const G4double  yInput[],
const G4double  dydx[],
const G4double  Step 
)
inline

Definition at line 91 of file G4DormandPrinceRK56.hh.

93  {
94  SetupInterpolate_low( yInput, dydx, Step);
95  }
G4DormandPrinceRK56::SetupInterpolate_low
void SetupInterpolate_low(const G4double yInput[], const G4double dydx[], const G4double Step)
Definition: G4DormandPrinceRK56.cc:417

Here is the call graph for this function:

Here is the caller graph for this function:

void G4DormandPrinceRK56::SetupInterpolate_high ( const G4double  yInput[],
const G4double  dydx[],
const G4double  Step 
)

Definition at line 510 of file G4DormandPrinceRK56.cc.

512  {
513 
514  //Coefficients for the additional stages :
515 
516  G4double
517  b101 = 33797.0/460800.0 ,
518  b102 = 0.0 ,
519  b103 = 0.0 ,
520  b104 = 27757.0/70785.0 ,
521  b105 = 7923501.0/26329600.0 ,
522  b106 = -927.0/3760.0 ,
523  b107 = -3314760575.0/23165835264.0 ,
524  b108 = 2479.0/23040.0 ,
525  b109 = 1.0/64.0 ,
526 
527  b111 = 5843.0/76800.0 ,
528  b112 = 0.0 ,
529  b113 = 0.0 ,
530  b114 = 464.0/2673.0 ,
531  b115 = 353997.0/1196800.0 ,
532  b116 = -15068.0/57105.0 ,
533  b117 = -282475249.0/3644974080.0 ,
534  b118 = 8678831.0/156245760.0 ,
535  b119 = 116113.0/11718432.0 ,
536  b1110 = -25.0/243.0 ,
537 
538  b121 = 15088049.0/199065600.0 ,
539  b122 = 0.0 ,
540  b123 = 0.0 ,
541  b124 = 2.0/5.0 ,
542  b125 = 92222037.0/268083200.0 ,
543  b126 = -433420501.0/1528586640.0 ,
544  b127 = -11549242677007.0/83630285291520.0 ,
545  b128 = 2725085557.0/26167173120.0 ,
546  b129 = 235429367.0/16354483200.0 ,
547  b1210 = -90924917.0/1040739840.0 ,
548  b1211 = -271149.0/21414400.0 ;
549 
550 
551  const G4int numberOfVariables= this->GetNumberOfVariables();
552 
553  // Saving yInput because yInput and yOut can be aliases for same array
554  for(int i=0;i<numberOfVariables;i++)
555  {
556  yIn[i]=yInput[i];
557  }
558 
559  yTemp[7] = yIn[7];
560 
561 
562  // New memory for Additional stages
563 
564  if(ak10 == NULL)
565  ak10 = new G4double[numberOfVariables];
566  if(ak11 == NULL)
567  ak11 = new G4double[numberOfVariables];
568  if(ak12 == NULL)
569  ak12 = new G4double[numberOfVariables];
570 
571  //Evaluate the extra stages :
572 
573  for(int i=0;i<numberOfVariables;i++)
574  {
575  yTemp[i] = yIn[i] + Step*(b101*dydx[i] + b102*ak2[i] + b103*ak3[i] +
576  b104*ak4[i] + b105*ak5[i] + b106*ak6[i] +
577  b107*ak7[i] + b108*ak8[i] + b109*ak9[i]);
578  }
579  RightHandSide(yTemp, ak10); //10th Stage
580 
581  for(int i=0;i<numberOfVariables;i++)
582  {
583  yTemp[i] = yIn[i] + Step*(b111*dydx[i] + b112*ak2[i] + b113*ak3[i] +
584  b114*ak4[i] + b115*ak5[i] + b116*ak6[i] +
585  b117*ak7[i] + b118*ak8[i] + b119*ak9[i] +
586  b1110*ak10[i]);
587  }
588  RightHandSide(yTemp, ak11); //11th Stage
589 
590  for(int i=0;i<numberOfVariables;i++)
591  {
592  yTemp[i] = yIn[i] + Step*(b121*dydx[i] + b122*ak2[i] + b123*ak3[i] +
593  b124*ak4[i] + b125*ak5[i] + b126*ak6[i] +
594  b127*ak7[i] + b128*ak8[i] + b129*ak9[i] +
595  b1210*ak10[i] + b1211*ak11[i]);
596  }
597  RightHandSide(yTemp, ak12); //12th Stage
598 
599 }
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 G4DormandPrinceRK56::SetupInterpolate_low ( const G4double  yInput[],
const G4double  dydx[],
const G4double  Step 
)

Definition at line 417 of file G4DormandPrinceRK56.cc.

419  {
420  const G4int numberOfVariables= this->GetNumberOfVariables();
421 
422  G4double
423  b_101 = 33797.0/460800.0 ,
424  b_102 = 0. ,
425  b_103 = 0. ,
426  b_104 = 27757.0/70785.0 ,
427  b_105 = 7923501.0/26329600.0 ,
428  b_106 = -927.0/3760.0 ,
429  b_107 = -3314760575.0/23165835264.0 ,
430  b_108 = 2479.0/23040.0 ,
431  b_109 = 1.0/64.0 ;
432 
433  ak10_low = new G4double[numberOfVariables];
434 
435  for(int i=0;i<numberOfVariables;i++)
436  {
437  yIn[i]=yInput[i];
438  }
439 
440 
441  for(int i=0;i<numberOfVariables;i++)
442  {
443  yTemp[i] = yIn[i] + Step*(b_101*dydx[i] + b_102*ak2[i] + b_103*ak3[i] +
444  b_104*ak4[i] + b_105*ak5[i] + b_106*ak6[i] +
445  b_107*ak7[i] + b_108*ak8[i] + b_109*ak9[i]);
446  }
447  RightHandSide(yTemp, ak10_low); //10th Stage
448 
449 }
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:

void G4DormandPrinceRK56::SetupInterpolation ( )
inline

Definition at line 88 of file G4DormandPrinceRK56.hh.

89  { SetupInterpolate( fLastInitialVector, fLastDyDx, fLastStepLength); }
G4DormandPrinceRK56::SetupInterpolate
void SetupInterpolate(const G4double yInput[], const G4double dydx[], const G4double Step)
Definition: G4DormandPrinceRK56.hh:91

Here is the call graph for this function:

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

Implements G4MagIntegratorStepper.

Definition at line 119 of file G4DormandPrinceRK56.cc.

126 {
127  G4int i;
128 
129  //The various constants defined on the basis of butcher tableu
130  const G4double //G4double - only once
131 
132  // Old Coefficients from
133  // Table 1. RK6(5)8M
134 //---Ref---
135 //[P. J. Prince and J. R. Dormand, “High order embedded Runge-Kutta formulae,”
136 // Journal of Computational and Applied Mathematics, vol. 7, no. 1, pp. 67–75,
137 // Dec. 1980.
138 //----------------
139 
140  b21 = 1.0/10.0 ,
141 
142  b31 = -2.0/81.0 ,
143  b32 = 20.0/81.0 ,
144 
145  b41 = 615.0/1372.0 ,
146  b42 = -270.0/343.0 ,
147  b43 = 1053.0/1372.0 ,
148 
149  b51 = 3243.0/5500.0 ,
150  b52 = -54.0/55.0 ,
151  b53 = 50949.0/71500.0 ,
152  b54 = 4998.0/17875.0 ,
153 
154  b61 = -26492.0/37125.0 ,
155  b62 = 72.0/55.0 ,
156  b63 = 2808.0/23375.0 ,
157  b64 = -24206.0/37125.0 ,
158  b65 = 338.0/459.0 ,
159 
160  b71 = 5561.0/2376.0 ,
161  b72 = -35.0/11.0 ,
162  b73 = -24117.0/31603.0 ,
163  b74 = 899983.0/200772.0 ,
164  b75 = -5225.0/1836.0 ,
165  b76 = 3925.0/4056.0 ,
166 
167  b81 = 465467.0/266112.0 ,
168  b82 = -2945.0/1232.0 ,
169  b83 = -5610201.0/14158144.0 ,
170  b84 = 10513573.0/3212352.0 ,
171  b85 = -424325.0/205632.0 ,
172  b86 = 376225.0/454272.0 ,
173  b87 = 0.0 ,
174 
175  c1 = 61.0/864.0 ,
176  c2 = 0.0 ,
177  c3 = 98415.0/321776.0 ,
178  c4 = 16807.0/146016.0 ,
179  c5 = 1375.0/7344.0 ,
180  c6 = 1375.0/5408.0 ,
181  c7 = -37.0/1120.0 ,
182  c8 = 1.0/10.0 ,
183 
184  b91 = 61.0/864.0 ,
185  b92 = 0.0 ,
186  b93 = 98415.0/321776.0 ,
187  b94 = 16807.0/146016.0 ,
188  b95 = 1375.0/7344.0 ,
189  b96 = 1375.0/5408.0 ,
190  b97 = -37.0/1120.0 ,
191  b98 = 1.0/10.0 ,
192 
193  dc1 = c1 - 821.0/10800.0 ,
194  dc2 = c2 - 0.0 ,
195  dc3 = c3 - 19683.0/71825,
196  dc4 = c4 - 175273.0/912600.0 ,
197  dc5 = c5 - 395.0/3672.0 ,
198  dc6 = c6 - 785.0/2704.0 ,
199  dc7 = c7 - 3.0/50.0 ,
200  dc8 = c8 - 0.0 ,
201  dc9 = 0.0;
202 
203 
204 // New Coefficients obtained from
205  // Table 3 RK6(5)9FM with corrected coefficients
206 //---Ref---
207 // T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
208 // “Continuous approximation with embedded Runge-Kutta methods,”
209 // Applied Numerical Mathematics, vol. 22, no. 1, pp. 51–62, 1996.
210 //------------------------------------
211 
212 // b21 = 1.0/9.0 ,
213 //
214 // b31 = 1.0/24.0 ,
215 // b32 = 1.0/8.0 ,
216 //
217 // b41 = 1.0/16.0 ,
218 // b42 = 0.0 ,
219 // b43 = 3.0/16.0 ,
220 //
221 // b51 = 280.0/729.0 ,
222 // b52 = 0.0 ,
223 // b53 = -325.0/243.0 ,
224 // b54 = 1100.0/729.0 ,
225 //
226 // b61 = 6127.0/14680.0 ,
227 // b62 = 0.0 ,
228 // b63 = -1077.0/734.0 ,
229 // b64 = 6494.0/4037.0 ,
230 // b65 = -9477.0/161480.0 ,
231 //
232 // b71 = -13426273320.0/14809773769.0 ,
233 // b72 = 0.0 ,
234 // b73 = 4192558704.0/2115681967.0 ,
235 // b74 = 14334750144.0/14809773769.0 ,
236 // b75 = 117092732328.0/14809773769.0 ,
237 // b76 = -361966176.0/40353607.0 ,
238 //
239 // b81 = -2340689.0/1901060.0 ,
240 // b82 = 0.0 ,
241 // b83 = 31647.0/13579.0 ,
242 // b84 = 253549596.0/149518369.0 ,
243 // b85 = 10559024082.0/977620105.0 ,
244 // b86 = -152952.0/12173.0 ,
245 // b87 = -5764801.0/186010396.0 ,
246 //
247 // b91 = 203.0/2880.0 ,
248 // b92 = 0.0 ,
249 // b93 = 0.0 ,
250 // b94 = 30208.0/70785.0 ,
251 // b95 = 177147.0/164560.0 ,
252 // b96 = -536.0/705.0 ,
253 // b97 = 1977326743.0/3619661760.0 ,
254 // b98 = -259.0/720.0 ,
255 //
256 //
257 // dc1 = 36567.0/458800.0 - b91,
258 // dc2 = 0.0 - b92,
259 // dc3 = 0.0 - b93,
260 // dc4 = 9925984.0/27063465.0 - b94,
261 // dc5 = 85382667.0/117968950.0 - b95,
262 // dc6 = - 310378.0/808635.0 - b96 ,
263 // dc7 = 262119736669.0/345979336560.0 - b97,
264 // dc8 = - 1.0/2.0 - b98 ,
265 // dc9 = -101.0/2294.0 ;
266 
267 
268  //end of declaration
269 
270 
271  const G4int numberOfVariables= this->GetNumberOfVariables();
272 
273  // The number of variables to be integrated over
274  yOut[7] = yTemp[7] = yIn[7];
275  // Saving yInput because yInput and yOut can be aliases for same array
276 
277  for(i=0;i<numberOfVariables;i++)
278  {
279  yIn[i]=yInput[i];
280  }
281 
282 
283 
284  // RightHandSide(yIn, dydx) ;
285  // 1st Step - Not doing, getting passed
286 
287  for(i=0;i<numberOfVariables;i++)
288  {
289  yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
290  }
291  RightHandSide(yTemp, ak2) ; // 2nd Stage
292 
293  for(i=0;i<numberOfVariables;i++)
294  {
295  yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
296  }
297  RightHandSide(yTemp, ak3) ; // 3rd Stage
298 
299  for(i=0;i<numberOfVariables;i++)
300  {
301  yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
302  }
303  RightHandSide(yTemp, ak4) ; // 4th Stage
304 
305  for(i=0;i<numberOfVariables;i++)
306  {
307  yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
308  b54*ak4[i]) ;
309  }
310  RightHandSide(yTemp, ak5) ; // 5th Stage
311 
312  for(i=0;i<numberOfVariables;i++)
313  {
314  yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
315  b64*ak4[i] + b65*ak5[i]) ;
316  }
317  RightHandSide(yTemp, ak6) ; // 6th Stage
318 
319  for(i=0;i<numberOfVariables;i++)
320  {
321  yTemp[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
322  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
323  }
324  RightHandSide(yTemp, ak7); //7th Stage
325 
326  for(i=0;i<numberOfVariables;i++)
327  {
328  yTemp[i] = yIn[i] + Step*(b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
329  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
330  b87*ak7[i]);
331  }
332  RightHandSide(yTemp, ak8); //8th Stage
333 
334  for(i=0;i<numberOfVariables;i++)
335  {
336  yOut[i] = yIn[i] + Step*(b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
337  b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
338  b97*ak7[i] + b98*ak8[i] );
339  }
340  RightHandSide(yOut, ak9); //9th Stage
341 
342 
343  for(i=0;i<numberOfVariables;i++)
344  {
345 
346  // Estimate error as difference between 5th and
347  // 6th order methods
348 
349  yErr[i] = Step*( dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i]
350  + dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]
351  + dc9*ak9[i] ) ;
352 
353  // - Saving 'estimated' derivative at end-point
354  // nextDydx[i] = ak9[i];
355 
356  // Store Input and Final values, for possible use in calculating chord
357  fLastInitialVector[i] = yIn[i] ;
358  fLastFinalVector[i] = yOut[i];
359  fLastDyDx[i] = dydx[i];
360 
361  }
362 
363  fLastStepLength = Step;
364 
365  return ;
366 }
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: