G4TsitourasRK45 Class Reference

#include <G4TsitourasRK45.hh>

Inheritance diagram for G4TsitourasRK45:

Collaboration diagram for G4TsitourasRK45:

Public Member Functions
	G4TsitourasRK45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
	~G4TsitourasRK45 ()
void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
void	SetupInterpolation ()
void	Interpolate (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
void	interpolate (const G4double yInput[], const G4double dydx[], G4double yOut[], G4double Step, G4double tau)
G4double	DistChord () const
G4int	IntegratorOrder () const
Public Member Functions inherited from G4MagIntegratorStepper
	G4MagIntegratorStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, bool isFSAL=false)
virtual	~G4MagIntegratorStepper ()
virtual void	ComputeRightHandSide (const G4double y[], G4double dydx[])
void	NormaliseTangentVector (G4double vec[6])
void	NormalisePolarizationVector (G4double vec[12])
void	RightHandSide (const double y[], double dydx[])
G4int	GetNumberOfVariables () const
G4int	GetNumberOfStateVariables () const
G4int	IntegrationOrder ()
G4EquationOfMotion *	GetEquationOfMotion ()
void	SetEquationOfMotion (G4EquationOfMotion *newEquation)
unsigned long	GetfNoRHSCalls ()
void	ResetfNORHSCalls ()
bool	IsFSAL ()

Additional Inherited Members
Protected Member Functions inherited from G4MagIntegratorStepper
void	SetIntegrationOrder (int order)
void	SetFSAL (bool flag=true)

Detailed Description

Definition at line 39 of file G4TsitourasRK45.hh.

Constructor & Destructor Documentation

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

Definition at line 51 of file G4TsitourasRK45.cc.

  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
    fLastStepLength(0.), fAuxStepper(0)
{
  const G4int numberOfVariables = noIntegrationVariables;
  G4cout << "G4TsitourasRK45 constructor called." << G4endl;
  
  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] ;


  // Must ensure space extra 'state' variables exists - i.e. yIn[7]
  const G4int numStateMax  = std::max(GetNumberOfStateVariables(), 8);  
  const G4int numStateVars = std::max(noIntegrationVariables,
                                      numStateMax );
                                   // GetNumberOfStateVariables() ); 
                                      
  yTemp = new G4double[numStateVars] ;
  yIn = new G4double[numStateVars] ;

  fLastInitialVector = new G4double[numberOfVariables] ;
  fLastFinalVector = new G4double[numberOfVariables] ;

  fLastDyDx = new G4double[numberOfVariables];

  fMidVector = new G4double[numberOfVariables];
  fMidError =  new G4double[numberOfVariables];
  if( primary )
  { 
    fAuxStepper = new G4TsitourasRK45(EqRhs, numberOfVariables, !primary);
  }
}

Here is the call graph for this function:

G4TsitourasRK45::~G4TsitourasRK45

(

)

Definition at line 95 of file G4TsitourasRK45.cc.

{
  delete[] ak2;
  delete[] ak3;
  delete[] ak4;
  delete[] ak5;
  delete[] ak6;
  delete[] ak7;
  delete[] ak8;

  delete[] yTemp;
  delete[] yIn;

  delete[] fLastInitialVector;
  delete[] fLastFinalVector;
  delete[] fLastDyDx;
  delete[] fMidVector;
  delete[] fMidError; 

  delete fAuxStepper;
}

Member Function Documentation

G4double G4TsitourasRK45::DistChord ( ) const

virtual

Implements G4MagIntegratorStepper.

Definition at line 309 of file G4TsitourasRK45.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]); 

  // Do half a step using StepNoErr

  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength, 
           fMidVector,   fMidError );

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

Here is the call graph for this function:

G4int G4TsitourasRK45::IntegratorOrder ( ) const

inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 72 of file G4TsitourasRK45.hh.

{return 4; }

void G4TsitourasRK45::Interpolate	(	const G4double	yInput[],
		const G4double	dydx[],
		const G4double	Step,
		G4double	yOut[],
		G4double	tau
	)

Definition at line 266 of file G4TsitourasRK45.cc.

                                                                                                                                {

    
    G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
    // Coefficients for all the seven stages.
    
    const G4int numberOfVariables= this->GetNumberOfVariables();
    
    G4double tau0 = tau;
    
    for(int i=0;i<numberOfVariables;i++)
    {
        yIn[i]=yInput[i];
    }
    
    G4double
    tau_2 = tau0*tau0 ;
//    tau_3 = tau0*tau_2,
//    tau_4 = tau_2*tau_2;

    bf1 = -1.0530884977290216*tau*(tau - 1.3299890189751412)*(tau_2 -
                1.4364028541716351*tau + 0.7139816917074209),
    bf2 = 0.1017*tau_2*(tau_2 - 2.1966568338249754*tau +
                        1.2949852507374631),
    bf3 = 2.490627285651252793*tau_2*(tau_2 - 2.38535645472061657*tau
                                      + 1.57803468208092486) ,
    bf4 = -16.54810288924490272*(tau - 1.21712927295533244)*
                    (tau - 0.61620406037800089)*tau_2,
    bf5 = 47.37952196281928122*(tau - 1.203071208372362603)*
                    (tau - 0.658047292653547382)*tau_2,
    bf6 = -34.87065786149660974*(tau - 1.2)*(tau -
                                0.666666666666666667)*tau_2,
    bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2;
    
    //Putting together the coefficients calculated as the respective stage coefficients
    for( int i=0; i<numberOfVariables; i++){
        yOut[i] = yIn[i] + Step*(  bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] + bf4*ak4[i]
                                 + bf5*ak5[i]  + bf6*ak6[i] + bf7*ak7[i]  ) ;
    }
}

Here is the call graph for this function:

void G4TsitourasRK45::interpolate	(	const G4double	yInput[],
		const G4double	dydx[],
		G4double	yOut[],
		G4double	Step,
		G4double	tau
	)

void G4TsitourasRK45::SetupInterpolation

(

)

Definition at line 260 of file G4TsitourasRK45.cc.

{
    //Nothing to be done
}

void G4TsitourasRK45::Stepper ( const G4double  y[], const G4double  dydx[], G4double  h, G4double  yout[], G4double  yerr[]  )

virtual

Implements G4MagIntegratorStepper.

Definition at line 129 of file G4TsitourasRK45.cc.

{
    G4int i;

    const G4double
    b21 = 0.161 ,
    
    b31 = -0.00848065549235698854 ,
    b32 = 0.335480655492356989 ,
    
    b41 =  2.89715305710549343 ,
    b42 = -6.35944848997507484 ,
    b43 = 4.36229543286958141 ,

    b51 = 5.325864828439257,
    b52 = -11.748883564062828,
    b53 = 7.49553934288983621 ,
    b54 = -0.09249506636175525,

    b61 = 5.8614554429464200,
    b62 = -12.9209693178471093 ,
    b63 = 8.1593678985761586 ,
    b64 = -0.071584973281400997,
    b65 = -0.0282690503940683829,

    
    b71 = 0.0964607668180652295 ,
    b72 = 0.01, 
    b73 = 0.479889650414499575,
    b74 = 1.37900857410374189,
    b75 = -3.2900695154360807,
    b76 = 2.32471052409977398,
    
//    c1 = 0.001780011052226 ,
//    c2 = 0.000816434459657 ,
//    c3 = -0.007880878010262 ,
//    c4 = 0.144711007173263 ,
//    c5 = -0.582357165452555 ,
//    c6 = 0.458082105929187 ,
//    c7 = 1.0/66.0 ;

    dc1 = 0.0935237485818927066 - b71 , // - 0.001780011052226,
    dc2 = 0.00865288314156636761 - b72, // - 0.000816434459657,
    dc3 = 0.492893099131431868 - b73 ,  // + 0.007880878010262,
    dc4 = 1.14023541226785810 - b74 ,   //   0.144711007173263,
    dc5 = - 2.3291801924393646 - b75,   // + 0.582357165452555,
    dc6 = 1.56887504931661552 - b76 ,   // - 0.458082105929187,
    dc7 = 0.025; //- 1.0/66.0 ;
    
//    dc1 = -3.0/1280.0,
//    dc2 = 0.0,
//    dc3 = 6561.0/632320.0,
//    dc4 = -343.0/20800.0,
//    dc5 = 243.0/12800.0,
//    dc6 = -1.0/95.0,
//    dc7 = 0.0   ;
    
    const G4int numberOfVariables= this->GetNumberOfVariables();
    
    // The number of variables to be integrated over
    yOut[7] = yTemp[7]  = yIn[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 Stage
    
    for(i=0;i<numberOfVariables;i++)
    {
        yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
    }
    RightHandSide(yTemp, ak3) ;              // 3rd Stage
    
    for(i=0;i<numberOfVariables;i++)
    {
        yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
    }
    RightHandSide(yTemp, ak4) ;              // 4th Stage
    
    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 Stage
    
    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 Stage
    
    for(i=0;i<numberOfVariables;i++)
    {
        yOut[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
                                 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
    }
    RightHandSide(yOut, ak7);       //7th Stage

    //Calculate the error in the step:
    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] ) ;
        
        // 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;
    
    return ;
}

Here is the call graph for this function:

Here is the caller graph for this function:

---

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

• source/geant4.10.03.p03/source/geometry/magneticfield/include/G4TsitourasRK45.hh
• source/geant4.10.03.p03/source/geometry/magneticfield/src/G4TsitourasRK45.cc