G4ConstRK4 Class Reference

#include <G4ConstRK4.hh>

Inheritance diagram for G4ConstRK4:

G4MagErrorStepper G4MagIntegratorStepper

Public Member Functions

 G4ConstRK4 (G4Mag_EqRhs *EquationMotion, G4int numberOfStateVariables=8)
 ~G4ConstRK4 ()
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
void DumbStepper (const G4double yIn[], const G4double dydx[], G4double h, G4double yOut[])
G4double DistChord () const
void RightHandSideConst (const G4double y[], G4double dydx[]) const
void GetConstField (const G4double y[], G4double Field[])
G4int IntegratorOrder () const

Detailed Description

Definition at line 51 of file G4ConstRK4.hh.

Constructor & Destructor Documentation

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

Definition at line 42 of file G4ConstRK4.cc.

References FatalException, G4endl, and G4Exception().

00043   : G4MagErrorStepper(EqRhs, 6, numStateVariables)
00044 {
00045   // const G4int numberOfVariables= 6;
00046   if( numStateVariables < 8 ) 
00047   {
00048     std::ostringstream message;
00049     message << "The number of State variables at least 8 " << G4endl
00050             << "Instead it is - numStateVariables= " << numStateVariables;
00051     G4Exception("G4ConstRK4::G4ConstRK4()", "GeomField0002",
00052                 FatalException, message, "Use another Stepper!");
00053   }
00054 
00055   fEq = EqRhs;
00056   yMiddle  = new G4double[8];
00057   dydxMid  = new G4double[8];
00058   yInitial = new G4double[8];
00059   yOneStep = new G4double[8];
00060 
00061   dydxm = new G4double[8];
00062   dydxt = new G4double[8]; 
00063   yt    = new G4double[8]; 
00064   Field[0]=0.; Field[1]=0.; Field[2]=0.;
00065 }

G4ConstRK4::~G4ConstRK4 (  ) 

Definition at line 71 of file G4ConstRK4.cc.

00072 {
00073    delete [] yMiddle;
00074    delete [] dydxMid;
00075    delete [] yInitial;
00076    delete [] yOneStep;
00077    delete [] dydxm;
00078    delete [] dydxt;
00079    delete [] yt;
00080 }

Member Function Documentation

G4double G4ConstRK4::DistChord (  )  const [virtual]

Reimplemented from G4MagErrorStepper.

Definition at line 213 of file G4ConstRK4.cc.

References G4LineSection::Distline().

00214 {
00215   G4double distLine, distChord; 
00216 
00217   if (fInitialPoint != fFinalPoint)
00218   {
00219      distLine= G4LineSection::Distline( fMidPoint, fInitialPoint, fFinalPoint );
00220        // This is a class method that gives distance of Mid 
00221        // from the Chord between the Initial and Final points
00222      distChord = distLine;
00223   }
00224   else
00225   {
00226      distChord = (fMidPoint-fInitialPoint).mag();
00227   }
00228   return distChord;
00229 }

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

Implements G4MagErrorStepper.

Definition at line 92 of file G4ConstRK4.cc.

References RightHandSideConst().

Referenced by Stepper().

00096 {
00097    G4double  hh = h*0.5 , h6 = h/6.0  ;
00098    
00099    // 1st Step K1=h*dydx
00100    yt[5] = yIn[5] + hh*dydx[5] ;
00101    yt[4] = yIn[4] + hh*dydx[4] ;
00102    yt[3] = yIn[3] + hh*dydx[3] ;
00103    yt[2] = yIn[2] + hh*dydx[2] ;
00104    yt[1] = yIn[1] + hh*dydx[1] ;
00105    yt[0] = yIn[0] + hh*dydx[0] ;
00106    RightHandSideConst(yt,dydxt) ;        
00107 
00108    // 2nd Step K2=h*dydxt
00109    yt[5] = yIn[5] + hh*dydxt[5] ;
00110    yt[4] = yIn[4] + hh*dydxt[4] ;
00111    yt[3] = yIn[3] + hh*dydxt[3] ;
00112    yt[2] = yIn[2] + hh*dydxt[2] ;
00113    yt[1] = yIn[1] + hh*dydxt[1] ;
00114    yt[0] = yIn[0] + hh*dydxt[0] ;
00115    RightHandSideConst(yt,dydxm) ;     
00116 
00117    // 3rd Step K3=h*dydxm
00118    // now dydxm=(K2+K3)/h
00119    yt[5] = yIn[5] + h*dydxm[5] ;
00120    dydxm[5] += dydxt[5] ;  
00121    yt[4] = yIn[4] + h*dydxm[4] ;
00122    dydxm[4] += dydxt[4] ;  
00123    yt[3] = yIn[3] + h*dydxm[3] ;
00124    dydxm[3] += dydxt[3] ;  
00125    yt[2] = yIn[2] + h*dydxm[2] ;
00126    dydxm[2] += dydxt[2] ;  
00127    yt[1] = yIn[1] + h*dydxm[1] ;
00128    dydxm[1] += dydxt[1] ;  
00129    yt[0] = yIn[0] + h*dydxm[0] ;
00130    dydxm[0] += dydxt[0] ;  
00131    RightHandSideConst(yt,dydxt) ;   
00132 
00133    // 4th Step K4=h*dydxt
00134    yOut[5] = yIn[5]+h6*(dydx[5]+dydxt[5]+2.0*dydxm[5]);
00135    yOut[4] = yIn[4]+h6*(dydx[4]+dydxt[4]+2.0*dydxm[4]);
00136    yOut[3] = yIn[3]+h6*(dydx[3]+dydxt[3]+2.0*dydxm[3]);
00137    yOut[2] = yIn[2]+h6*(dydx[2]+dydxt[2]+2.0*dydxm[2]);
00138    yOut[1] = yIn[1]+h6*(dydx[1]+dydxt[1]+2.0*dydxm[1]);
00139    yOut[0] = yIn[0]+h6*(dydx[0]+dydxt[0]+2.0*dydxm[0]);
00140    
00141 }  // end of DumbStepper ....................................................

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

Definition at line 114 of file G4ConstRK4.hh.

Referenced by Stepper().

00115 {
00116   G4double  PositionAndTime[4];
00117 
00118   PositionAndTime[0] = y[0];
00119   PositionAndTime[1] = y[1];
00120   PositionAndTime[2] = y[2];
00121   // Global Time
00122   PositionAndTime[3] = y[7];  
00123   fEq -> GetFieldValue(PositionAndTime, B) ;
00124 }

G4int G4ConstRK4::IntegratorOrder (  )  const [inline, virtual]

Implements G4MagIntegratorStepper.

Definition at line 76 of file G4ConstRK4.hh.

Referenced by Stepper().

00076 { return 4; }

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

Definition at line 96 of file G4ConstRK4.hh.

References G4Mag_EqRhs::FCof().

Referenced by DumbStepper(), and Stepper().

00098 {
00099   
00100   G4double momentum_mag_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5];
00101   G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square );
00102     
00103   G4double cof =fEq->FCof()*inv_momentum_magnitude;
00104 
00105   dydx[0] = y[3]*inv_momentum_magnitude;       //  (d/ds)x = Vx/V
00106   dydx[1] = y[4]*inv_momentum_magnitude;       //  (d/ds)y = Vy/V
00107   dydx[2] = y[5]*inv_momentum_magnitude;       //  (d/ds)z = Vz/V
00108  
00109   dydx[3] = cof*(y[4]*Field[2] - y[5]*Field[1]) ;   // Ax = a*(Vy*Bz - Vz*By)
00110   dydx[4] = cof*(y[5]*Field[0] - y[3]*Field[2]) ;   // Ay = a*(Vz*Bx - Vx*Bz)
00111   dydx[5] = cof*(y[3]*Field[1] - y[4]*Field[0]) ;   // Az = a*(Vx*By - Vy*Bx)
00112 }

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

Reimplemented from G4MagErrorStepper.

Definition at line 148 of file G4ConstRK4.cc.

References DumbStepper(), GetConstField(), G4MagIntegratorStepper::GetNumberOfStateVariables(), IntegratorOrder(), and RightHandSideConst().

00153 {
00154    const G4int nvar = 6;  // number of variables integrated
00155    const G4int maxvar= GetNumberOfStateVariables();
00156 
00157    // Correction for Richardson extrapolation
00158    G4double  correction = 1. / ( (1 << IntegratorOrder()) -1 );
00159 
00160    G4int i;
00161    
00162    // Saving yInput because yInput and yOutput can be aliases for same array
00163    for (i=0;    i<maxvar; i++) { yInitial[i]= yInput[i]; }
00164  
00165    // Must copy the part of the state *not* integrated to the output
00166    for (i=nvar; i<maxvar; i++) { yOutput[i]=  yInput[i]; }
00167 
00168    // yInitial[7]= yInput[7];  //  The time is typically needed
00169    yMiddle[7]  = yInput[7];   // Copy the time from initial value 
00170    yOneStep[7] = yInput[7];   // As it contributes to final value of yOutput ?
00171    // yOutput[7] = yInput[7];  // -> dumb stepper does it too for RK4
00172    yError[7] = 0.0;         
00173 
00174    G4double halfStep = hstep * 0.5; 
00175 
00176    // Do two half steps
00177    //
00178    GetConstField(yInitial,Field);
00179    DumbStepper  (yInitial,  dydx,   halfStep, yMiddle);
00180    RightHandSideConst(yMiddle, dydxMid);    
00181    DumbStepper  (yMiddle, dydxMid, halfStep, yOutput); 
00182 
00183    // Store midpoint, chord calculation
00184    //
00185    fMidPoint = G4ThreeVector( yMiddle[0],  yMiddle[1],  yMiddle[2]); 
00186 
00187    // Do a full Step
00188    //
00189    DumbStepper(yInitial, dydx, hstep, yOneStep);
00190    for(i=0;i<nvar;i++)
00191    {
00192       yError [i] = yOutput[i] - yOneStep[i] ;
00193       yOutput[i] += yError[i]*correction ;
00194         // Provides accuracy increased by 1 order via the 
00195         // Richardson extrapolation  
00196    }
00197 
00198    fInitialPoint = G4ThreeVector( yInitial[0], yInitial[1], yInitial[2]); 
00199    fFinalPoint   = G4ThreeVector( yOutput[0],  yOutput[1],  yOutput[2]); 
00200 
00201    return;
00202 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:42 2013 for Geant4 by  doxygen 1.4.7