Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4DormandPrince745.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 // $Id: G4DormandPrince745.cc 107470 2017-11-15 07:14:28Z gcosmo $
27 //
28 // Class description:
29 //
30 // DormandPrince7 - 5(4) non-FSAL
31 //
32 // This is the source file of G4DormandPrince745 class containing the
33 // definition of the stepper() method that evaluates one step in
34 // field propagation.
35 // The coefficients and the algorithm have been adapted from
36 //
37 // Table 2 : Coefficients of RK5(4)7M
38 // ---Ref---
39 // J. R. Dormand and P. J. Prince, “A family of embedded Runge-Kutta formulae,”
40 // Journal of computational and applied …, vol. 6, no. 1, pp. 19–26, 1980.
41 // ------------------
42 //
43 // The Butcher table of the Dormand-Prince-7-4-5 method is as follows :
44 //
45 // 0 |
46 // 1/5 | 1/5
47 // 3/10| 3/40 9/40
48 // 4/5 | 44/45 −56/15 32/9
49 // 8/9 | 19372/6561 −25360/2187 64448/6561 −212/729
50 // 1 | 9017/3168 −355/33 46732/5247 49/176 −5103/18656
51 // 1 | 35/384 0 500/1113 125/192 −2187/6784 11/84
52 // ------------------------------------------------------------------------
53 // 35/384 0 500/1113 125/192 −2187/6784 11/84 0
54 // 5179/57600 0 7571/16695 393/640 −92097/339200 187/2100 1/40
55 //
56 //
57 // Implementation by Somnath Banerjee - GSoC 2015
58 // Work supported by Google as part of Google Summer of Code 2015.
59 // Supervision / code review: John Apostolakis
60 //
61 // First version: 25 May 2015 - Somnath Banerjee
62 //
63 // Note: Current version includes 3 versions of 'DistChord' method.
64 // Default is hard-coded interpolation.
65 //
66 #include "G4DormandPrince745.hh"
67 #include "G4LineSection.hh"
68 #include <cmath>
69 
70 //Constructor
72  G4int noIntegrationVariables,
73  G4bool primary)
74 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
75  fAuxStepper(0)
76 {
77  const G4int numberOfVariables = // noIntegrationVariables;
78  std::max( noIntegrationVariables,
79  ( ( (noIntegrationVariables-1)/4 + 1 ) * 4 ) );
80  // For better alignment with cache-line
81 
82  //New Chunk of memory being created for use by the stepper
83 
84  //ak_i - for storing intermediate RHS
85  ak2 = new G4double[numberOfVariables];
86  ak3 = new G4double[numberOfVariables];
87  ak4 = new G4double[numberOfVariables];
88  ak5 = new G4double[numberOfVariables];
89  ak6 = new G4double[numberOfVariables];
90  ak7 = new G4double[numberOfVariables];
91  // Also always allocate arrays for interpolation stages
92  ak8 = new G4double[numberOfVariables];
93  ak9 = new G4double[numberOfVariables];
94 
95  // Must ensure space for extra 'state' variables exists - i.e. yIn[7]
96  const G4int numStateVars =
97  std::max(noIntegrationVariables,
99  );
100  yTemp = new G4double[numStateVars];
101  yIn = new G4double[numStateVars];
102 
103  fLastInitialVector = new G4double[numStateVars] ;
104  fLastFinalVector = new G4double[numStateVars] ;
105 
106  // fLastDyDx = new G4double[numberOfVariables];
107 
108  fMidVector = new G4double[numStateVars];
109  fMidError = new G4double[numStateVars];
110 
111  yTemp = new G4double[numberOfVariables] ;
112  yIn = new G4double[numberOfVariables] ;
113 
114  fLastInitialVector = new G4double[numberOfVariables] ;
115  fLastFinalVector = new G4double[numberOfVariables] ;
116  fInitialDyDx = new G4double[numberOfVariables];
117 
118  fMidVector = new G4double[numberOfVariables];
119  fMidError = new G4double[numberOfVariables];
120  fAuxStepper = nullptr;
121  if( primary )
122  {
123  fAuxStepper = new G4DormandPrince745(EqRhs, numberOfVariables,
124  !primary);
125  }
126  fLastStepLength = -1.0;
127 }
128 
129 //Destructor
131 {
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  // Used only for interpolation
140  delete[] ak8;
141  delete[] ak9;
142 
143  delete[] yTemp;
144  delete[] yIn;
145 
146  delete[] fLastInitialVector;
147  delete[] fLastFinalVector;
148  delete[] fInitialDyDx;
149  delete[] fMidVector;
150  delete[] fMidError;
151 
152  delete fAuxStepper;
153 }
154 
155 
156 // The coefficients and the algorithm have been adapted from
157 // Table 2 : Coefficients of RK5(4)7M
158 // ---Ref---
159 // J. R. Dormand and P. J. Prince, “A family of embedded Runge-Kutta formulae,”
160 // Journal of computational and applied …, vol. 6, no. 1, pp. 19–26, 1980.
161 // ------------------
162 
163 // The Butcher table of the Dormand-Prince-7-4-5 method is as follows :
164 //
165 // 0 |
166 // 1/5 | 1/5
167 // 3/10| 3/40 9/40
168 // 4/5 | 44/45 −56/15 32/9
169 // 8/9 | 19372/6561 −25360/2187 64448/6561 −212/729
170 // 1 | 9017/3168 −355/33 46732/5247 49/176 −5103/18656
171 // 1 | 35/384 0 500/1113 125/192 −2187/6784 11/84
172 // ------------------------------------------------------------------------
173 // 35/384 0 500/1113 125/192 −2187/6784 11/84 0
174 // 5179/57600 0 7571/16695 393/640 −92097/339200 187/2100 1/40
175 
176 
177 //Stepper :
178 
179 // Passing in the value of yInput[],the first time dydx[] and Step length
180 // Giving back yOut and yErr arrays for output and error respectively
181 
183  const G4double DyDx[],
184  G4double Step,
185  G4double yOut[],
186  G4double yErr[] )
187 {
188  G4int i;
189 
190  //The various constants defined on the basis of butcher tableu
191  const G4double //G4double - only once
192  b21 = 0.2 ,
193 
194  b31 = 3.0/40.0, b32 = 9.0/40.0 ,
195 
196  b41 = 44.0/45.0, b42 = -56.0/15.0, b43 = 32.0/9.0,
197 
198  b51 = 19372.0/6561.0, b52 = -25360.0/2187.0, b53 = 64448.0/6561.0,
199  b54 = -212.0/729.0 ,
200 
201  b61 = 9017.0/3168.0 , b62 = -355.0/33.0,
202  b63 = 46732.0/5247.0 , b64 = 49.0/176.0 ,
203  b65 = -5103.0/18656.0 ,
204 
205  b71 = 35.0/384.0, b72 = 0.,
206  b73 = 500.0/1113.0, b74 = 125.0/192.0,
207  b75 = -2187.0/6784.0, b76 = 11.0/84.0,
208 
209  //Sum of columns, sum(bij) = ei
210 // e1 = 0. ,
211 // e2 = 1.0/5.0 ,
212 // e3 = 3.0/10.0 ,
213 // e4 = 4.0/5.0 ,
214 // e5 = 8.0/9.0 ,
215 // e6 = 1.0 ,
216 // e7 = 1.0 ,
217 
218 // Difference between the higher and the lower order method coeff. :
219  // b7j are the coefficients of higher order
220 
221  dc1 = -( b71 - 5179.0/57600.0),
222  dc2 = -( b72 - .0),
223  dc3 = -( b73 - 7571.0/16695.0),
224  dc4 = -( b74 - 393.0/640.0),
225  dc5 = -( b75 + 92097.0/339200.0),
226  dc6 = -( b76 - 187.0/2100.0),
227  dc7 = -( - 1.0/40.0 ); //end of declaration
228 
229  const G4int numberOfVariables= this->GetNumberOfVariables();
230 
231  // The number of variables to be integrated over
232  yOut[7] = yTemp[7] = yInput[7];
233  // Saving yInput because yInput and yOut can be aliases for same array
234 
235  for(i=0;i<numberOfVariables;i++)
236  {
237  yIn[i]=yInput[i];
238  }
239 
240  // RightHandSide(yIn, DyDx) ;
241  // 1st Step - Not doing, getting passed
242 
243  for(i=0;i<numberOfVariables;i++)
244  {
245  yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
246  }
247  RightHandSide(yTemp, ak2) ; // 2nd stage
248 
249  for(i=0;i<numberOfVariables;i++)
250  {
251  yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
252  }
253  RightHandSide(yTemp, ak3) ; // 3rd stage
254 
255  for(i=0;i<numberOfVariables;i++)
256  {
257  yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
258  }
259  RightHandSide(yTemp, ak4) ; // 4th stage
260 
261  for(i=0;i<numberOfVariables;i++)
262  {
263  yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
264  b54*ak4[i]) ;
265  }
266  RightHandSide(yTemp, ak5) ; // 5th stage
267 
268  for(i=0;i<numberOfVariables;i++)
269  {
270  yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
271  b64*ak4[i] + b65*ak5[i]) ;
272  }
273  RightHandSide(yTemp, ak6) ; // 6th stage
274 
275  for(i=0;i<numberOfVariables;i++)
276  {
277  yOut[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
278  b74*ak4[i] + b75*ak5[i] + b76*ak6[i] );
279  }
280  RightHandSide(yOut, ak7); //7th and Final stage
281 
282  for(i=0;i<numberOfVariables;i++)
283  {
284  yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
285  dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) + 1.5e-18 ;
286 
287  // Store Input and Final values, for possible use in calculating chord
288  fLastInitialVector[i] = yIn[i] ;
289  fLastFinalVector[i] = yOut[i];
290  fInitialDyDx[i] = DyDx[i];
291  }
292 
293  fLastStepLength = Step;
294 
295  return ;
296 }
297 
298 
299 // Calculate DistChord given start, mid and end-point of step
301 {
302  G4double distLine, distChord;
303  G4ThreeVector initialPoint, finalPoint, midPoint;
304 
305  initialPoint = G4ThreeVector( yStart[0], yStart[1], yStart[2]);
306  finalPoint = G4ThreeVector( yEnd[0], yEnd[1], yEnd[2]);
307  midPoint = G4ThreeVector( yMid[0], yMid[1], yMid[2]);
308 
309  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
310  // distance of Chord
311  if (initialPoint != finalPoint)
312  {
313  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
314  distChord = distLine;
315  }
316  else
317  {
318  distChord = (midPoint-initialPoint).mag();
319  }
320  return distChord;
321 }
322 
323 // (New) DistChord function using interpolation
325 {
326  // Copy the values of stages from this (original) into the Aux Stepper
327  *fAuxStepper = *this;
328 
329  //Preparing for the interpolation
330  fAuxStepper->SetupInterpolation(); // (fLastInitialVector, fInitialDyDx, fLastStepLength);
331  //Interpolate to half step
332  fAuxStepper->Interpolate( /*fLastInitialVector, fInitialDyDx, fLastStepLength,*/ 0.5, fAuxStepper->fMidVector);
333 
335 }
336 
338 {
339  //Coefficients for halfway interpolation
340  const G4double
341  hf1 = 5783653.0/57600000.0 ,
342  hf2 = 0. ,
343  hf3 = 466123.0/1192500.0 ,
344  hf4 = -41347.0/1920000.0 ,
345  hf5 = 16122321.0/339200000.0 ,
346  hf6 = -7117.0/20000.0,
347  hf7 = 183.0/10000.0 ;
348 
349  for(int i=0; i<3; i++){
351  hf1*fInitialDyDx[i] + hf2*ak2[i] + hf3*ak3[i] + hf4*ak4[i] +
352  hf5*ak5[i] + hf6*ak6[i] + hf7*ak7[i] );
353  }
354 
355  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
356  // distance of Chord
357 
359 }
360 
361 //The original DistChord() function for the class
363 {
364  // Do half a step using StepNoErr
368 }
369 
370 // The lower (4th) order interpolant given by Dormand and prince
371 // "An RK 5(4) triple"
372 //---Ref---
373 // J. R. Dormand and P. J. Prince, “Runge-Kutta triples,”
374 // Computers & Mathematics with Applications, vol. 12, no. 9,
375 // pp. 1007–1017, 1986.
376 //---------------------------
377 
378 void G4DormandPrince745::SetupInterpolation_low() // const G4double *yInput, const G4double *dydx, const G4double Step)
379 {
380  //Nothing to be done
381 }
382 
383 void G4DormandPrince745::Interpolate_low( /* const G4double yInput[],
384  const G4double dydx[],
385  const G4double Step, */
386  G4double yOut[],
387  G4double tau )
388 {
389  G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
390  // Coefficients for all the seven stages.
392  const G4double *dydx= fInitialDyDx;
393 
394  const G4int numberOfVariables= this->GetNumberOfVariables();
395 
396  // for(int i=0;i<numberOfVariables;i++) { yIn[i]=yInput[i]; }
397 
398  const G4double
399  tau_2 = tau * tau,
400  tau_3 = tau * tau_2,
401  tau_4 = tau_2 * tau_2;
402 
403  bf1 = (157015080.0*tau_4 - 13107642775.0*tau_3+ 34969693132.0*tau_2- 32272833064.0*tau
404  + 11282082432.0)/11282082432.0,
405  bf2 = 0.0 ,
406  bf3 = - 100.0*tau*(15701508.0*tau_3 - 914128567.0*tau_2 + 2074956840.0*tau
407  - 1323431896.0)/32700410799.0,
408  bf4 = 25.0*tau*(94209048.0*tau_3- 1518414297.0*tau_2+ 2460397220.0*tau - 889289856.0)/5641041216.0 ,
409  bf5 = -2187.0*tau*(52338360.0*tau_3 - 451824525.0*tau_2 + 687873124.0*tau - 259006536.0)/199316789632.0 ,
410  bf6 = 11.0*tau*(106151040.0*tau_3- 661884105.0*tau_2 + 946554244.0*tau - 361440756.0)/2467955532.0 ,
411  bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 82437520.0*tau + 44764047.0)/ 29380423.0 ;
412 
413  //Putting together the coefficients calculated as the respective stage coefficients
414  for( int i=0; i<numberOfVariables; i++){
415  yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] + bf4*ak4[i]
416  + bf5*ak5[i] + bf6*ak6[i] + bf7*ak7[i] ) ;
417  }
418 }
419 
420 
421 // Following interpolant of order 5 was given by Baker,Dormand,Gilmore, Prince :
422 //---Ref---
423 // T. S. Baker, J. R. Dormand, J. P. Gilmore, and P. J. Prince,
424 // “Continuous approximation with embedded Runge-Kutta methods,”
425 // Applied Numerical Mathematics, vol. 22, no. 1, pp. 51–62, 1996.
426 //---------------------
427 
428 // Calculating the extra stages for the interpolant :
429 void G4DormandPrince745::SetupInterpolation_high( /* const G4double yInput[],
430  const G4double dydx[],
431  const G4double Step */ ){
432 
433  //Coefficients for the additional stages :
434  const G4double
435  b81 = 6245.0/62208.0 ,
436  b82 = 0.0 ,
437  b83 = 8875.0/103032.0 ,
438  b84 = -125.0/1728.0 ,
439  b85 = 801.0/13568.0 ,
440  b86 = -13519.0/368064.0 ,
441  b87 = 11105.0/368064.0 ,
442 
443  b91 = 632855.0/4478976.0 ,
444  b92 = 0.0 ,
445  b93 = 4146875.0/6491016.0 ,
446  b94 = 5490625.0/14183424.0 ,
447  b95 = -15975.0/108544.0 ,
448  b96 = 8295925.0/220286304.0 ,
449  b97 = -1779595.0/62938944.0 ,
450  b98 = -805.0/4104.0 ;
451 
452  const G4int numberOfVariables= this->GetNumberOfVariables();
453  const G4double *dydx = fInitialDyDx;
454  const G4double Step = fLastStepLength;
455 
456  // Saving yInput because yInput and yOut can be aliases for same array
457  // for(int i=0;i<numberOfVariables;i++) { yIn[i]=yInput[i]; }
458  // yTemp[7] = yIn[7];
459 
460  //Evaluate the extra stages :
461  for(int i=0;i<numberOfVariables;i++)
462  {
463  yTemp[i] = yIn[i] + Step*(b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
464  b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
465  b87*ak7[i]);
466  }
467  RightHandSide(yTemp, ak8); //8th Stage
468 
469  for(int i=0;i<numberOfVariables;i++)
470  {
471  yTemp[i] = yIn[i] + Step*(b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
472  b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
473  b97*ak7[i] + b98*ak8[i] );
474  }
475  RightHandSide(yTemp, ak9); //9th Stage
476 }
477 
478 
479 // Calculating the interpolated result yOut with the coefficients
480 void G4DormandPrince745::Interpolate_high( /* const G4double yInput[],
481  const G4double dydx[],
482  const G4double Step, */
483  G4double yOut[],
484  G4double tau ){
485  //Define the coefficients for the polynomials
486  G4double bi[10][5], b[10];
487  const G4int numberOfVariables = this->GetNumberOfVariables();
488  const G4double *dydx = fInitialDyDx;
489  // const G4double fullStep = fLastStepLength;
490 
491  // If given requestedStep in argument:
492  // G4double tau = requestedStep / fLastStepLength;
493 
494  // COEFFICIENTS OF bi[1]
495  bi[1][0] = 1.0 ,
496  bi[1][1] = -38039.0/7040.0 ,
497  bi[1][2] = 125923.0/10560.0 ,
498  bi[1][3] = -19683.0/1760.0 ,
499  bi[1][4] = 3303.0/880.0 ,
500  // --------------------------------------------------------
501  //
502  // COEFFICIENTS OF bi[2]
503  bi[2][0] = 0.0 ,
504  bi[2][1] = 0.0 ,
505  bi[2][2] = 0.0 ,
506  bi[2][3] = 0.0 ,
507  bi[2][4] = 0.0 ,
508  // --------------------------------------------------------
509  //
510  // COEFFICIENTS OF bi[3]
511  bi[3][0] = 0.0 ,
512  bi[3][1] = -12500.0/4081.0 ,
513  bi[3][2] = 205000.0/12243.0 ,
514  bi[3][3] = -90000.0/4081.0 ,
515  bi[3][4] = 36000.0/4081.0 ,
516  // --------------------------------------------------------
517  //
518  // COEFFICIENTS OF bi[4]
519  bi[4][0] = 0.0 ,
520  bi[4][1] = -3125.0/704.0 ,
521  bi[4][2] = 25625.0/1056.0 ,
522  bi[4][3] = -5625.0/176.0 ,
523  bi[4][4] = 1125.0/88.0 ,
524  // --------------------------------------------------------
525  //
526  // COEFFICIENTS OF bi[5]
527  bi[5][0] = 0.0 ,
528  bi[5][1] = 164025.0/74624.0 ,
529  bi[5][2] = -448335.0/37312.0 ,
530  bi[5][3] = 295245.0/18656.0 ,
531  bi[5][4] = -59049.0/9328.0 ,
532  // --------------------------------------------------------
533  //
534  // COEFFICIENTS OF bi[6]
535  bi[6][0] = 0.0 ,
536  bi[6][1] = -25.0/28.0 ,
537  bi[6][2] = 205.0/42.0 ,
538  bi[6][3] = -45.0/7.0 ,
539  bi[6][4] = 18.0/7.0 ,
540  // --------------------------------------------------------
541  //
542  // COEFFICIENTS OF bi[7]
543  bi[7][0] = 0.0 ,
544  bi[7][1] = -2.0/11.0 ,
545  bi[7][2] = 73.0/55.0 ,
546  bi[7][3] = -171.0/55.0 ,
547  bi[7][4] = 108.0/55.0 ,
548  // --------------------------------------------------------
549  //
550  // COEFFICIENTS OF bi[8]
551  bi[8][0] = 0.0 ,
552  bi[8][1] = 189.0/22.0 ,
553  bi[8][2] = -1593.0/55.0 ,
554  bi[8][3] = 3537.0/110.0 ,
555  bi[8][4] = -648.0/55.0 ,
556  // --------------------------------------------------------
557  //
558  // COEFFICIENTS OF bi[9]
559  bi[9][0] = 0.0 ,
560  bi[9][1] = 351.0/110.0 ,
561  bi[9][2] = -999.0/55.0 ,
562  bi[9][3] = 2943.0/110.0 ,
563  bi[9][4] = -648.0/55.0 ;
564  // --------------------------------------------------------
565 
566  // for(G4int i = 0; i< numberOfVariables; i++) { yIn[i] = yInput[i]; }
567 
568  // Calculating the polynomials :
569 #if 1
570  for(int iStage=1; iStage<=9; iStage++){
571  b[iStage] = 0;
572  }
573 
574  for(int j=0; j<=4; j++){
575  G4double tauPower = 1.0;
576  for(int iStage=1; iStage<=9; iStage++){
577  b[iStage] += bi[iStage][j]*tauPower;
578  }
579  tauPower *= tau;
580  }
581 #else
582  G4double tau0 = tau;
583 
584  for(int i=1; i<=9; i++){ //Here i is NOT the coordinate no. , it's stage no.
585  b[i] = 0;
586  tau = 1.0;
587  for(int j=0; j<=4; j++){
588  b[i] += bi[i][j]*tau;
589  tau*=tau0;
590  }
591  }
592 #endif
593 
594  G4double stepLen = fLastStepLength * tau;
595  for(int i=0; i<numberOfVariables; i++){ //Here i IS the cooridnate no.
596  yOut[i] = yIn[i] + stepLen *(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
597  b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
598  b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] );
599  }
600 
601 }
602 
603 //G4DormandPrince745::G4DormandPrince745(G4DormandPrince745& DP_Obj){
604 //
605 //}
606 
607 // Overloaded = operator
609 {
610 // this->G4DormandPrince745(right.GetEquationOfMotion(),right.GetNumberOfVariables(), false);
611 
612  int noVars = right.GetNumberOfVariables();
613  for(int i =0; i< noVars; i++)
614  {
615  this->ak2[i] = right.ak2[i];
616  this->ak3[i] = right.ak3[i];
617  this->ak4[i] = right.ak4[i];
618  this->ak5[i] = right.ak5[i];
619  this->ak6[i] = right.ak6[i];
620  this->ak7[i] = right.ak7[i];
621  this->ak8[i] = right.ak8[i];
622  this->ak9[i] = right.ak9[i];
623 
624  this->fInitialDyDx[i] = right.fInitialDyDx[i];
625  this->fLastInitialVector[i] = right.fLastInitialVector[i];
626  this->fMidVector[i] = right.fMidVector[i];
627  this->fMidError[i] = right.fMidError[i];
628  }
629 
630  this->fLastStepLength = right.fLastStepLength;
631 
632  return *this;
633 }
T max(const T t1, const T t2)
brief Return the largest of the two arguments
CLHEP::Hep3Vector G4ThreeVector
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
G4DormandPrince745 * fAuxStepper
G4double DistChord2() const
G4double DistChord3() const
G4DormandPrince745 & operator=(const G4DormandPrince745 &)
void Interpolate_low(G4double yOut[], G4double tau)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4DormandPrince745(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
G4int GetNumberOfStateVariables() const
void RightHandSide(const double y[], double dydx[]) const
void Interpolate_high(G4double yOut[], G4double tau)
G4double DistChord() const
G4double DistLine(G4double yStart[], G4double yMid[], G4double yEnd[]) const
int G4int
Definition: G4Types.hh:78
void Interpolate(G4double tau, G4double yOut[])
G4int GetNumberOfVariables() const
Definition: Step.hh:41