G4NystromRK4 Class Reference

G4NystromRK4 integrates the equations of the motion of a particle in a magnetic field using 4th Runge-Kutta-Nystrom method with errors estimation. The current form can be used only for 'pure' magnetic field. More...

#include <G4NystromRK4.hh>

Inheritance diagram for G4NystromRK4:

Public Member Functions
	G4NystromRK4 (G4Mag_EqRhs *EquationMotion, G4double distanceConstField=0.0)
	~G4NystromRK4 () override=default
void	Stepper (const G4double y[], const G4double dydx[], G4double hstep, G4double yOut[], G4double yError[]) override
void	SetDistanceForConstantField (G4double length)
G4double	GetDistanceForConstantField () const
G4int	IntegratorOrder () const override
G4double	DistChord () const override
G4StepperType	StepperType () 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

G4NystromRK4 integrates the equations of the motion of a particle in a magnetic field using 4th Runge-Kutta-Nystrom method with errors estimation. The current form can be used only for 'pure' magnetic field.

Definition at line 56 of file G4NystromRK4.hh.

Constructor & Destructor Documentation

G4NystromRK4()

G4NystromRK4::G4NystromRK4	(	G4Mag_EqRhs *	EquationMotion,
		G4double	distanceConstField = 0.0 )

Constructor for G4NystromRK4. Can be used only for Magnetic Fields and for 6 variables (x,p).

Parameters

[in]	EquationMotion	Pointer to the provided equation of motion.
[in]	distanceConstField	Distance value for constant field.

Definition at line 53 of file G4NystromRK4.cc.

  : G4MagIntegratorStepper(equation, INTEGRATED_COMPONENTS)
{
    if (distanceConstField > 0)
    {
        SetDistanceForConstantField(distanceConstField);
    }
}

~G4NystromRK4()

G4NystromRK4::~G4NystromRK4 ( )

overridedefault

Default Destructor.

Member Function Documentation

DistChord()

G4double G4NystromRK4::DistChord ( ) const

overridevirtual

Returns the distance from chord line.

Implements G4MagIntegratorStepper.

Definition at line 185 of file G4NystromRK4.cc.

{
    return G4LineSection::Distline(fMidPoint, fInitialPoint, fEndPoint);
}

GetDistanceForConstantField()

G4double G4NystromRK4::GetDistanceForConstantField

(

)

const

Definition at line 207 of file G4NystromRK4.cc.

{
  if (GetField() == nullptr)
  {
    return 0.0;
  }
  return GetField()->GetConstDistance(); 
}

IntegratorOrder()

G4int G4NystromRK4::IntegratorOrder ( ) const

inlineoverridevirtual

Returns the order, 4, of integration.

Implements G4MagIntegratorStepper.

SetDistanceForConstantField()

void G4NystromRK4::SetDistanceForConstantField

(

G4double

length

)

Setter and getter for the distance value for constant field.

Definition at line 190 of file G4NystromRK4.cc.

{
  if (GetField() == nullptr)
  {
    G4Exception("G4NystromRK4::SetDistanceForConstantField",
                "Nystrom 001", JustWarning,
        "Provided field is not G4CachedMagneticField. Changing field type.");
 
    fCachedField = std::make_unique<G4CachedMagneticField>(
           dynamic_cast<G4MagneticField*>(GetEquationOfMotion()->GetFieldObj()),
           length);
 
    GetEquationOfMotion()->SetFieldObj(fCachedField.get());
  }
  GetField()->SetConstDistance(length);
}

Referenced by G4NystromRK4().

Stepper()

void G4NystromRK4::Stepper ( const G4double y[], const G4double dydx[], G4double hstep, G4double yOut[], G4double yError[] )

overridevirtual

The stepper for the Runge Kutta integration. The stepsize is fixed, with the step size given by 'hstep'. Integrates ODE starting values y[0 to 6]. Outputs yOut[] and its estimated error yError[]. Provides error via analytical method.

Parameters

[in]	y	Starting values array of integration variables.
[in]	dydx	Derivatives array.
[in]	hstep	The given step size.
[out]	yOut	Integration output.
[out]	yError	The estimated error.

Implements G4MagIntegratorStepper.

Definition at line 62 of file G4NystromRK4.cc.

{
    fInitialPoint = { y[0], y[1], y[2] };
 
    G4double field[3];
 
    constexpr G4double one_sixth= 1./6.;
    const G4double S5 = 0.5 * hstep;
    const G4double S4 = .25 * hstep;
    const G4double S6 = hstep * one_sixth;
   
    const G4double momentum2 = getValue2(y, Value3D::Momentum);
    if (notEquals(momentum2, fMomentum2))
    {
        fMomentum = std::sqrt(momentum2);
        fMomentum2 = momentum2;
        fInverseMomentum  = 1. / fMomentum;
        fCoefficient = GetFCof() * fInverseMomentum;
    }
    
    // Point 1
    const G4double K1[3] = { 
        fInverseMomentum * dydx[3], 
        fInverseMomentum * dydx[4], 
        fInverseMomentum * dydx[5]
    };
    
    // Point2
    G4double p[4] = {
        y[0] + S5 * (dydx[0] + S4 * K1[0]),
        y[1] + S5 * (dydx[1] + S4 * K1[1]),
        y[2] + S5 * (dydx[2] + S4 * K1[2]),
        y[7]
    };
 
    GetFieldValue(p, field);
  
    const G4double A2[3] = {
        dydx[0] + S5 * K1[0],
        dydx[1] + S5 * K1[1],
        dydx[2] + S5 * K1[2]
    };
 
    const G4double K2[3] = {
        (A2[1] * field[2] - A2[2] * field[1]) * fCoefficient,
        (A2[2] * field[0] - A2[0] * field[2]) * fCoefficient,
        (A2[0] * field[1] - A2[1] * field[0]) * fCoefficient
    };
 
    fMidPoint = { p[0], p[1], p[2] };
 
    // Point 3 with the same magnetic field
    const G4double A3[3] = {
        dydx[0] + S5 * K2[0],
        dydx[1] + S5 * K2[1],
        dydx[2] + S5 * K2[2]
    };
 
    const G4double K3[3] = {
        (A3[1] * field[2] - A3[2] * field[1]) * fCoefficient,
        (A3[2] * field[0] - A3[0] * field[2]) * fCoefficient,
        (A3[0] * field[1] - A3[1] * field[0]) * fCoefficient
    };
 
    // Point 4
    p[0] = y[0] + hstep * (dydx[0] + S5 * K3[0]);
    p[1] = y[1] + hstep * (dydx[1] + S5 * K3[1]);
    p[2] = y[2] + hstep * (dydx[2] + S5 * K3[2]);             
 
    GetFieldValue(p, field);
  
    const G4double A4[3] = {
        dydx[0] + hstep * K3[0],
        dydx[1] + hstep * K3[1],
        dydx[2] + hstep * K3[2]
    };
 
    const G4double K4[3] = {
        (A4[1] * field[2] - A4[2] * field[1]) * fCoefficient,
        (A4[2] * field[0] - A4[0] * field[2]) * fCoefficient,
        (A4[0] * field[1] - A4[1] * field[0]) * fCoefficient
    };
  
    // New position
    yOut[0] = y[0] + hstep * (dydx[0] + S6 * (K1[0] + K2[0] + K3[0]));
    yOut[1] = y[1] + hstep * (dydx[1] + S6 * (K1[1] + K2[1] + K3[1]));
    yOut[2] = y[2] + hstep * (dydx[2] + S6 * (K1[2] + K2[2] + K3[2]));
    // New direction
    yOut[3] = dydx[0] + S6 * (K1[0] + K4[0] + 2. * (K2[0] + K3[0]));
    yOut[4] = dydx[1] + S6 * (K1[1] + K4[1] + 2. * (K2[1] + K3[1]));
    yOut[5] = dydx[2] + S6 * (K1[2] + K4[2] + 2. * (K2[2] + K3[2]));
    // Pass Energy, time unchanged -- time is not integrated !!
    yOut[6] = y[6]; 
    yOut[7] = y[7];
 
    fEndPoint = { yOut[0], yOut[1], yOut[2] };
 
    // Errors
    yError[3] = hstep * std::fabs(K1[0] - K2[0] - K3[0] + K4[0]);
    yError[4] = hstep * std::fabs(K1[1] - K2[1] - K3[1] + K4[1]);
    yError[5] = hstep * std::fabs(K1[2] - K2[2] - K3[2] + K4[2]);
    yError[0] = hstep * yError[3];
    yError[1] = hstep * yError[4];
    yError[2] = hstep * yError[5];
    yError[3] *= fMomentum;
    yError[4] *= fMomentum;
    yError[5] *= fMomentum;
 
    // Normalize momentum
    const G4double normF = fMomentum / getValue(yOut, Value3D::Momentum);
    yOut[3] *= normF; 
    yOut[4] *= normF; 
    yOut[5] *= normF; 
 
    // My trial code:
    // G4double endMom2 = yOut[3]*yOut[3]+yOut[4]*yOut[4]+yOut[5]*yOut[5];
    // G4double normF = std::sqrt( startMom2 / endMom2 );  
}

StepperType()

G4StepperType G4NystromRK4::StepperType ( ) const

inlineoverridevirtual

Returns the stepper type-ID, "kNystromRK4".

Reimplemented from G4MagIntegratorStepper.

---

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

• G4NystromRK4.hh
• G4NystromRK4.cc