G4BogackiShampine45 Class Reference

#include <G4BogackiShampine45.hh>

Inheritance diagram for G4BogackiShampine45:

Public Member Functions
	G4BogackiShampine45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
	~G4BogackiShampine45 () override
	G4BogackiShampine45 (const G4BogackiShampine45 &)=delete
G4BogackiShampine45 &	operator= (const G4BogackiShampine45 &)=delete
void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) override
void	SetupInterpolation ()
void	Interpolate (G4double tau, G4double yOut[])
void	SetupInterpolationHigh ()
void	InterpolateHigh (G4double tau, G4double yOut[]) const
G4double	DistChord () const override
G4int	IntegratorOrder () const override
void	GetLastDydx (G4double dyDxLast[])
void	PrepareConstants ()
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
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
G4int	IntegrationOrder ()
G4EquationOfMotion *	GetEquationOfMotion ()
const G4EquationOfMotion *	GetEquationOfMotion () const
void	SetEquationOfMotion (G4EquationOfMotion *newEquation)
unsigned long	GetfNoRHSCalls ()
void	ResetfNORHSCalls ()
G4bool	IsFSAL () const
G4bool	isQSS () const
void	SetIsQSS (G4bool val)

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

Detailed Description

Definition at line 51 of file G4BogackiShampine45.hh.

Constructor & Destructor Documentation

G4BogackiShampine45() [1/2]

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

Definition at line 70 of file G4BogackiShampine45.cc.

   : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
{
    const G4int numberOfVariables = noIntegrationVariables;
    
    // New Chunk of memory being created for use by the stepper
    
    // aki - for storing intermediate RHS
    ak2 = new G4double[numberOfVariables];
    ak3 = new G4double[numberOfVariables];
    ak4 = new G4double[numberOfVariables];
    ak5 = new G4double[numberOfVariables];
    ak6 = new G4double[numberOfVariables];
    ak7 = new G4double[numberOfVariables];
    ak8 = new G4double[numberOfVariables];
    ak9 = new G4double[numberOfVariables];
    ak10 = new G4double[numberOfVariables];
    ak11 = new G4double[numberOfVariables];    
 
    for (auto & i : p)
    {
       i= new G4double[numberOfVariables];
    }
 
    assert ( GetNumberOfStateVariables() >= 8 );    
    const G4int numStateVars = std::max(noIntegrationVariables,
                                        GetNumberOfStateVariables() );  
 
    // Must ensure space extra 'state' variables exists - i.e. yIn[7]
    yTemp = new G4double[numStateVars];
    yIn = new G4double[numStateVars] ;
    
    fLastInitialVector = new G4double[numStateVars] ;
    fLastFinalVector = new G4double[numStateVars] ;
    fLastDyDx = new G4double[numberOfVariables];  // Only derivatives
 
    fMidVector = new G4double[numberOfVariables];
    fMidError =  new G4double[numberOfVariables];
 
    if( ! fPreparedConstants )
    {
       PrepareConstants();
    }
    
    if( primary )
    {
       fAuxStepper = new G4BogackiShampine45(EqRhs, numberOfVariables, false);
    }
}

Referenced by G4BogackiShampine45().

~G4BogackiShampine45()

G4BogackiShampine45::~G4BogackiShampine45 ( )

override

Definition at line 124 of file G4BogackiShampine45.cc.

{
    // Clear all previously allocated memory for stepper and DistChord
    //
    delete [] ak2;
    delete [] ak3;
    delete [] ak4;
    delete [] ak5;
    delete [] ak6;
    delete [] ak7;
    delete [] ak8;
    delete [] ak9;
    delete [] ak10;
    delete [] ak11;
 
    for (auto & i : p)
    {
       delete [] i;
    }
 
    delete [] yTemp;
    delete [] yIn;
    
    delete [] fLastInitialVector;
    delete [] fLastFinalVector;
    delete [] fLastDyDx;
    delete [] fMidVector;
    delete [] fMidError;
    
    delete fAuxStepper;
}

G4BogackiShampine45() [2/2]

G4BogackiShampine45::G4BogackiShampine45 ( const G4BogackiShampine45 & )

delete

Member Function Documentation

DistChord()

G4double G4BogackiShampine45::DistChord ( ) const

overridevirtual

Implements G4MagIntegratorStepper.

Definition at line 310 of file G4BogackiShampine45.cc.

{
    G4double distLine, distChord;
    G4ThreeVector initialPoint, finalPoint, midPoint;
    
    // Store last initial and final points
    // (they will be overwritten in self-Stepper call!)
    //
    initialPoint = G4ThreeVector(fLastInitialVector[0],
                                 fLastInitialVector[1], fLastInitialVector[2]);
    finalPoint   = G4ThreeVector(fLastFinalVector[0],
                                 fLastFinalVector[1],  fLastFinalVector[2]);
 
#if 1
    // Old method -- Do half a step using StepNoErr
    //
    fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5*fLastStepLength,
                          fMidVector, fMidError);
#else
    // New method -- Using interpolation,
    // requires only 3 extra stages (ie 3 extra field evaluations )
 
    // Use Interpolation, instead of auxiliary stepper to evaluate midpoint
    //
    if( ! fPreparedInterpolation )
    {
       G4BogackiShampine45* cThis = const_cast<G4BogackiShampine45 *>(this);
       cThis-> SetupInterpolationHigh();
    }
    // For calculating the output at the tau fraction of Step
    //
    G4double tau = 0.5;
    InterpolateHigh( tau, fMidVector );    
#endif
   
    midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
    
    // Use stored values of Initial and Endpoint + new Midpoint to evaluate
    // distance of Chord
    
    if (initialPoint != finalPoint)
    {
        distLine  = G4LineSection::Distline( midPoint,initialPoint,finalPoint );
        distChord = distLine;
    }
    else
    {
        distChord = (midPoint-initialPoint).mag();
    }
    return distChord;
}

GetLastDydx()

void G4BogackiShampine45::GetLastDydx

(

G4double

dyDxLast[]

)

Definition at line 156 of file G4BogackiShampine45.cc.

{
   const G4int numberOfVariables = GetNumberOfVariables();
   
   for(G4int i=0; i < numberOfVariables; ++i )
   {
      dyDxLast[i] = ak9[i];
   }
}       

IntegratorOrder()

G4int G4BogackiShampine45::IntegratorOrder ( ) const

inlineoverridevirtual

Implements G4MagIntegratorStepper.

Definition at line 97 of file G4BogackiShampine45.hh.

{ return 4; }

Interpolate()

void G4BogackiShampine45::Interpolate ( G4double tau, G4double yOut[] )

inline

Definition at line 82 of file G4BogackiShampine45.hh.

    {
      InterpolateHigh( tau, yOut);       
      // InterpolateHigh( yInput, dydx, Step, yOut, tau);
    }

InterpolateHigh()

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

Definition at line 563 of file G4BogackiShampine45.cc.

{
    G4int numberOfVariables = GetNumberOfVariables();
    
    G4Exception("G4BogackiShampine45::InterpolateHigh()", "GeomField0001",
                FatalException, "Method is not yet validated.");
 
    // const G4double *yIn=  fLastInitialVector;
    // const G4double *dydx= fLastDyDx;
    const G4double Step = fLastStepLength;
    
#if 1
    G4int nwant = numberOfVariables;
    const G4int norder= 6;
    G4int l, k;
 
    for (l = 0; l < nwant; ++l)
    {
      yOut[l] = p[norder-1][l] * tau;
    }
    for (k = norder - 2; k >= 1; --k)
    {
      for (l = 0; l < nwant; ++l)
      {
         yOut[l] = ( yOut[l] + p[k][l] ) * tau;
      }
    }
    for (l = 0; l < nwant; ++l)
    {
      yOut[l] = ( yOut[l] + Step * ak8[l] ) * tau + yIn[l];
    }
    // The derivative at the end-point is nextDydx[i] = ak8[i];
#else
    // The scheme tries to do the same as the DormandPrince745 routine,
    // but fails 
 
    G4double b[12];
    const G4double* dydx = fLastDyDx;
    
    G4double tau0 = tau;
    
    for(G4int iStage=1; iStage<=11; ++iStage)   // iStage = stage number
    {
        b[iStage] = 0.0;
        tau = tau0;
        for(G4int j=6; j>=1; --j)               // j reversed
        {
            b[iStage] += bi[iStage][j] * tau;
            tau *= tau0;
        }
    }
 
    for(G4int i=0; i<numberOfVariables; ++i)
    {
        yOut[i] = yIn[i] + Step*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
                                 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
                                 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
                                 b[10]*ak10[i] + b[11]*ak11[i] );
    }
#endif    
}

Referenced by DistChord(), and Interpolate().

operator=()

G4BogackiShampine45 & G4BogackiShampine45::operator= ( const G4BogackiShampine45 & )

delete

PrepareConstants()

void G4BogackiShampine45::PrepareConstants

(

)

Definition at line 487 of file G4BogackiShampine45.cc.

{
    for(auto i=1; i<= 11; ++i)
    {
        bi[i][1] = 0.0 ;
    }
    
    for(auto i=1; i<=6; ++i)
    {
        bi[2][i] = 0.0 ;
    }
 
    bi[1][6] = -12134338393.0 / 1050809760.0 ,
    bi[1][5] =  -1620741229.0 / 50038560.0 ,
    bi[1][4] =  -2048058893.0 / 59875200.0 ,
    bi[1][3] = -87098480009.0 / 5254048800.0 ,
    bi[1][2] = -11513270273.0 / 3502699200.0 ,
    // 
    bi[3][6] =  -33197340367.0 / 1218433216.0 ,
    bi[3][5] = -539868024987.0 / 6092166080.0 ,
    bi[3][4] =  -39991188681.0 / 374902528.0 ,
    bi[3][3] =  -69509738227.0 / 1218433216.0 ,
    bi[3][2] =  -29327744613.0 / 2436866432.0 ,
    //
    bi[4][6] =   -284800997201.0 /  19905339168.0 ,
    bi[4][5] =  -7896875450471.0 / 165877826400.0 ,
    bi[4][4] =   -333945812879.0 /   5671036800.0 ,
    bi[4][3] = -16209923456237.0 / 497633479200.0 ,
    bi[4][2] =  -2382590741699.0 / 331755652800.0 ,
    //
    bi[5][6] = -540919.0 / 741312.0 ,
    bi[5][5] = -103626067.0 / 43243200.0 ,
    bi[5][4] = -633779.0 / 211200.0 ,
    bi[5][3] = -32406787.0 / 18532800.0 ,
    bi[5][2] = -36591193.0 / 86486400.0 ,
    //
    bi[6][6] = 7157998304.0 / 374350977.0 ,
    bi[6][5] = 30405842464.0 / 623918295.0 ,
    bi[6][4] = 183022264.0 / 5332635.0 ,
    bi[6][3] = -3357024032.0 / 1871754885.0 ,
    bi[6][2] = -611586736.0 / 89131185.0 ,
    //
    bi[7][6] =  -138073.0 /  9408.0 ,
    bi[7][5] =  -719433.0 / 15680.0 ,
    bi[7][4] = -1620541.0 / 31360.0 ,
    bi[7][3] =  -385151.0 / 15680.0 ,
    bi[7][2] =  -65403.0  / 15680.0 ,
    //
    bi[8][6] = 1245.0 / 64.0 ,
    bi[8][5] = 3991.0 / 64.0 ,
    bi[8][4] = 4715.0 / 64.0 ,
    bi[8][3] = 2501.0 / 64.0 ,
    bi[8][2] =  149.0 / 16.0 ,
    bi[8][1] = 1.0 ,
    //
    bi[9][6] = 55.0 / 3.0 ,
    bi[9][5] = 71.0 ,
    bi[9][4] = 103.0 ,
    bi[9][3] = 199.0 / 3.0 ,
    bi[9][2] = 16.0 ,
    //
    bi[10][6] =  -1774004627.0  /  75810735.0 ,
    bi[10][5] =  -1774004627.0  /  25270245.0 ,
    bi[10][4] =    -26477681.0  /    359975.0 ,
    bi[10][3] = -11411880511.0  / 379053675.0 ,
    bi[10][2] =   -423642896.0  / 126351225.0 ,
    //
    bi[11][6] =  35.0 ,
    bi[11][5] = 105.0 ,
    bi[11][4] = 117.0 ,
    bi[11][3] =  59.0 ,
    bi[11][2] =  12.0 ;
 
    fPreparedConstants = true;
}

Referenced by G4BogackiShampine45().

SetupInterpolation()

void G4BogackiShampine45::SetupInterpolation ( )

inline

Definition at line 75 of file G4BogackiShampine45.hh.

    {
      SetupInterpolationHigh(); // ( yInput, dydx, Step);
    }

SetupInterpolationHigh()

void G4BogackiShampine45::SetupInterpolationHigh

(

)

Definition at line 362 of file G4BogackiShampine45.cc.

{
    // Coefficients for the additional stages
    //
    const G4double
    a91 = 455.0/6144.0 ,
    a92 = 0.0 ,
    a93 = 10256301.0/35409920.0 ,
    a94 = 2307361.0/17971200.0 ,
    a95 = -387.0/102400.0 ,
    a96 = 73.0/5130.0 ,
    a97 = -7267.0/215040.0 ,
    a98 = 1.0/32.0 ,
    
    a101 = -837888343715.0/13176988637184.0 ,
    a102 = 30409415.0/52955362.0 ,
    a103 = -48321525963.0/759168069632.0 ,
    a104 = 8530738453321.0/197654829557760.0 ,
    a105 = 1361640523001.0/1626788720640.0 ,
    a106 = -13143060689.0/38604458898.0 ,
    a107 = 18700221969.0/379584034816.0 ,
    a108 = -5831595.0/847285792.0 ,
    a109 = -5183640.0/26477681.0 ,
    
    a111 = 98719073263.0/1551965184000.0 ,
    a112 = 1307.0/123552.0 ,
    a113 = 4632066559387.0/70181753241600.0 ,
    a114 = 7828594302389.0/382182512025600.0 ,
    a115 = 40763687.0/11070259200.0 ,
    a116 = 34872732407.0/224610586200.0 ,
    a117 = -2561897.0/30105600.0 ,
    a118 = 1.0/10.0 ,
    a119 = -1.0/10.0 ,
    a1110 = -1403317093.0/11371610250.0 ;
    
    const G4int numberOfVariables= this->GetNumberOfVariables();
    const G4double* dydx= fLastDyDx;
    const G4double Step = fLastStepLength;
    
    yTemp[7]  = yIn[7];
 
    // Evaluate the extra stages
    //
    for(G4int i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
                                  a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
                                  a97*ak7[i] + a98*ak8[i] );
    }
    
    RightHandSide(yTemp, ak9);   // 9th stage
    
    for(G4int i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
                                  a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
                                  a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
    }
    
    RightHandSide(yTemp, ak10);  // 10th stage
    
    for(G4int i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
                                  a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
                                  a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
                                  a1110*ak10[i] );
    }
    RightHandSide(yTemp, ak11);  // 11th stage
 
    //  In future we can restrict the number of variables interpolated
    //
    G4int nwant = numberOfVariables; 
 
    //  Form the coefficients of the interpolating polynomial in its shifted
    //  and scaled form.  The terms are grouped to minimize the errors 
    //  of the transformation, to cope with ill-conditioning. ( From RKSUITE )
    //   
    for (G4int l = 0; l < nwant; ++l)
    {
        //  Coefficient of tau^6        
        p[5][l] =   bi[5][6]*ak5[l] +
                  ((bi[10][6]*ak10[l] + bi[8][6]*ak8[l]) + 
                   (bi[7][6]*ak7[l] + bi[6][6]*ak6[l]))  + 
                  ((bi[4][6]*ak4[l] + bi[9][6]*ak9[l]) + 
                   (bi[3][6]*ak3[l] + bi[11][6]*ak11[l]) + 
                    bi[1][6]*dydx[l]);
        //  Coefficient of tau^5
        p[4][l] = (bi[10][5]*ak10[l] + bi[9][5]*ak9[l])  + 
                 ((bi[7][5]*ak7[l] + bi[6][5]*ak6[l]) + 
                   bi[5][5]*ak5[l])  +  ((bi[4][5]*ak4[l] + 
                   bi[8][5]*ak8[l]) + (bi[3][5]*ak3[l] + 
                   bi[11][5]*ak11[l]) + bi[1][5]*dydx[l]);
        //  Coefficient of tau^4
        p[3][l] = ((bi[4][4]*ak4[l] + bi[8][4]*ak8[l]) + 
                   (bi[7][4]*ak7[l] + bi[6][4]*ak6[l]) + 
                    bi[5][4]*ak5[l]) + ((bi[10][4]*ak10[l] + 
                    bi[9][4]*ak9[l]) +  (bi[3][4]*ak3[l] + 
                    bi[11][4]*ak11[l]) + bi[1][4]*dydx[l]);
        //  Coefficient of tau^3
        p[2][l] =  bi[5][3]*ak5[l] + bi[6][3]*ak6[l]  + 
                 ((bi[3][3]*ak3[l] + bi[9][3]*ak9[l]) + 
                 (bi[10][3]*ak10[l]+ bi[8][3]*ak8[l]) + bi[1][3]*dydx[l]) + 
                 ((bi[4][3]*ak4[l] + bi[11][3]*ak11[l]) + bi[7][3]*ak7[l]);
        //  Coefficient of tau^2
        p[1][l] = bi[5][2]*ak5[l]  + ((bi[6][2]*ak6[l] + 
                  bi[8][2]*ak8[l]) +   bi[1][2]*dydx[l])  + 
                ((bi[3][2]*ak3[l]  +   bi[9][2]*ak9[l]) + 
                 bi[10][2]*ak10[l])+ ((bi[4][2]*ak4[l] + 
                 bi[11][2]*ak2[l]) +   bi[7][2]*ak7[l]);
      }
 
      //  Scale all the coefficients by the step size.
      //
      for (auto & i : p)
      {
         for (G4int l = 0; l < nwant; ++l)
         {  
            i[l] *= Step;
         }
      }
 
    fPreparedInterpolation = true;
}

Referenced by DistChord(), and SetupInterpolation().

Stepper()

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

overridevirtual

Implements G4MagIntegratorStepper.

Definition at line 171 of file G4BogackiShampine45.cc.

{
    G4int i;
    
    // Constants from the Butcher tableu
    //
    const G4double 
       b21 = 1.0/6.0 ,
       b31 = 2.0/27.0 ,     b32 = 4.0/27.0,
       
       b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
       
       b51 = 68.0/297.0,    b52 = -4.0/11.0,
       b53 = 42.0/143.0,    b54 = 1960.0/3861.0,
       
       b61 = 597.0/22528.0,    b62 = 81.0/352.0,
       b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
       b65 = 4617.0/20480.0,
       
       b71 = 174197.0/959244.0,    b72 = -30942.0/79937.0,
       b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
       b75 = -29421.0/29068.0,     b76 = 482048.0/414219.0,
       
       b81 = 587.0/8064.0,         b82 = 0.0,
       b83 = 4440339.0/15491840.0, b84 = 24353.0/124800.0,
       b85 = 387.0/44800.0,        b86 = 2152.0/5985.0,
       b87 = 7267.0/94080.0;
    
//    c1 = 2479.0/34992.0,
//    c2 = 0.0,
//    c3 = 123.0/416.0,
//    c4 = 612941.0/3411720.0,
//    c5 = 43.0/1440.0,
//    c6 = 2272.0/6561.0,
//    c7 = 79937.0/1113912.0,
//    c8 = 3293.0/556956.0,
    
    // For the embedded higher order method only the difference of values
    // taken and is used directly later (instead of defining the last row
    // of Butcher table in separate constants and taking the
    // difference)
    //
    const G4double 
       dc1 = b81 -   2479.0 /   34992.0 ,
       dc2 = 0.0,
       dc3 = b83 -    123.0 /     416.0 ,
       dc4 = b84 - 612941.0 / 3411720.0,
       dc5 = b85 -     43.0 /    1440.0,
       dc6 = b86 -   2272.0 /    6561.0,
       dc7 = b87 -  79937.0 / 1113912.0,
       dc8 =     -   3293.0 /  556956.0;
 
    const G4int numberOfVariables = GetNumberOfVariables();
    
    // The number of variables to be integrated over
    //
    yOut[7] = yTemp[7]  = yIn[7] = yInput[7];
 
    // Saving yInput because yInput and yOut can be aliases for same array
    //
    for(i=0; i<numberOfVariables; ++i)
    {
        yIn[i]=yInput[i];
    }
    
    // RightHandSide(yIn, dydx) ;
    // 1st Step - Not doing, getting passed
    //
    for(i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
    }
    RightHandSide(yTemp, ak2) ;              // 2nd Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
    }
    RightHandSide(yTemp, ak3) ;              // 3rd Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
    }
    RightHandSide(yTemp, ak4) ;              // 4th Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
                                  b54*ak4[i]) ;
    }
    RightHandSide(yTemp, ak5) ;              // 5th Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
                                  b64*ak4[i] + b65*ak5[i]) ;
    }
    RightHandSide(yTemp, ak6) ;              // 6th Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
                                  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
    }
    RightHandSide(yTemp, ak7);               // 7th Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] +
                                  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
                                  b87*ak7[i]);
    }
    RightHandSide(yOut, ak8);                // 8th Step - Final one Using FSAL
 
    for(i=0; i<numberOfVariables; ++i)
    {
        yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
                        dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
        
        // Store Input and Final values, for possible use in calculating chord
        //
        fLastInitialVector[i] = yIn[i] ;
        fLastFinalVector[i]   = yOut[i];
        fLastDyDx[i]          = DyDx[i];
    }
    
    fLastStepLength = Step;
    fPreparedInterpolation= false;
 
    return ;
}

Referenced by DistChord().

---

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

• G4BogackiShampine45.hh
• G4BogackiShampine45.cc