G4NystromRK4 Class Reference

#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
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 G4NystromRK4.hh.

Constructor & Destructor Documentation

G4NystromRK4()

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

Definition at line 53 of file G4NystromRK4.cc.

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

~G4NystromRK4()

G4NystromRK4::~G4NystromRK4 ( )

overridedefault

Member Function Documentation

DistChord()

G4double G4NystromRK4::DistChord ( ) const

overridevirtual

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

Implements G4MagIntegratorStepper.

Definition at line 71 of file G4NystromRK4.hh.

{ return 4; }

SetDistanceForConstantField()

void G4NystromRK4::SetDistanceForConstantField

(

G4double

length

)

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

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

---

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

• G4NystromRK4.hh
• G4NystromRK4.cc