G4NystromRK4 Class Reference

#include <G4NystromRK4.hh>

Inheritance diagram for G4NystromRK4:

Public Member Functions
	G4NystromRK4 (G4Mag_EqRhs *EquationMotion, G4double distanceConstField=0.0)
	~G4NystromRK4 ()
virtual void	Stepper (const G4double y[], const G4double dydx[], G4double hstep, G4double yOut[], G4double yError[]) override
void	SetDistanceForConstantField (G4double length)
G4double	GetDistanceForConstantField () const
virtual G4int	IntegratorOrder () const override
virtual 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
virtual void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])=0
virtual G4double	DistChord () const =0
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
virtual G4int	IntegratorOrder () const =0
G4int	IntegrationOrder ()
G4EquationOfMotion *	GetEquationOfMotion ()
const G4EquationOfMotion *	GetEquationOfMotion () const
void	SetEquationOfMotion (G4EquationOfMotion *newEquation)
unsigned long	GetfNoRHSCalls ()
void	ResetfNORHSCalls ()
G4bool	IsFSAL () const

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 51 of file G4NystromRK4.cc.

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

~G4NystromRK4()

G4NystromRK4::~G4NystromRK4 ( )

inline

Definition at line 58 of file G4NystromRK4.hh.

{}

Member Function Documentation

DistChord()

G4double G4NystromRK4::DistChord ( ) const

overridevirtual

Implements G4MagIntegratorStepper.

Definition at line 183 of file G4NystromRK4.cc.

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

GetDistanceForConstantField()

G4double G4NystromRK4::GetDistanceForConstantField

(

)

const

Definition at line 206 of file G4NystromRK4.cc.

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

IntegratorOrder()

virtual 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 188 of file G4NystromRK4.cc.

{
  if (GetField() == nullptr)
  {
    G4Exception("G4NystromRK4::SetDistanceForConstantField",
                "Nystrom 001", JustWarning,
        "Provided field is not G4CachedMagneticField. Changing field type.");
 
    fCachedField = std::unique_ptr<G4CachedMagneticField>(
       new 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 60 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