G4BogackiShampine45 Class Reference

G4BogackiShampine45 is an integrator of particle's equation of motion based on the Bogacki-Shampine method with FSAL property, allowing the reuse of the last derivative in the next step. This Stepper provides 'dense output'. After a successful step, it is possible to obtain an estimate of the value of the function at an intermediate point of the interval. This requires only two additional evaluations of the derivative (and thus the field). More...

#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	SetupInterpolationHigh ()
void	SetupInterpolation ()
void	InterpolateHigh (G4double tau, G4double yOut[]) const
void	Interpolate (G4double tau, G4double yOut[])
G4double	DistChord () const override
G4int	IntegratorOrder () const override
G4StepperType	StepperType () 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

G4BogackiShampine45 is an integrator of particle's equation of motion based on the Bogacki-Shampine method with FSAL property, allowing the reuse of the last derivative in the next step. This Stepper provides 'dense output'. After a successful step, it is possible to obtain an estimate of the value of the function at an intermediate point of the interval. This requires only two additional evaluations of the derivative (and thus the field).

Definition at line 61 of file G4BogackiShampine45.hh.

Constructor & Destructor Documentation

G4BogackiShampine45() [1/2]

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

Constructor for G4BogackiShampine45.

Parameters

[in]	EqRhs	Pointer to the provided equation of motion.
[in]	numberOfVariables	The number of integration variables.
[in]	primary	Flag for initialisation of the auxiliary stepper.

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 DistChord(), G4BogackiShampine45(), G4BogackiShampine45(), and operator=().

~G4BogackiShampine45()

G4BogackiShampine45::~G4BogackiShampine45 ( )

override

Destructor.

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

Copy constructor and assignment operator not allowed.

Member Function Documentation

DistChord()

G4double G4BogackiShampine45::DistChord ( ) const

overridevirtual

Returns the distance from chord line.

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[]

)

Acccessor for dydx array.

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

Returns the order, 4, of integration.

Implements G4MagIntegratorStepper.

Definition at line 126 of file G4BogackiShampine45.hh.

{ return 4; }

Interpolate()

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

inline

Definition at line 115 of file G4BogackiShampine45.hh.

                                               { InterpolateHigh( tau, yOut); }

InterpolateHigh()

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

Calculates the output at the tau fraction of step.

Parameters

[in]	tau	The tau fraction of the step.
[out]	yOut	Interpolation output.

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

(

)

Initialises the values of the bi[][] array.

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 107 of file G4BogackiShampine45.hh.

{ SetupInterpolationHigh(); }

SetupInterpolationHigh()

void G4BogackiShampine45::SetupInterpolationHigh

(

)

Setup all coefficients for interpolation.

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

The stepper for the Runge Kutta integration. The stepsize is fixed, with the step size given by 'h'. Integrates ODE starting values y[0 to 6]. Outputs yout[] and its estimated error yerr[].

Parameters

[in]	y	Starting values array of integration variables.
[in]	dydx	Derivatives array.
[in]	h	The given step size.
[out]	yout	Integration output.
[out]	yerr	The estimated error.

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 ;
}

StepperType()

G4StepperType G4BogackiShampine45::StepperType ( ) const

inlineoverridevirtual

Returns the stepper type-ID, "kBogackiShampine45".

Reimplemented from G4MagIntegratorStepper.

Definition at line 131 of file G4BogackiShampine45.hh.

{ return kBogackiShampine45; }

---

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

• G4BogackiShampine45.hh
• G4BogackiShampine45.cc