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