Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DoLoMcPriRK34.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 // Dormand-Lockyer-McGorrigan-Prince-6-3-4 non-FSAL implementation
27 // RK4(3)6FD - forced Non-FSAL
28 //
29 // Design/implementation by Somnath Banerjee
30 // Sponsored by Google in Google Summer of Code 2015.
31 // Supervision / code review: John Apostolakis
32 //
33 // First version: 7 July 2015
34 //
35 // G4DoLoMcPriRK34.cc
36 // Geant4
37 //
38 // History
39 // -----------------------------
40 // Created by Somnath on 7 July 2015
41 //
42 // This is the source file of G4DoLoMcPriRK34 class containing the
43 // definition of the Stepper() method that evaluates one Step in
44 // field propagation.
45 //
46 // The Butcher table of the Dormand-Lockyer-McGorrigan-Prince-6-3-4 method is as follows :
47 // [ to be added here ]
48 
49 #include "G4DoLoMcPriRK34.hh"
50 #include "G4LineSection.hh"
51 
52 // Constructor
54  G4int noIntegrationVariables,
55  G4bool primary)
56  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
57  fAuxStepper(0)
58 {
59  const G4int numberOfVariables = noIntegrationVariables;
60 
61  //New Chunk of memory being created for use by the Stepper
62 
63  //aki - for storing intermediate RHS
64  ak2 = new G4double[numberOfVariables];
65  ak3 = new G4double[numberOfVariables];
66  ak4 = new G4double[numberOfVariables];
67  ak5 = new G4double[numberOfVariables];
68  ak6 = new G4double[numberOfVariables];
69 
70  yTemp = new G4double[numberOfVariables] ;
71  yIn = new G4double[numberOfVariables] ;
72 
73  fLastInitialVector = new G4double[numberOfVariables] ;
74  fLastFinalVector = new G4double[numberOfVariables] ;
75  fLastDyDx = new G4double[numberOfVariables];
76 
77  fMidVector = new G4double[numberOfVariables];
78  fMidError = new G4double[numberOfVariables];
79  if( primary )
80  {
81  fAuxStepper = new G4DoLoMcPriRK34(EqRhs, numberOfVariables,
82  !primary);
83  }
84 }
85 
86 
87 //Destructor
89 {
90  //clear all previously allocated memory for Stepper and DistChord
91  delete[] ak2;
92  delete[] ak3;
93  delete[] ak4;
94  delete[] ak5;
95  delete[] ak6;
96 
97  delete[] yTemp;
98  delete[] yIn;
99 
100  delete[] fLastInitialVector;
101  delete[] fLastFinalVector;
102  delete[] fLastDyDx;
103  delete[] fMidVector;
104  delete[] fMidError;
105 
106  delete fAuxStepper;
107 
108 
109 }
110 
111 
112 //Stepper :
113 
114 // Passing in the value of yInput[],the first time dydx[] and Step length
115 // Giving back yOut and yErr arrays for output and error respectively
116 
117 void G4DoLoMcPriRK34::Stepper(const G4double yInput[],
118  const G4double DyDx[],
119  G4double Step,
120  G4double yOut[],
121  G4double yErr[] )
122 {
123  G4int i;
124 
125  //The various constants defined on the basis of butcher tableu
126  const G4double //G4double - only once
127 
128  b21 = 7.0/27.0 ,
129 
130 
131  b31 = 7.0/72.0 ,
132  b32 = 7.0/24.0 ,
133 
134  b41 = 3043.0/3528.0 ,
135  b42 = -3757.0/1176.0 ,
136  b43 = 1445.0/441.0,
137 
138  b51 = 17617.0/11662.0 ,
139  b52 = -4023.0/686.0 ,
140  b53 = 9372.0/1715.0 ,
141  b54 = -66.0/595.0 ,
142 
143  b61 = 29.0/238.0 ,
144  b62 = 0.0 ,
145  b63 = 216.0/385.0 ,
146  b64 = 54.0/85.0 ,
147  b65 = -7.0/22.0 ,
148 
149 
150 
151  dc1 = 363.0/2975.0 - b61 ,
152  dc2 = 0.0 - b62 ,
153  dc3 = 981.0/1750.0 - b63,
154  dc4 = 2709.0/4250.0 - b64 ,
155  dc5 = -3.0/10.0 - b65 ,
156  dc6 = -1.0/50.0 ; //end of declaration
157 
158 
159  const G4int numberOfVariables= this->GetNumberOfVariables();
160 
161  // The number of variables to be integrated over
162  yOut[7] = yTemp[7] = yIn[7];
163  // Saving yInput because yInput and yOut can be aliases for same array
164 
165  for(i=0;i<numberOfVariables;i++)
166  {
167  yIn[i]=yInput[i];
168  }
169 
170 
171 
172  // RightHandSide(yIn, DyDx) ;
173  // 1st stage - Not doing, getting passed
174 
175  for(i=0;i<numberOfVariables;i++)
176  {
177  yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
178  }
179  RightHandSide(yTemp, ak2) ; // 2nd stage
180 
181  for(i=0;i<numberOfVariables;i++)
182  {
183  yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
184  }
185  RightHandSide(yTemp, ak3) ; // 3rd stage
186 
187  for(i=0;i<numberOfVariables;i++)
188  {
189  yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
190  }
191  RightHandSide(yTemp, ak4) ; // 4th stage
192 
193  for(i=0;i<numberOfVariables;i++)
194  {
195  yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
196  b54*ak4[i]) ;
197  }
198  RightHandSide(yTemp, ak5) ; // 5th stage
199 
200  for(i=0;i<numberOfVariables;i++)
201  {
202  yOut[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
203  b64*ak4[i] + b65*ak5[i]) ;
204  }
205  RightHandSide(yOut, ak6) ; // 6th and Final stage
206 
207 
208 
209  for(i=0;i<numberOfVariables;i++)
210  {
211 
212  yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
213  dc5*ak5[i] + dc6*ak6[i] ) ;
214 
215 
216  // Store Input and Final values, for possible use in calculating chord
217  fLastInitialVector[i] = yIn[i] ;
218  fLastFinalVector[i] = yOut[i];
219  fLastDyDx[i] = DyDx[i];
220 
221 
222  }
223 
224  fLastStepLength = Step;
225 
226  return ;
227 }
228 
229 
230 //The following has not been tested
231 
232 //The DistChord() function fot the class - must define it here.
234 {
235  G4double distLine, distChord;
236  G4ThreeVector initialPoint, finalPoint, midPoint;
237 
238  // Store last initial and final points (they will be overwritten in self-Stepper call!)
239  initialPoint = G4ThreeVector( fLastInitialVector[0],
240  fLastInitialVector[1], fLastInitialVector[2]);
241  finalPoint = G4ThreeVector( fLastFinalVector[0],
242  fLastFinalVector[1], fLastFinalVector[2]);
243 
244  // Do half a Step using StepNoErr
245 
246  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
247  fMidVector, fMidError );
248 
249  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
250 
251  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
252  // distance of Chord
253 
254 
255  if (initialPoint != finalPoint)
256  {
257  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
258  distChord = distLine;
259  }
260  else
261  {
262  distChord = (midPoint-initialPoint).mag();
263  }
264  return distChord;
265 }
266 
268 {}
269 
270 void G4DoLoMcPriRK34::SetupInterpolate( const G4double /* yInput */ [] ,
271  const G4double /* dydx */ [] ,
272  const G4double /* Step */ )
273 {
274  //Do Nothing
275 }
276 
277 
279  G4double yOut[])
280 {
281  Interpolate( fLastInitialVector, fLastDyDx, fLastStepLength, yOut, tau );
282 }
283 
284 // Function to evaluate the interpolation at tau fraction of the step
286  const G4double dydx[],
287  const G4double Step,
288  G4double yOut[],
289  G4double tau ){
290  G4double
291  bf1, bf2, bf3, bf4, bf5, bf6;
292 
293 
294  const G4int numberOfVariables= this->GetNumberOfVariables();
295 
296  for(int i=0;i<numberOfVariables;i++)
297  {
298  yIn[i]=yInput[i];
299  }
300 
301  G4double
302  tau_2 = tau*tau ,
303  tau_3 = tau*tau_2;
304 
305  //Calculating the polynomials (coefficients for the respective stages)
306  bf1 = -(162.0*tau_3 - 504.0*tau_2 + 551.0*tau - 238.0)/238.0 ,
307  bf2 = 0.0 ,
308  bf3 = 27.0*tau*(27.0*tau_2 - 70.0*tau + 51.0 )/385.0 ,
309  bf4 = -27*tau*(27.0*tau_2 - 50.0*tau + 21.0)/85.0 ,
310  bf5 = 7.0*tau*(2232.0*tau_2 - 4166.0*tau + 1785.0 )/3278.0 ,
311  bf6 = tau*(tau - 1.0)*(387.0*tau - 238.0)/149.0 ;
312 
313  for( int i=0; i<numberOfVariables; i++){
314  yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] +
315  bf4*ak4[i] + bf5*ak5[i] + bf6*ak6[i] ) ;
316  }
317 
318 
319 
320 }
321 
322 
323 //-------Verified------- - hackabot
324 
void SetupInterpolate(const G4double yInput[], const G4double dydx[], const G4double Step)
CLHEP::Hep3Vector G4ThreeVector
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
G4double DistChord() const
void Interpolate(const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
int G4int
Definition: G4Types.hh:78
G4int GetNumberOfVariables() const
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
bool G4bool
Definition: G4Types.hh:79
Definition: Step.hh:41
void RightHandSide(const double y[], double dydx[])
double G4double
Definition: G4Types.hh:76
G4DoLoMcPriRK34(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)