G4ConstRK4 Class Reference

G4ConstRK4 performs the integration of one step with error calculation in constant magnetic field. The integration method is the same as in ClassicalRK4. The field value is assumed constant for the step. This field evaluation is called only once per step. G4ConstRK4 can be used only for magnetic fields. More...

#include <G4ConstRK4.hh>

Inheritance diagram for G4ConstRK4:

Public Member Functions
	G4ConstRK4 (G4Mag_EqRhs *EquationMotion, G4int numberOfStateVariables=8)
	~G4ConstRK4 () override
	G4ConstRK4 (const G4ConstRK4 &)=delete
G4ConstRK4 &	operator= (const G4ConstRK4 &)=delete
void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) override
void	DumbStepper (const G4double yIn[], const G4double dydx[], G4double h, G4double yOut[]) override
G4double	DistChord () const override
void	RightHandSideConst (const G4double y[], G4double dydx[]) const
void	GetConstField (const G4double y[], G4double Field[])
G4int	IntegratorOrder () const override
G4StepperType	StepperType () const override
Public Member Functions inherited from G4MagErrorStepper
	G4MagErrorStepper (G4EquationOfMotion *EqRhs, G4int numberOfVariables, G4int numStateVariables=12)
	~G4MagErrorStepper () override
	G4MagErrorStepper (const G4MagErrorStepper &)=delete
G4MagErrorStepper &	operator= (const G4MagErrorStepper &)=delete
void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) override
G4double	DistChord () const override
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

G4ConstRK4 performs the integration of one step with error calculation in constant magnetic field. The integration method is the same as in ClassicalRK4. The field value is assumed constant for the step. This field evaluation is called only once per step. G4ConstRK4 can be used only for magnetic fields.

Definition at line 53 of file G4ConstRK4.hh.

Constructor & Destructor Documentation

G4ConstRK4() [1/2]

G4ConstRK4::G4ConstRK4	(	G4Mag_EqRhs *	EquationMotion,
		G4int	numberOfStateVariables = 8 )

Constructor for G4ConstRK4.

Parameters

[in]	EqRhs	Pointer to the provided equation of motion.
[in]	numberOfVariables	The number of integration variables.

Definition at line 40 of file G4ConstRK4.cc.

  : G4MagErrorStepper(EqRhs, 6, numStateVariables)
{
  // const G4int numberOfVariables= 6;
  if( numStateVariables < 8 ) 
  {
    std::ostringstream message;
    message << "The number of State variables at least 8 " << G4endl
            << "Instead it is - numStateVariables= " << numStateVariables;
    G4Exception("G4ConstRK4::G4ConstRK4()", "GeomField0002",
                FatalException, message, "Use another Stepper!");
  }
 
  fEq = EqRhs;
  yMiddle  = new G4double[8];
  dydxMid  = new G4double[8];
  yInitial = new G4double[8];
  yOneStep = new G4double[8];
 
  dydxm = new G4double[8];
  dydxt = new G4double[8]; 
  yt    = new G4double[8]; 
  Field[0]=0.; Field[1]=0.; Field[2]=0.;
}

Referenced by G4ConstRK4(), and operator=().

~G4ConstRK4()

G4ConstRK4::~G4ConstRK4 ( )

override

Destructor.

Definition at line 69 of file G4ConstRK4.cc.

{
   delete [] yMiddle;
   delete [] dydxMid;
   delete [] yInitial;
   delete [] yOneStep;
   delete [] dydxm;
   delete [] dydxt;
   delete [] yt;
}

G4ConstRK4() [2/2]

G4ConstRK4::G4ConstRK4 ( const G4ConstRK4 & )

delete

Copy constructor and assignment operator not allowed.

Member Function Documentation

DistChord()

G4double G4ConstRK4::DistChord ( ) const

overridevirtual

Returns the distance from chord line.

Implements G4MagIntegratorStepper.

Definition at line 211 of file G4ConstRK4.cc.

{
  G4double distLine, distChord; 
 
  if (fInitialPoint != fFinalPoint)
  {
     distLine= G4LineSection::Distline( fMidPoint, fInitialPoint, fFinalPoint );
       // This is a class method that gives distance of Mid 
       // from the Chord between the Initial and Final points
     distChord = distLine;
  }
  else
  {
     distChord = (fMidPoint-fInitialPoint).mag();
  }
  return distChord;
}

DumbStepper()

void G4ConstRK4::DumbStepper ( const G4double yIn[], const G4double dydx[], G4double h, G4double yOut[] )

overridevirtual

Given values for the variables y[0,..,n-1] and their derivatives dydx[0,...,n-1] known at x, uses the classical 4th Runge-Kutta method to advance the solution over an interval h and returns the incremented variables as yout[0,...,n-1]. The user supplies the function RightHandSide(x,y,dydx), which returns derivatives dydx at x. The source is routine rk4 from NRC p.712-713.

Parameters

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

Implements G4MagErrorStepper.

Definition at line 90 of file G4ConstRK4.cc.

{
   G4double  hh = h*0.5 , h6 = h/6.0  ;
   
   // 1st Step K1=h*dydx
   yt[5] = yIn[5] + hh*dydx[5] ;
   yt[4] = yIn[4] + hh*dydx[4] ;
   yt[3] = yIn[3] + hh*dydx[3] ;
   yt[2] = yIn[2] + hh*dydx[2] ;
   yt[1] = yIn[1] + hh*dydx[1] ;
   yt[0] = yIn[0] + hh*dydx[0] ;
   RightHandSideConst(yt,dydxt) ;        
 
   // 2nd Step K2=h*dydxt
   yt[5] = yIn[5] + hh*dydxt[5] ;
   yt[4] = yIn[4] + hh*dydxt[4] ;
   yt[3] = yIn[3] + hh*dydxt[3] ;
   yt[2] = yIn[2] + hh*dydxt[2] ;
   yt[1] = yIn[1] + hh*dydxt[1] ;
   yt[0] = yIn[0] + hh*dydxt[0] ;
   RightHandSideConst(yt,dydxm) ;     
 
   // 3rd Step K3=h*dydxm
   // now dydxm=(K2+K3)/h
   yt[5] = yIn[5] + h*dydxm[5] ;
   dydxm[5] += dydxt[5] ;  
   yt[4] = yIn[4] + h*dydxm[4] ;
   dydxm[4] += dydxt[4] ;  
   yt[3] = yIn[3] + h*dydxm[3] ;
   dydxm[3] += dydxt[3] ;  
   yt[2] = yIn[2] + h*dydxm[2] ;
   dydxm[2] += dydxt[2] ;  
   yt[1] = yIn[1] + h*dydxm[1] ;
   dydxm[1] += dydxt[1] ;  
   yt[0] = yIn[0] + h*dydxm[0] ;
   dydxm[0] += dydxt[0] ;  
   RightHandSideConst(yt,dydxt) ;   
 
   // 4th Step K4=h*dydxt
   yOut[5] = yIn[5]+h6*(dydx[5]+dydxt[5]+2.0*dydxm[5]);
   yOut[4] = yIn[4]+h6*(dydx[4]+dydxt[4]+2.0*dydxm[4]);
   yOut[3] = yIn[3]+h6*(dydx[3]+dydxt[3]+2.0*dydxm[3]);
   yOut[2] = yIn[2]+h6*(dydx[2]+dydxt[2]+2.0*dydxm[2]);
   yOut[1] = yIn[1]+h6*(dydx[1]+dydxt[1]+2.0*dydxm[1]);
   yOut[0] = yIn[0]+h6*(dydx[0]+dydxt[0]+2.0*dydxm[0]);
   
}  // end of DumbStepper ....................................................

Referenced by Stepper().

GetConstField()

void G4ConstRK4::GetConstField ( const G4double y[], G4double Field[] )

inline

Returns the field values, at position and time 'y'.

Parameters

[in]	y	The position vector plus time (x,y,z,t).
[out]	Field	The field value in output.

Definition at line 169 of file G4ConstRK4.hh.

{
  G4double PositionAndTime[4];
 
  PositionAndTime[0] = y[0];
  PositionAndTime[1] = y[1];
  PositionAndTime[2] = y[2];
  // Global Time
  PositionAndTime[3] = y[7];
  fEq -> GetFieldValue(PositionAndTime, B);
}

Referenced by Stepper().

IntegratorOrder()

G4int G4ConstRK4::IntegratorOrder ( ) const

inlineoverridevirtual

Returns the order, 4, of integration.

Implements G4MagIntegratorStepper.

Definition at line 132 of file G4ConstRK4.hh.

{ return 4; }

Referenced by Stepper().

operator=()

G4ConstRK4 & G4ConstRK4::operator= ( const G4ConstRK4 & )

delete

RightHandSideConst()

void G4ConstRK4::RightHandSideConst ( const G4double y[], G4double dydx[] ) const

inline

Returns the derivatives value, at position and time 'y'.

Parameters

[in]	y	The position vector plus time (x,y,z,t).
[out]	dydx	The derivatives array.

Definition at line 151 of file G4ConstRK4.hh.

{
  
  G4double momentum_mag_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5];
  G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
    
  G4double cof = fEq->FCof()*inv_momentum_magnitude;
 
  dydx[0] = y[3]*inv_momentum_magnitude;       //  (d/ds)x = Vx/V
  dydx[1] = y[4]*inv_momentum_magnitude;       //  (d/ds)y = Vy/V
  dydx[2] = y[5]*inv_momentum_magnitude;       //  (d/ds)z = Vz/V
 
  dydx[3] = cof*(y[4]*Field[2] - y[5]*Field[1]) ;   // Ax = a*(Vy*Bz - Vz*By)
  dydx[4] = cof*(y[5]*Field[0] - y[3]*Field[2]) ;   // Ay = a*(Vz*Bx - Vx*Bz)
  dydx[5] = cof*(y[3]*Field[1] - y[4]*Field[0]) ;   // Az = a*(Vx*By - Vy*Bx)
}

Referenced by DumbStepper(), and Stepper().

Stepper()

void G4ConstRK4::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 146 of file G4ConstRK4.cc.

{
   const G4int nvar = 6;  // number of variables integrated
   const G4int maxvar = GetNumberOfStateVariables();
 
   // Correction for Richardson extrapolation
   G4double  correction = 1. / ( (1 << IntegratorOrder()) -1 );
 
   G4int i;
   
   // Saving yInput because yInput and yOutput can be aliases for same array
   for (i=0;    i<maxvar; ++i) { yInitial[i]= yInput[i]; }
 
   // Must copy the part of the state *not* integrated to the output
   for (i=nvar; i<maxvar; ++i) { yOutput[i]=  yInput[i]; }
 
   // yInitial[7]= yInput[7];  //  The time is typically needed
   yMiddle[7]  = yInput[7];   // Copy the time from initial value 
   yOneStep[7] = yInput[7];   // As it contributes to final value of yOutput ?
   // yOutput[7] = yInput[7];  // -> dumb stepper does it too for RK4
   yError[7] = 0.0;         
 
   G4double halfStep = hstep * 0.5; 
 
   // Do two half steps
   //
   GetConstField(yInitial,Field);
   DumbStepper  (yInitial,  dydx,   halfStep, yMiddle);
   RightHandSideConst(yMiddle, dydxMid);    
   DumbStepper  (yMiddle, dydxMid, halfStep, yOutput); 
 
   // Store midpoint, chord calculation
   //
   fMidPoint = G4ThreeVector( yMiddle[0],  yMiddle[1],  yMiddle[2]); 
 
   // Do a full Step
   //
   DumbStepper(yInitial, dydx, hstep, yOneStep);
   for(i=0; i<nvar; ++i)
   {
      yError [i] = yOutput[i] - yOneStep[i] ;
      yOutput[i] += yError[i]*correction ;
        // Provides accuracy increased by 1 order via the 
        // Richardson extrapolation  
   }
 
   fInitialPoint = G4ThreeVector( yInitial[0], yInitial[1], yInitial[2]); 
   fFinalPoint   = G4ThreeVector( yOutput[0],  yOutput[1],  yOutput[2]); 
 
   return;
}

StepperType()

G4StepperType G4ConstRK4::StepperType ( ) const

inlineoverridevirtual

Returns the stepper type-ID, "kConstRK4".

Reimplemented from G4MagIntegratorStepper.

Definition at line 137 of file G4ConstRK4.hh.

{ return kConstRK4; }

---

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

• G4ConstRK4.hh
• G4ConstRK4.cc