Geant4  10.01.p01
G4ConstRK4.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4ConstRK4.cc 66356 2012-12-18 09:02:32Z gcosmo $
28 //
29 //
30 // - 18.09.2008 - J.Apostolakis, T.Nikitina - Created
31 // -------------------------------------------------------------------
32 
33 #include "G4ConstRK4.hh"
34 #include "G4ThreeVector.hh"
35 #include "G4LineSection.hh"
36 
38 //
39 // Constructor sets the number of *State* variables (default = 8)
40 // The number of variables integrated is always 6
41 
42 G4ConstRK4::G4ConstRK4(G4Mag_EqRhs* EqRhs, G4int numStateVariables)
43  : G4MagErrorStepper(EqRhs, 6, numStateVariables)
44 {
45  // const G4int numberOfVariables= 6;
46  if( numStateVariables < 8 )
47  {
48  std::ostringstream message;
49  message << "The number of State variables at least 8 " << G4endl
50  << "Instead it is - numStateVariables= " << numStateVariables;
51  G4Exception("G4ConstRK4::G4ConstRK4()", "GeomField0002",
52  FatalException, message, "Use another Stepper!");
53  }
54 
55  fEq = EqRhs;
56  yMiddle = new G4double[8];
57  dydxMid = new G4double[8];
58  yInitial = new G4double[8];
59  yOneStep = new G4double[8];
60 
61  dydxm = new G4double[8];
62  dydxt = new G4double[8];
63  yt = new G4double[8];
64  Field[0]=0.; Field[1]=0.; Field[2]=0.;
65 }
66 
68 //
69 // Destructor
70 
72 {
73  delete [] yMiddle;
74  delete [] dydxMid;
75  delete [] yInitial;
76  delete [] yOneStep;
77  delete [] dydxm;
78  delete [] dydxt;
79  delete [] yt;
80 }
81 
83 //
84 // Given values for the variables y[0,..,n-1] and their derivatives
85 // dydx[0,...,n-1] known at x, use the classical 4th Runge-Kutta
86 // method to advance the solution over an interval h and return the
87 // incremented variables as yout[0,...,n-1], which is not a distinct
88 // array from y. The user supplies the routine RightHandSide(x,y,dydx),
89 // which returns derivatives dydx at x. The source is routine rk4 from
90 // NRC p. 712-713 .
91 
93  const G4double dydx[],
94  G4double h,
95  G4double yOut[])
96 {
97  G4double hh = h*0.5 , h6 = h/6.0 ;
98 
99  // 1st Step K1=h*dydx
100  yt[5] = yIn[5] + hh*dydx[5] ;
101  yt[4] = yIn[4] + hh*dydx[4] ;
102  yt[3] = yIn[3] + hh*dydx[3] ;
103  yt[2] = yIn[2] + hh*dydx[2] ;
104  yt[1] = yIn[1] + hh*dydx[1] ;
105  yt[0] = yIn[0] + hh*dydx[0] ;
107 
108  // 2nd Step K2=h*dydxt
109  yt[5] = yIn[5] + hh*dydxt[5] ;
110  yt[4] = yIn[4] + hh*dydxt[4] ;
111  yt[3] = yIn[3] + hh*dydxt[3] ;
112  yt[2] = yIn[2] + hh*dydxt[2] ;
113  yt[1] = yIn[1] + hh*dydxt[1] ;
114  yt[0] = yIn[0] + hh*dydxt[0] ;
116 
117  // 3rd Step K3=h*dydxm
118  // now dydxm=(K2+K3)/h
119  yt[5] = yIn[5] + h*dydxm[5] ;
120  dydxm[5] += dydxt[5] ;
121  yt[4] = yIn[4] + h*dydxm[4] ;
122  dydxm[4] += dydxt[4] ;
123  yt[3] = yIn[3] + h*dydxm[3] ;
124  dydxm[3] += dydxt[3] ;
125  yt[2] = yIn[2] + h*dydxm[2] ;
126  dydxm[2] += dydxt[2] ;
127  yt[1] = yIn[1] + h*dydxm[1] ;
128  dydxm[1] += dydxt[1] ;
129  yt[0] = yIn[0] + h*dydxm[0] ;
130  dydxm[0] += dydxt[0] ;
131  RightHandSideConst(yt,dydxt) ;
132 
133  // 4th Step K4=h*dydxt
134  yOut[5] = yIn[5]+h6*(dydx[5]+dydxt[5]+2.0*dydxm[5]);
135  yOut[4] = yIn[4]+h6*(dydx[4]+dydxt[4]+2.0*dydxm[4]);
136  yOut[3] = yIn[3]+h6*(dydx[3]+dydxt[3]+2.0*dydxm[3]);
137  yOut[2] = yIn[2]+h6*(dydx[2]+dydxt[2]+2.0*dydxm[2]);
138  yOut[1] = yIn[1]+h6*(dydx[1]+dydxt[1]+2.0*dydxm[1]);
139  yOut[0] = yIn[0]+h6*(dydx[0]+dydxt[0]+2.0*dydxm[0]);
140 
141 } // end of DumbStepper ....................................................
142 
144 //
145 // Stepper
146 
147 void
149  const G4double dydx[],
150  G4double hstep,
151  G4double yOutput[],
152  G4double yError [] )
153 {
154  const G4int nvar = 6; // number of variables integrated
155  const G4int maxvar= GetNumberOfStateVariables();
156 
157  // Correction for Richardson extrapolation
158  G4double correction = 1. / ( (1 << IntegratorOrder()) -1 );
159 
160  G4int i;
161 
162  // Saving yInput because yInput and yOutput can be aliases for same array
163  for (i=0; i<maxvar; i++) { yInitial[i]= yInput[i]; }
164 
165  // Must copy the part of the state *not* integrated to the output
166  for (i=nvar; i<maxvar; i++) { yOutput[i]= yInput[i]; }
167 
168  // yInitial[7]= yInput[7]; // The time is typically needed
169  yMiddle[7] = yInput[7]; // Copy the time from initial value
170  yOneStep[7] = yInput[7]; // As it contributes to final value of yOutput ?
171  // yOutput[7] = yInput[7]; // -> dumb stepper does it too for RK4
172  yError[7] = 0.0;
173 
174  G4double halfStep = hstep * 0.5;
175 
176  // Do two half steps
177  //
179  DumbStepper (yInitial, dydx, halfStep, yMiddle);
181  DumbStepper (yMiddle, dydxMid, halfStep, yOutput);
182 
183  // Store midpoint, chord calculation
184  //
186 
187  // Do a full Step
188  //
189  DumbStepper(yInitial, dydx, hstep, yOneStep);
190  for(i=0;i<nvar;i++)
191  {
192  yError [i] = yOutput[i] - yOneStep[i] ;
193  yOutput[i] += yError[i]*correction ;
194  // Provides accuracy increased by 1 order via the
195  // Richardson extrapolation
196  }
197 
199  fFinalPoint = G4ThreeVector( yOutput[0], yOutput[1], yOutput[2]);
200 
201  return;
202 }
203 
205 //
206 // Estimate the maximum distance from the curve to the chord
207 //
208 // We estimate this using the distance of the midpoint to chord.
209 // The method below is good only for angle deviations < 2 pi;
210 // this restriction should not be a problem for the Runge Kutta methods,
211 // which generally cannot integrate accurately for large angle deviations
212 
214 {
215  G4double distLine, distChord;
216 
217  if (fInitialPoint != fFinalPoint)
218  {
220  // This is a class method that gives distance of Mid
221  // from the Chord between the Initial and Final points
222  distChord = distLine;
223  }
224  else
225  {
226  distChord = (fMidPoint-fInitialPoint).mag();
227  }
228  return distChord;
229 }
G4ConstRK4(G4Mag_EqRhs *EquationMotion, G4int numberOfStateVariables=8)
Definition: G4ConstRK4.cc:42
CLHEP::Hep3Vector G4ThreeVector
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
G4double DistChord() const
Definition: G4ConstRK4.cc:213
G4double * yInitial
Definition: G4ConstRK4.hh:89
G4double Field[3]
Definition: G4ConstRK4.hh:91
G4double * dydxm
Definition: G4ConstRK4.hh:88
int G4int
Definition: G4Types.hh:78
G4ThreeVector fFinalPoint
Definition: G4ConstRK4.hh:86
G4ThreeVector fInitialPoint
Definition: G4ConstRK4.hh:86
G4double * dydxMid
Definition: G4ConstRK4.hh:89
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
Definition: G4ConstRK4.cc:148
G4double * yMiddle
Definition: G4ConstRK4.hh:89
G4Mag_EqRhs * fEq
Definition: G4ConstRK4.hh:90
G4ThreeVector fMidPoint
Definition: G4ConstRK4.hh:86
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double * yt
Definition: G4ConstRK4.hh:88
void GetConstField(const G4double y[], G4double Field[])
Definition: G4ConstRK4.hh:114
void RightHandSideConst(const G4double y[], G4double dydx[]) const
Definition: G4ConstRK4.hh:96
G4int GetNumberOfStateVariables() const
G4int IntegratorOrder() const
Definition: G4ConstRK4.hh:76
G4double * yOneStep
Definition: G4ConstRK4.hh:89
G4double * dydxt
Definition: G4ConstRK4.hh:88
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void DumbStepper(const G4double yIn[], const G4double dydx[], G4double h, G4double yOut[])
Definition: G4ConstRK4.cc:92