Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4BogackiShampine45.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 // Bogacki-Shampine's RK 5(4) non-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: 25 May 2015
32 //
33 // History
34 // -----------------------------
35 // Created by Somnath Banerjee on May-August 2015
36 // Improvements by John Apostolakis, May 2016
38 //
39 // This is the source file of G4BogackiShampine45 class containing the
40 // definition of the stepper() method that evaluates one step in
41 // field propagation.
42 //
43 // The Butcher table of the Bogacki-Shampine-8-4-5 method is:
44 //
45 // 0 |
46 // 1/6 | 1/6
47 // 2/9 | 2/27 4/27
48 // 3/7 | 183/1372 -162/343 1053/1372
49 // 2/3 | 68/297 -4/11 42/143 1960/3861
50 // 3/4 | 597/22528 81/352 63099/585728 58653/366080 4617/20480
51 // 1 | 174197/959244 -30942/79937 8152137/19744439 666106/1039181 -29421/29068 482048/414219
52 // 1 | 587/8064 0 4440339/15491840 24353/124800 387/44800 2152/5985 7267/94080
53 //-------------------------------------------------------------------------------------------------------------------
54 // 587/8064 0 4440339/15491840 24353/124800 387/44800 2152/5985 7267/94080 0
55 // 2479/34992 0 123/416 612941/3411720 43/1440 2272/6561 79937/1113912 3293/556956
56 //
57 // Do NOT re-indent the lines above - their meaning becomes lost
58 // ********************************
59 // Coefficients have been obtained from rksuite.f : http://www.netlib.org/ode/rksuite/
60 
61 // Note on meaning of label "non-FSAL version":
62 // This method calculates the deriviative dy/dx at the endpoint of the integration interval at each step.
63 // as part of its evaluation of the endpoint and its error.
64 // So this value is available to be returned, for re-use in case of a successful step.
65 // ( This is done in a 'later' version using a refined interface. )
66 
67 #include <cassert>
68 
69 #include "G4BogackiShampine45.hh"
70 #include "G4LineSection.hh"
71 
74 
75 //Constructor
77  G4int noIntegrationVariables,
78  G4bool primary)
79  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
80  fLastStepLength(-1.0),
81  fAuxStepper(nullptr),
82  fPreparedInterpolation(false)
83 {
84 
85  const G4int numberOfVariables = noIntegrationVariables;
86 
87  //New Chunk of memory being created for use by the stepper
88 
89  //aki - for storing intermediate RHS
90  ak2 = new G4double[numberOfVariables];
91  ak3 = new G4double[numberOfVariables];
92  ak4 = new G4double[numberOfVariables];
93  ak5 = new G4double[numberOfVariables];
94  ak6 = new G4double[numberOfVariables];
95  ak7 = new G4double[numberOfVariables];
96  ak8 = new G4double[numberOfVariables];
97  ak9 = new G4double[numberOfVariables];
98  ak10 = new G4double[numberOfVariables];
99  ak11 = new G4double[numberOfVariables];
100 
101  for (int i = 0; i < 6; i++) {
102  p[i]= new G4double[numberOfVariables];
103  }
104 
105  assert ( GetNumberOfStateVariables() >= 8 );
106  const G4int numStateVars = std::max(noIntegrationVariables,
108 
109  // Must ensure space extra 'state' variables exists - i.e. yIn[7]
110  yTemp = new G4double[numStateVars];
111  yIn = new G4double[numStateVars] ;
112 
113  fLastInitialVector = new G4double[numStateVars] ;
114  fLastFinalVector = new G4double[numStateVars] ;
115  fLastDyDx = new G4double[numberOfVariables]; // Only derivatives
116 
117  fMidVector = new G4double[numberOfVariables]; // new G4double[numStateVars];
118  fMidError = new G4double[numberOfVariables]; // new G4double[numStateVars];
119 
120  if( ! fPreparedConstants )
122 
123  if( primary )
124  {
125  fAuxStepper = new G4BogackiShampine45(EqRhs, numberOfVariables, false);
126  }
127 }
128 
129 
130 //Destructor
132  //clear all previously allocated memory for stepper and DistChord
133  delete[] ak2;
134  delete[] ak3;
135  delete[] ak4;
136  delete[] ak5;
137  delete[] ak6;
138  delete[] ak7;
139  delete[] ak8;
140  delete[] ak9;
141  delete[] ak10;
142  delete[] ak11;
143 
144  for (int i = 0; i < 6; i++) {
145  delete[] p[i];
146  }
147 
148  delete[] yTemp;
149  delete[] yIn;
150 
151  delete[] fLastInitialVector;
152  delete[] fLastFinalVector;
153  delete[] fLastDyDx;
154  delete[] fMidVector;
155  delete[] fMidError;
156 
157  delete fAuxStepper;
158 }
159 
160 // G4double* G4BogackiShampine45::getLastDydx(){
161 // return ak8;
162 // }
163 
164 void
166 {
167  const G4int numberOfVariables= this->GetNumberOfVariables();
168 
169  for(G4int i=0; i < numberOfVariables; i++ ){
170  dyDxLast[i] = ak9[i];
171  }
172 }
173 
174 //Stepper :
175 
176 // Passing in the value of yInput[],the first time dydx[] and Step length
177 // Giving back yOut and yErr arrays for output and error respectively
178 
180  const G4double DyDx[],
181  G4double Step,
182  G4double yOut[],
183  G4double yErr[] )
184 {
185  G4int i;
186 
187  // Constants from the Butcher tableu
188  const G4double
189  b21 = 1.0/6.0 ,
190  b31 = 2.0/27.0 , b32 = 4.0/27.0,
191 
192  b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
193 
194  b51 = 68.0/297.0, b52 = -4.0/11.0,
195  b53 = 42.0/143.0, b54 = 1960.0/3861.0,
196 
197  b61 = 597.0/22528.0, b62 = 81.0/352.0,
198  b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
199  b65 = 4617.0/20480.0,
200 
201  b71 = 174197.0/959244.0, b72 = -30942.0/79937.0,
202  b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
203  b75 = -29421.0/29068.0, b76 = 482048.0/414219.0,
204 
205  b81 = 587.0/8064.0, b82 = 0.0,
206  b83 = 4440339.0/15491840.0, b84 = 24353.0/124800.0,
207  b85 = 387.0/44800.0, b86 = 2152.0/5985.0,
208  b87 = 7267.0/94080.0;
209 
210 // c1 = 2479.0/34992.0,
211 // c2 = 0.0,
212 // c3 = 123.0/416.0,
213 // c4 = 612941.0/3411720.0,
214 // c5 = 43.0/1440.0,
215 // c6 = 2272.0/6561.0,
216 // c7 = 79937.0/1113912.0,
217 // c8 = 3293.0/556956.0,
218 
219  // For the embedded higher order method only the difference of values
220  // taken and is used directly later (instead of defining the last row
221  // of Butcher table in separate constants and taking the
222  // difference)
223  const G4double
224  dc1 = b81 - 2479.0 / 34992.0 ,
225  dc2 = 0.0,
226  dc3 = b83 - 123.0 / 416.0 ,
227  dc4 = b84 - 612941.0 / 3411720.0,
228  dc5 = b85 - 43.0 / 1440.0,
229  dc6 = b86 - 2272.0 / 6561.0,
230  dc7 = b87 - 79937.0 / 1113912.0,
231  dc8 = - 3293.0 / 556956.0;
232 
233  const G4int numberOfVariables= this->GetNumberOfVariables();
234 
235  // The number of variables to be integrated over
236  yOut[7] = yTemp[7] = yIn[7];
237  // Saving yInput because yInput and yOut can be aliases for same array
238 
239  for(i=0;i<numberOfVariables;i++)
240  {
241  yIn[i]=yInput[i];
242  }
243 
244  // RightHandSide(yIn, dydx) ;
245  // 1st Step - Not doing, getting passed
246 
247  for(i=0;i<numberOfVariables;i++)
248  {
249  yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
250  }
251  RightHandSide(yTemp, ak2) ; // 2nd Step
252 
253  for(i=0;i<numberOfVariables;i++)
254  {
255  yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
256  }
257  RightHandSide(yTemp, ak3) ; // 3rd Step
258 
259  for(i=0;i<numberOfVariables;i++)
260  {
261  yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
262  }
263  RightHandSide(yTemp, ak4) ; // 4th Step
264 
265  for(i=0;i<numberOfVariables;i++)
266  {
267  yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
268  b54*ak4[i]) ;
269  }
270  RightHandSide(yTemp, ak5) ; // 5th Step
271 
272  for(i=0;i<numberOfVariables;i++)
273  {
274  yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
275  b64*ak4[i] + b65*ak5[i]) ;
276  }
277  RightHandSide(yTemp, ak6) ; // 6th Step
278 
279  for(i=0;i<numberOfVariables;i++)
280  {
281  yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
282  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
283  }
284  RightHandSide(yTemp, ak7); //7th Step
285 
286  for(i=0;i<numberOfVariables;i++)
287  {
288  yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] +
289  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
290  b87*ak7[i]);
291  }
292  RightHandSide(yOut, ak8); //8th Step - Final one Using FSAL
293 
294  for(i=0;i<numberOfVariables;i++)
295  {
296  yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
297  dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
298 
299  // Store Input and Final values, for possible use in calculating chord
300  fLastInitialVector[i] = yIn[i] ;
301  fLastFinalVector[i] = yOut[i];
302  fLastDyDx[i] = DyDx[i];
303  }
304 
305  fLastStepLength = Step;
306  fPreparedInterpolation= false;
307 
308  return ;
309 }
310 
311 
312 //The following has not been tested
313 
314 //The DistChord() function fot the class - must define it here.
316 {
317  G4double distLine, distChord;
318  G4ThreeVector initialPoint, finalPoint, midPoint;
319 
320  // Store last initial and final points (they will be overwritten in self-Stepper call!)
321  initialPoint = G4ThreeVector( fLastInitialVector[0],
323  finalPoint = G4ThreeVector( fLastFinalVector[0],
325 
326 #if 1
327  // Old method -- Do half a step using StepNoErr
330 #else
331  // New method -- Using interpolation, requires only 3 extra stages (ie 3 extra field evaluations )
332 
333  // Use Interpolation, instead of auxiliary stepper to evaluate midpoint
334  if( ! fPreparedInterpolation ) {
335  G4BogackiShampine45 *cThis= const_cast<G4BogackiShampine45 *>(this);
336  cThis-> SetupInterpolationHigh(); // ( fLastInitialVector, fLastDyDx, fLastStepLength );
337  }
338  //For calculating the output at the tau fraction of Step
339  G4double tau = 0.5;
340  // cThis->InterpolateHigh( /* fLastInitialVector, fLastDyDx, fLastStepLength, */, fMidVector, tau);
341  // Old arguments: ( /*yInput, dydx, step,*/ yOut, tau );
342  this->InterpolateHigh( tau, fMidVector );
343 #endif
344 
345  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
346 
347  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
348  // distance of Chord
349 
350  if (initialPoint != finalPoint)
351  {
352  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
353  distChord = distLine;
354  }
355  else
356  {
357  distChord = (midPoint-initialPoint).mag();
358  }
359  return distChord;
360 }
361 
363  // ( const G4double *yInput, const G4double *dydx, const G4double Step)
364 {
365  //Coefficients for the additional stages :
366  const G4double
367  a91 = 455.0/6144.0 ,
368  a92 = 0.0 ,
369  a93 = 10256301.0/35409920.0 ,
370  a94 = 2307361.0/17971200.0 ,
371  a95 = -387.0/102400.0 ,
372  a96 = 73.0/5130.0 ,
373  a97 = -7267.0/215040.0 ,
374  a98 = 1.0/32.0 ,
375 
376  a101 = -837888343715.0/13176988637184.0 ,
377  a102 = 30409415.0/52955362.0 ,
378  a103 = -48321525963.0/759168069632.0 ,
379  a104 = 8530738453321.0/197654829557760.0 ,
380  a105 = 1361640523001.0/1626788720640.0 ,
381  a106 = -13143060689.0/38604458898.0 ,
382  a107 = 18700221969.0/379584034816.0 ,
383  a108 = -5831595.0/847285792.0 ,
384  a109 = -5183640.0/26477681.0 ,
385 
386  a111 = 98719073263.0/1551965184000.0 ,
387  a112 = 1307.0/123552.0 ,
388  a113 = 4632066559387.0/70181753241600.0 ,
389  a114 = 7828594302389.0/382182512025600.0 ,
390  a115 = 40763687.0/11070259200.0 ,
391  a116 = 34872732407.0/224610586200.0 ,
392  a117 = -2561897.0/30105600.0 ,
393  a118 = 1.0/10.0 ,
394  a119 = -1.0/10.0 ,
395  a1110 = -1403317093.0/11371610250.0 ;
396 
397  const G4int numberOfVariables= this->GetNumberOfVariables();
398 
399  // const G4double *yIn= fLastInitialVector;
400  const G4double *dydx= fLastDyDx;
401  const G4double Step = fLastStepLength;
402 
403  yTemp[7] = yIn[7];
404 
405  //Evaluate the extra stages :
406  for(int i=0; i<numberOfVariables; i++){
407  yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
408  a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
409  a97*ak7[i] + a98*ak8[i] );
410  }
411 
412  RightHandSide(yTemp, ak9); //9th stage
413 
414  for(int i=0; i<numberOfVariables; i++){
415  yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
416  a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
417  a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
418  }
419 
420  RightHandSide(yTemp, ak10); //10th stage
421 
422  for(int i=0; i<numberOfVariables; i++){
423  yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
424  a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
425  a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
426  a1110*ak10[i] );
427  }
428  RightHandSide(yTemp, ak11); //11th stage
429 
430  // In future we can restrict the number of variables interpolated
431  int nwant = numberOfVariables;
432 
433 // Form the coefficients of the interpolating polynomial in its shifted
434 // and scaled form. The terms are grouped to minimize the errors
435 // of the transformation, to cope with ill-conditioning. ( From RKSUITE )
436 //
437  for (int l = 0; l < nwant; l++) {
438  // Coefficient of tau^6
439  p[5][l] = bi[5][6]*ak5[l] +
440  ((bi[10][6]*ak10[l] + bi[8][6]*ak8[l]) +
441  (bi[7][6]*ak7[l] + bi[6][6]*ak6[l])) +
442  ((bi[4][6]*ak4[l] + bi[9][6]*ak9[l]) +
443  (bi[3][6]*ak3[l] + bi[11][6]*ak11[l]) +
444  bi[1][6]*dydx[l]);
445  // Coefficient of tau^5
446  p[4][l] = (bi[10][5]*ak10[l] + bi[9][5]*ak9[l]) +
447  ((bi[7][5]*ak7[l] + bi[6][5]*ak6[l]) +
448  bi[5][5]*ak5[l]) + ((bi[4][5]*ak4[l] +
449  bi[8][5]*ak8[l]) + (bi[3][5]*ak3[l] +
450  bi[11][5]*ak11[l]) + bi[1][5]*dydx[l]);
451  // Coefficient of tau^4
452  p[3][l] = ((bi[4][4]*ak4[l] + bi[8][4]*ak8[l]) +
453  (bi[7][4]*ak7[l] + bi[6][4]*ak6[l]) +
454  bi[5][4]*ak5[l]) + ((bi[10][4]*ak10[l] +
455  bi[9][4]*ak9[l]) + (bi[3][4]*ak3[l] +
456  bi[11][4]*ak11[l]) + bi[1][4]*dydx[l]);
457  // Coefficient of tau^3
458  p[2][l] = bi[5][3]*ak5[l] + bi[6][3]*ak6[l] +
459  ((bi[3][3]*ak3[l] + bi[9][3]*ak9[l]) +
460  (bi[10][3]*ak10[l]+ bi[8][3]*ak8[l]) + bi[1][3]*dydx[l]) +
461  ((bi[4][3]*ak4[l] + bi[11][3]*ak11[l]) + bi[7][3]*ak7[l]);
462  // Coefficient of tau^2
463  p[1][l] = bi[5][2]*ak5[l] + ((bi[6][2]*ak6[l] +
464  bi[8][2]*ak8[l]) + bi[1][2]*dydx[l]) +
465  ((bi[3][2]*ak3[l] + bi[9][2]*ak9[l]) +
466  bi[10][2]*ak10[l])+ ((bi[4][2]*ak4[l] +
467  bi[11][2]*ak2[l]) + bi[7][2]*ak7[l]);
468  }
469 //
470 // Scale all the coefficients by the step size.
471 //
472  for (int i = 0; i < 6; i++) {
473  for (int l = 0; l < nwant; l++) {
474  p[i][l] *= Step;
475  }
476  }
477 
479 }
480 
482 {
483  for(int i=1; i<= 11; i++)
484  bi[i][1] = 0.0 ;
485 
486  for(int i=1; i<=6; i++)
487  bi[2][i] = 0.0 ;
488 
489  bi[1][6] = -12134338393.0 / 1050809760.0 ,
490  bi[1][5] = -1620741229.0 / 50038560.0 ,
491  bi[1][4] = -2048058893.0 / 59875200.0 ,
492  bi[1][3] = -87098480009.0 / 5254048800.0 ,
493  bi[1][2] = -11513270273.0 / 3502699200.0 ,
494  //
495  bi[3][6] = -33197340367.0 / 1218433216.0 ,
496  bi[3][5] = -539868024987.0 / 6092166080.0 ,
497  bi[3][4] = -39991188681.0 / 374902528.0 ,
498  bi[3][3] = -69509738227.0 / 1218433216.0 ,
499  bi[3][2] = -29327744613.0 / 2436866432.0 ,
500  //
501  bi[4][6] = -284800997201.0 / 19905339168.0 ,
502  bi[4][5] = -7896875450471.0 / 165877826400.0 ,
503  bi[4][4] = -333945812879.0 / 5671036800.0 ,
504  bi[4][3] = -16209923456237.0 / 497633479200.0 ,
505  bi[4][2] = -2382590741699.0 / 331755652800.0 ,
506  //
507  bi[5][6] = -540919.0 / 741312.0 ,
508  bi[5][5] = -103626067.0 / 43243200.0 ,
509  bi[5][4] = -633779.0 / 211200.0 ,
510  bi[5][3] = -32406787.0 / 18532800.0 ,
511  bi[5][2] = -36591193.0 / 86486400.0 ,
512  //
513  bi[6][6] = 7157998304.0 / 374350977.0 ,
514  bi[6][5] = 30405842464.0 / 623918295.0 ,
515  bi[6][4] = 183022264.0 / 5332635.0 ,
516  bi[6][3] = -3357024032.0 / 1871754885.0 ,
517  bi[6][2] = -611586736.0 / 89131185.0 ,
518  //
519  bi[7][6] = -138073.0 / 9408.0 ,
520  bi[7][5] = -719433.0 / 15680.0 ,
521  bi[7][4] = -1620541.0 / 31360.0 ,
522  bi[7][3] = -385151.0 / 15680.0 ,
523  bi[7][2] = -65403.0 / 15680.0 ,
524  //
525  bi[8][6] = 1245.0 / 64.0 ,
526  bi[8][5] = 3991.0 / 64.0 ,
527  bi[8][4] = 4715.0 / 64.0 ,
528  bi[8][3] = 2501.0 / 64.0 ,
529  bi[8][2] = 149.0 / 16.0 ,
530  bi[8][1] = 1.0 ,
531  //
532  bi[9][6] = 55.0 / 3.0 ,
533  bi[9][5] = 71.0 ,
534  bi[9][4] = 103.0 ,
535  bi[9][3] = 199.0 / 3.0 ,
536  bi[9][2] = 16.0 ,
537  //
538  bi[10][6] = -1774004627.0 / 75810735.0 ,
539  bi[10][5] = -1774004627.0 / 25270245.0 ,
540  bi[10][4] = -26477681.0 / 359975.0 ,
541  bi[10][3] = -11411880511.0 / 379053675.0 ,
542  bi[10][2] = -423642896.0 / 126351225.0 ,
543  //
544  bi[11][6] = 35.0 ,
545  bi[11][5] = 105.0 ,
546  bi[11][4] = 117.0 ,
547  bi[11][3] = 59.0 ,
548  bi[11][2] = 12.0 ;
549  fPreparedConstants= true;
550 }
551 
553 // ( const G4double *yInput, const G4double *dydx, const G4double Step, G4double *yOut, G4double tau)
554 {
555  G4int numberOfVariables = this->GetNumberOfVariables();
556  assert( fPreparedConstants);
557 
558  G4Exception("G4BogackiShampine45::InterpolateHigh()", "GeomField0001",
559  FatalException, "Method is not yet validated.");
560 
561  // const G4double *yIn= fLastInitialVector;
562  // const G4double *dydx= fLastDyDx;
563  const G4double Step = fLastStepLength;
564 
565  // for(G4int i = 0; i< numberOfVariables; i++) yIn[i] = yInput[i];
566 #if 1
567  G4int nwant = numberOfVariables;
568  const G4int norder= 6;
569  G4int l, k;
570 
571  for (l = 0; l < nwant; l++) {
572  yOut[l] = p[norder-1][l] * tau;
573  }
574  for (k = norder - 2; k >= 1; k--) {
575  for (l = 0; l < nwant; l++) {
576  yOut[l] = ( yOut[l] + p[k][l] ) * tau;
577  }
578  }
579  for (l = 0; l < nwant; l++) {
580  yOut[l] = ( yOut[l] + Step * ak8[l] ) * tau + yIn[l];
581  }
582  // The derivative at the end-point is nextDydx[i] = ak8[i];
583 #else
584  // The scheme tries to do the same as the DormandPrince745 routine, but fails
585  G4double b[12];
586  const G4double *dydx= fLastDyDx;
587 
588  G4double tau0 = tau;
589 
590  for(int iStage=1; iStage<=11; iStage++){ // iStage = stage number
591  b[iStage] = 0.0;
592  tau = tau0;
593  for(int j=6; j>=1; j--){ // j reversed
594  b[iStage] += bi[iStage][j] * tau;
595  tau *= tau0;
596  }
597  }
598 
599  for(int i=0; i<numberOfVariables; i++){
600  yOut[i] = yIn[i] + Step*(b[1] * dydx[i] + b[2] * ak2[i] + b[3] * ak3[i] +
601  b[4] * ak4[i] + b[5] * ak5[i] + b[6] * ak6[i] +
602  b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i] +
603  b[10] * ak10[i] + b[11] * ak11[i] );
604  }
605 #endif
606 }
T max(const T t1, const T t2)
brief Return the largest of the two arguments
CLHEP::Hep3Vector G4ThreeVector
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
void InterpolateHigh(G4double tau, G4double yOut[]) const
const char * p
Definition: xmltok.h:285
G4double DistChord() const
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
G4BogackiShampine45 * fAuxStepper
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
void GetLastDydx(G4double dyDxLast[])
G4int GetNumberOfStateVariables() const
void RightHandSide(const double y[], double dydx[]) const
static G4double bi[12][7]
G4BogackiShampine45(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
int G4int
Definition: G4Types.hh:78
static const G4double ak2
G4int GetNumberOfVariables() const
Definition: Step.hh:41