Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FSALBogackiShampine45.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 // Bogacki-Shampine - 8 - 5(4) FSAL implementation by Somnath Banerjee
27 // Supervision / code review: John Apostolakis
28 //
29 // Sponsored by Google in Google Summer of Code 2015.
30 //
31 // First version: 26 May 2015
32 //
33 // History
34 // -----------------------------
35 // Created by Somnath on 26 May 2015
36 //
38 // Renamed to G4 standard naming
39 // Plan is that this source file / class will be merged with the updated
40 // BogackiShampine45 class, which contains improvements (May 2016)
41 // J. Apostolakis, 31 May 2016
43 //
44 //
45 //This is the source file of BogackiShampine45 class containing the
46 //definition of the stepper() method that evaluates one step in
47 //field propagation.
48 //The Butcher table of the Bogacki-Shampine-8-4-5 method is as follows :
49 //
50 //0 |
51 //1/6 | 1/6
52 //2/9 | 2/27 4/27
53 //3/7 | 183/1372 -162/343 1053/1372
54 //2/3 | 68/297 -4/11 42/143 1960/3861
55 //3/4 | 597/22528 81/352 63099/585728 58653/366080 4617/20480
56 //1 | 174197/959244 -30942/79937 8152137/19744439 666106/1039181 -29421/29068 482048/414219
57 //1 | 587/8064 0 4440339/15491840 24353/124800 387/44800 2152/5985 7267/94080
58 //-------------------------------------------------------------------------------------------------------------------
59 // 587/8064 0 4440339/15491840 24353/124800 387/44800 2152/5985 7267/94080 0
60 // 2479/34992 0 123/416 612941/3411720 43/1440 2272/6561 79937/1113912 3293/556956
61 
62 #include <cassert>
63 
64 #include "G4FSALBogackiShampine45.hh"
65 #include "G4LineSection.hh"
66 
67 //Constructor
68 G4FSALBogackiShampine45::G4FSALBogackiShampine45(G4EquationOfMotion *EqRhs,
69  G4int noIntegrationVariables,
70  G4bool primary)
71 : G4VFSALIntegrationStepper(EqRhs, noIntegrationVariables)
72 {
73 
74  const G4int numberOfVariables = noIntegrationVariables;
75 
76  //New Chunk of memory being created for use by the stepper
77 
78  //aki - for storing intermediate RHS
79  ak2 = new G4double[numberOfVariables];
80  ak3 = new G4double[numberOfVariables];
81  ak4 = new G4double[numberOfVariables];
82  ak5 = new G4double[numberOfVariables];
83  ak6 = new G4double[numberOfVariables];
84  ak7 = new G4double[numberOfVariables];
85  ak8 = new G4double[numberOfVariables];
86 
87  ak9 = new G4double[numberOfVariables];
88  ak10 = new G4double[numberOfVariables];
89  ak11 = new G4double[numberOfVariables];
90 
91  assert ( GetNumberOfStateVariables() >= 8 );
92  const G4int numStateVars = std::max(noIntegrationVariables,
93  GetNumberOfStateVariables() );
94 
95  // Must ensure space extra 'state' variables exists - i.e. yIn[7]
96  yTemp = new G4double[numStateVars];
97  yIn = new G4double[numStateVars] ;
98 
99  fLastInitialVector = new G4double[numStateVars] ;
100  fLastFinalVector = new G4double[numStateVars] ;
101  fLastDyDx = new G4double[numberOfVariables]; // Only derivatives
102 
103  fMidVector = new G4double[numStateVars];
104  fMidError = new G4double[numStateVars];
105 
106  pseudoDydx_for_DistChord = new G4double[numberOfVariables];
107 
108  fMidVector = new G4double[numberOfVariables];
109  fMidError = new G4double[numberOfVariables];
110  if( primary )
111  {
112  fAuxStepper = new G4FSALBogackiShampine45(EqRhs, numberOfVariables,
113  !primary);
114  }
115 }
116 
117 
118 //Destructor
119 G4FSALBogackiShampine45::~G4FSALBogackiShampine45(){
120  //clear all previously allocated memory for stepper and DistChord
121  delete[] ak2;
122  delete[] ak3;
123  delete[] ak4;
124  delete[] ak5;
125  delete[] ak6;
126  delete[] ak7;
127  delete[] ak8;
128 
129  delete[] ak9;
130  delete[] ak10;
131  delete[] ak11;
132 
133  delete[] yTemp;
134  delete[] yIn;
135 
136  delete[] fLastInitialVector;
137  delete[] fLastFinalVector;
138  delete[] fLastDyDx;
139  delete[] fMidVector;
140  delete[] fMidError;
141 
142  delete fAuxStepper;
143 
144  delete[] pseudoDydx_for_DistChord;
145 }
146 
147 
148 //Stepper :
149 
150 // Passing in the value of yInput[],the first time dydx[] and Step length
151 // Giving back yOut and yErr arrays for output and error respectively
152 
153 void G4FSALBogackiShampine45::Stepper(const G4double yInput[],
154  const G4double dydx[],
155  G4double Step,
156  G4double yOut[],
157  G4double yErr[],
158  G4double nextDydx[])
159 {
160  G4int i;
161 
162  //The various constants defined on the basis of butcher tableu
163  const G4double //G4double - only once
164 
165  b21 = 1.0/6.0 ,
166  b31 = 2.0/27.0 , b32 = 4.0/27.0,
167 
168  b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
169 
170  b51 = 68.0/297.0, b52 = -4.0/11.0,
171  b53 = 42.0/143.0, b54 = 1960.0/3861.0,
172 
173  b61 = 597.0/22528.0, b62 = 81.0/352.0,
174  b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
175  b65 = 4617.0/20480.0,
176 
177  b71 = 174197.0/959244.0, b72 = -30942.0/79937.0,
178  b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
179  b75 = -29421.0/29068.0, b76 = 482048.0/414219.0,
180 
181  b81 = 587.0/8064.0, b82 = 0.0,
182  b83 = 4440339.0/15491840.0, b84 = 24353.0/124800.0,
183  b85 = 387.0/44800.0, b86 = 2152.0/5985.0,
184  b87 = 7267.0/94080.0,
185 
186 
187 // c1 = 2479.0/34992.0,
188 // c2 = 0.0,
189 // c3 = 123.0/416.0,
190 // c4 = 612941.0/3411720.0,
191 // c5 = 43.0/1440.0,
192 // c6 = 2272.0/6561.0,
193 // c7 = 79937.0/1113912.0,
194 // c8 = 3293.0/556956.0,
195 
196  //For the embedded higher order method only the difference of values
197  // taken and is used directly later instead of defining the last row
198  // of butcher table in a separate set of variables and taking the
199  // difference there
200 
201 
202  dc1 = b81 - 2479.0/34992.0 ,
203  dc2 = 0.0,
204  dc3 = b83 - 123.0/416.0 ,
205  dc4 = b84 - 612941.0/3411720.0,
206  dc5 = b85 - 43.0/1440.0,
207  dc6 = b86 - 2272.0/6561.0,
208  dc7 = b87 - 79937.0/1113912.0,
209  dc8 = -3293.0/556956.0; //end of declaration
210 
211 
212  const G4int numberOfVariables= this->GetNumberOfVariables();
213 
214  G4double *DyDx = new G4double[numberOfVariables];
215 
216  // The number of variables to be integrated over
217  yOut[7] = yTemp[7] = yIn[7];
218  // Saving yInput because yInput and yOut can be aliases for same array
219 
220  for(i=0;i<numberOfVariables;i++)
221  {
222  yIn[i]=yInput[i];
223  DyDx[i] = dydx[i];
224  }
225 
226 
227  // RightHandSide(yIn, dydx) ;
228  // 1st Step - Not doing, getting passed
229 
230  for(i=0;i<numberOfVariables;i++)
231  {
232  yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
233  }
234  RightHandSide(yTemp, ak2) ; // 2nd Step
235 
236  for(i=0;i<numberOfVariables;i++)
237  {
238  yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
239  }
240  RightHandSide(yTemp, ak3) ; // 3rd Step
241 
242  for(i=0;i<numberOfVariables;i++)
243  {
244  yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
245  }
246  RightHandSide(yTemp, ak4) ; // 4th Step
247 
248  for(i=0;i<numberOfVariables;i++)
249  {
250  yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
251  b54*ak4[i]) ;
252  }
253  RightHandSide(yTemp, ak5) ; // 5th Step
254 
255  for(i=0;i<numberOfVariables;i++)
256  {
257  yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
258  b64*ak4[i] + b65*ak5[i]) ;
259  }
260  RightHandSide(yTemp, ak6) ; // 6th Step
261 
262  for(i=0;i<numberOfVariables;i++)
263  {
264  yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
265  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
266  }
267  RightHandSide(yTemp, ak7); //7th Step
268 
269  for(i=0;i<numberOfVariables;i++)
270  {
271  yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] +
272  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
273  b87*ak7[i]);
274  }
275  RightHandSide(yOut, ak8); //8th Step - Final one Using FSAL
276 
277 
278  for(i=0;i<numberOfVariables;i++)
279  {
280 
281  yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
282  dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
283 
284 
285  //FSAL stepper : Must pass the last DyDx for the next step, here ak8
286  nextDydx[i] = ak8[i];
287 
288  // Store Input and Final values, for possible use in calculating chord
289  fLastInitialVector[i] = yIn[i] ;
290  fLastFinalVector[i] = yOut[i];
291  fLastDyDx[i] = DyDx[i];
292 
293  }
294 
295  fLastStepLength = Step;
296 
297  return ;
298 }
299 //
300 //G4double* G4FSALBogackiShampine45::getLastDydx(){
301 // return ak8;
302 //}
303 
304 //The following has not been tested
305 
306 //The DistChord() function fot the class - must define it here.
307 G4double G4FSALBogackiShampine45::DistChord() const
308 {
309  G4double distLine, distChord;
310  G4ThreeVector initialPoint, finalPoint, midPoint;
311 
312 
313  // Store last initial and final points (they will be overwritten in self-Stepper call!)
314  initialPoint = G4ThreeVector( fLastInitialVector[0],
315  fLastInitialVector[1], fLastInitialVector[2]);
316  finalPoint = G4ThreeVector( fLastFinalVector[0],
317  fLastFinalVector[1], fLastFinalVector[2]);
318 
319  // Do half a step using StepNoErr
320 
321  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
322  fMidVector, fMidError, pseudoDydx_for_DistChord );
323 
324  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
325 
326  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
327  // distance of Chord
328 
329 
330  if (initialPoint != finalPoint)
331  {
332  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
333  distChord = distLine;
334  }
335  else
336  {
337  distChord = (midPoint-initialPoint).mag();
338  }
339  return distChord;
340 }
341 
342 
343 
344 
345 void G4FSALBogackiShampine45::interpolate( const G4double yInput[],
346  const G4double dydx[],
347  G4double yOut[],
348  G4double Step,
349  G4double tau
350  ){
351 
352 
353 // G4double *ak9, *ak10, *ak11;
354 
355  G4double
356  a91 = 455.0/6144.0 ,
357  a92 = 0.0 ,
358  a93 = 10256301.0/35409920.0 ,
359  a94 = 2307361.0/17971200.0 ,
360  a95 = -387.0/102400.0 ,
361  a96 = 73.0/5130.0 ,
362  a97 = -7267.0/215040.0 ,
363  a98 = 1.0/32.0 ,
364 
365 
366  a101 = -837888343715.0/13176988637184.0 ,
367  a102 = 30409415.0/52955362.0 ,
368  a103 = -48321525963.0/759168069632.0 ,
369  a104 = 8530738453321.0/197654829557760.0 ,
370  a105 = 1361640523001.0/1626788720640.0 ,
371  a106 = -13143060689.0/38604458898.0 ,
372  a107 = 18700221969.0/379584034816.0 ,
373  a108 = -5831595.0/847285792.0 ,
374  a109 = -5183640.0/26477681.0 ,
375 
376  a111 = 98719073263.0/1551965184000.0 ,
377  a112 = 1307.0/123552.0 ,
378  a113 = 4632066559387.0/70181753241600.0 ,
379  a114 = 7828594302389.0/382182512025600.0 ,
380  a115 = 40763687.0/11070259200.0 ,
381  a116 = 34872732407.0/224610586200.0 ,
382  a117 = -2561897.0/30105600.0 ,
383  a118 = 1.0/10.0 ,
384  a119 = -1.0/10.0 ,
385  a1110 = -1403317093.0/11371610250.0 ;
386 
387 // a91 = 455.e0/6144.e0 ,
388 // a101 = -837888343715.e0/13176988637184.e0 ,
389 // a111 = 98719073263.e0/1551965184000.e0 ,
390 // a92 = 0.e0 ,
391 // a102 = 30409415.e0/52955362.e0 ,
392 // a112 = 1307.e0/123552.e0 ,
393 // a93 = 10256301.e0/35409920.e0 ,
394 // a103 = -48321525963.e0/759168069632.e0 ,
395 // a113 = 4632066559387.e0/70181753241600.e0 ,
396 // a94 = 2307361.e0/17971200.e0 ,
397 // a104 = 8530738453321.e0/197654829557760.e0 ,
398 // a114 = 7828594302389.e0/382182512025600.e0 ,
399 // a95 = -387.e0/102400.e0 ,
400 // a105 = 1361640523001.e0/1626788720640.e0 ,
401 // a115 = 40763687.e0/11070259200.e0 ,
402 // a96 = 73.e0/5130.e0 ,
403 // a106 = -13143060689.e0/38604458898.e0 ,
404 // a116 = 34872732407.e0/224610586200.e0 ,
405 // a97 = -7267.e0/215040.e0 ,
406 // a107 = 18700221969.e0/379584034816.e0 ,
407 // a117 = -2561897.e0/30105600.e0 ,
408 // a98 = 1.e0/32.e0 ,
409 // a108 = -5831595.e0/847285792.e0 ,
410 // a118 = 1.e0/10.e0 ,
411 // a109 = -5183640.e0/26477681.e0 ,
412 // a119 = -1.e0/10.e0 ,
413 // a1110 = -1403317093.e0/11371610250.e0 ;
414 
415 
416 
417  // --------------------------------------------------------
418  // COEFFICIENTS FOR INTERPOLANT bi WITH 11 STAGES
419  // --------------------------------------------------------
420 
421  G4double bi[13][7], b[13]; //For bi[1][1] to bi[11][6] and b[1] to b[11]
422 
423  for(int i=1; i<= 11; i++)
424  bi[i][1] = 0.0 ;
425 
426  for(int i=1; i<=6; i++)
427  bi[2][i] = 0.0 ;
428 
429  bi[1][6] = -12134338393.0/1050809760.0 ,
430  bi[1][5] = -1620741229.0/50038560.0 ,
431  bi[1][4] = -2048058893.0/59875200.0 ,
432  bi[1][3] = -87098480009.0/5254048800.0 ,
433  bi[1][2] = -11513270273.0/3502699200.0 ,
434  //
435  bi[3][6] = -33197340367.0/1218433216.0 ,
436  bi[3][5] = -539868024987.0/6092166080.0 ,
437  bi[3][4] = -39991188681.0/374902528.0 ,
438  bi[3][3] = -69509738227.0/1218433216.0 ,
439  bi[3][2] = -29327744613.0/2436866432.0 ,
440  //
441  bi[4][6] = -284800997201.0/19905339168.0 ,
442  bi[4][5] = -7896875450471.0/165877826400.0 ,
443  bi[4][4] = -333945812879.0/5671036800.0 ,
444  bi[4][3] = -16209923456237.0/497633479200.0 ,
445  bi[4][2] = -2382590741699.0/331755652800.0 ,
446  //
447  bi[5][6] = -540919.0/741312.0 ,
448  bi[5][5] = -103626067.0/43243200.0 ,
449  bi[5][4] = -633779.0/211200.0 ,
450  bi[5][3] = -32406787.0/18532800.0 ,
451  bi[5][2] = -36591193.0/86486400.0 ,
452  //
453  bi[6][6] = 7157998304.0/374350977.0 ,
454  bi[6][5] = 30405842464.0/623918295.0 ,
455  bi[6][4] = 183022264.0/5332635.0 ,
456  bi[6][3] = -3357024032.0/1871754885.0 ,
457  bi[6][2] = -611586736.0/89131185.0 ,
458  //
459  bi[7][6] = -138073.0/9408.0 ,
460  bi[7][5] = -719433.0/15680.0 ,
461  bi[7][4] = -1620541.0/31360.0 ,
462  bi[7][3] = -385151.0/15680.0 ,
463  bi[7][2] = -65403.0/15680.0 ,
464  //
465  bi[8][6] = 1245.0/64.0 ,
466  bi[8][5] = 3991.0/64.0 ,
467  bi[8][4] = 4715.0/64.0 ,
468  bi[8][3] = 2501.0/64.0 ,
469  bi[8][2] = 149.0/16.0 ,
470  bi[8][1] = 1.0 ,
471  //
472  bi[9][6] = 55.0/3.0 ,
473  bi[9][5] = 71.0 ,
474  bi[9][4] = 103.0 ,
475  bi[9][3] = 199.0/3.0 ,
476  bi[9][2] = 16.0 ,
477  //
478  bi[10][6] = -1774004627.0/75810735.0 ,
479  bi[10][5] = -1774004627.0/25270245.0 ,
480  bi[10][4] = -26477681.0/359975.0 ,
481  bi[10][3] = -11411880511.0/379053675.0 ,
482  bi[10][2] = -423642896.0/126351225.0 ,
483  //
484  bi[11][6] = 35.0 ,
485  bi[11][5] = 105.0 ,
486  bi[11][4] = 117.0 ,
487  bi[11][3] = 59.0 ,
488  bi[11][2] = 12.0 ;
489  //
490 
491 
492 // for (int i = 1; i <= 11; i++) {
493 // bi[i][1] = 0.e0;
494 // }
495 // for (int i = 1; i <= 6; i++) {
496 // bi[2][i] = 0.e0;
497 // }
498 // bi[1][6] = -12134338393.e0/1050809760.e0;
499 // bi[1][5] = -1620741229.e0/50038560.e0;
500 // bi[1][4] = -2048058893.e0/59875200.e0;
501 // bi[1][3] = -87098480009.e0/5254048800.e0;
502 // bi[1][2] = -11513270273.e0/3502699200.e0;
503 // //C
504 // bi[3][6] = -33197340367.e0/1218433216.e0;
505 // bi[3][5] = -539868024987.e0/6092166080.e0;
506 // bi[3][4] = -39991188681.e0/374902528.e0;
507 // bi[3][3] = -69509738227.e0/1218433216.e0;
508 // bi[3][2] = -29327744613.e0/2436866432.e0;
509 // //C
510 // bi[4][6] = -284800997201.e0/19905339168.e0;
511 // bi[4][5] = -7896875450471.e0/165877826400.e0;
512 // bi[4][4] = -333945812879.e0/5671036800.e0;
513 // bi[4][3] = -16209923456237.e0/497633479200.e0;
514 // bi[4][2] = -2382590741699.e0/331755652800.e0;
515 // //C
516 // bi[5][6] = -540919.e0/741312.e0;
517 // bi[5][5] = -103626067.e0/43243200.e0;
518 // bi[5][4] = -633779.e0/211200.e0;
519 // bi[5][3] = -32406787.e0/18532800.e0;
520 // bi[5][2] = -36591193.e0/86486400.e0;
521 // //C
522 // bi[6][6] = 7157998304.e0/374350977.e0;
523 // bi[6][5] = 30405842464.e0/623918295.e0;
524 // bi[6][4] = 183022264.e0/5332635.e0;
525 // bi[6][3] = -3357024032.e0/1871754885.e0;
526 // bi[6][2] = -611586736.e0/89131185.e0;
527 // //C
528 // bi[7][6] = -138073.e0/9408.e0;
529 // bi[7][5] = -719433.e0/15680.e0;
530 // bi[7][4] = -1620541.e0/31360.e0;
531 // bi[7][3] = -385151.e0/15680.e0;
532 // bi[7][2] = -65403.e0/15680.e0;
533 // //C
534 // bi[8][6] = 1245.e0/64.e0;
535 // bi[8][5] = 3991.e0/64.e0;
536 // bi[8][4] = 4715.e0/64.e0;
537 // bi[8][3] = 2501.e0/64.e0;
538 // bi[8][2] = 149.e0/16.e0;
539 // bi[8][1] = 1.e0;
540 // //C
541 // bi[9][6] = 55.e0/3.e0;
542 // bi[9][5] = 71.e0;
543 // bi[9][4] = 103.e0;
544 // bi[9][3] = 199.e0/3.e0;
545 // bi[9][2] = 16.0e0;
546 // //C
547 // bi[10][6] = -1774004627.e0/75810735.e0;
548 // bi[10][5] = -1774004627.e0/25270245.e0;
549 // bi[10][4] = -26477681.e0/359975.e0;
550 // bi[10][3] = -11411880511.e0/379053675.e0;
551 // bi[10][2] = -423642896.e0/126351225.e0;
552 // //C
553 // bi[11][6] = 35.e0;
554 // bi[11][5] = 105.e0;
555 // bi[11][4] = 117.e0;
556 // bi[11][3] = 59.e0;
557 // bi[11][2] = 12.e0;
558 
559 
560 
561  const G4int numberOfVariables= this->GetNumberOfVariables();
562 
563  // Saving yInput because yInput and yOut can be aliases for same array
564  for(int i=0;i<numberOfVariables;i++)
565  {
566  yIn[i]=yInput[i];
567  }
568 
569  // The number of variables to be integrated over
570  yOut[7] = yTemp[7] = yIn[7];
571 
572 
573 
574  // calculating extra stages
575 
576  for(int i=0; i<numberOfVariables; i++){
577  yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
578  a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
579  a97*ak7[i] + a98*ak8[i] );
580  }
581 
582  RightHandSide(yTemp, ak9);
583 
584  for(int i=0; i<numberOfVariables; i++){
585  yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
586  a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
587  a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
588  }
589 
590  RightHandSide(yTemp, ak10);
591 
592  for(int i=0; i<numberOfVariables; i++){
593  yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
594  a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
595  a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
596  a1110*ak10[i] );
597  }
598 
599  RightHandSide(yTemp, ak11);
600 
601 
602  G4double tau0 = tau;
603  // Calculating the polynomials :
604  for(int i=1; i<=11; i++){ //Here i is NOT the coordinate no. , it's stage no.
605  b[i] = 0.0;
606  tau = tau0;
607  for(int j=1; j<=6; j++){
608  b[i] += bi[i][j]*tau;
609  tau*=tau0;
610  }
611  }
612 
613  for(int i=0; i<numberOfVariables; i++){
614  yOut[i] = yIn[i] + Step*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
615  b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
616  b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
617  b[10]*ak10[i] + b[11]*ak11[i] );
618  }
619 
620 }
621 
622 
623 
624 
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition: G4LineSection.hh:84
G4VFSALIntegrationStepper::GetNumberOfStateVariables
G4int GetNumberOfStateVariables() const
G4LineSection.hh
G4FSALBogackiShampine45::G4FSALBogackiShampine45
G4FSALBogackiShampine45(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition: G4FSALBogackiShampine45.cc:68
G4int
int G4int
Definition: G4Types.hh:78
G4VFSALIntegrationStepper::RightHandSide
void RightHandSide(const double y[], double dydx[])
Definition: G4VFSALIntegrationStepper.cc:60
G4FSALBogackiShampine45::interpolate
void interpolate(const G4double yInput[], const G4double dydx[], G4double yOut[], G4double Step, G4double tau)
Definition: G4FSALBogackiShampine45.cc:345
test.b
tuple b
Definition: test.py:12
G4FSALBogackiShampine45::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[], G4double nextDydx[])
Definition: G4FSALBogackiShampine45.cc:153
G4bool
bool G4bool
Definition: G4Types.hh:79
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
G4VFSALIntegrationStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const
G4FSALBogackiShampine45.hh
Step
Definition: Step.hh:41
G4double
double G4double
Definition: G4Types.hh:76