72 G4int noIntegrationVariables,
77 const G4int numberOfVariables =
79 ( ( (noIntegrationVariables-1)/4 + 1 ) * 4 ) );
96 const G4int numStateVars =
194 b31 = 3.0/40.0, b32 = 9.0/40.0 ,
196 b41 = 44.0/45.0, b42 = -56.0/15.0, b43 = 32.0/9.0,
198 b51 = 19372.0/6561.0, b52 = -25360.0/2187.0, b53 = 64448.0/6561.0,
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 ,
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,
221 dc1 = -( b71 - 5179.0/57600.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 );
232 yOut[7] =
yTemp[7] = yInput[7];
235 for(i=0;i<numberOfVariables;i++)
243 for(i=0;i<numberOfVariables;i++)
245 yTemp[i] =
yIn[i] + b21*Step*DyDx[i] ;
249 for(i=0;i<numberOfVariables;i++)
251 yTemp[i] =
yIn[i] + Step*(b31*DyDx[i] + b32*
ak2[i]) ;
255 for(i=0;i<numberOfVariables;i++)
257 yTemp[i] =
yIn[i] + Step*(b41*DyDx[i] + b42*
ak2[i] + b43*
ak3[i]) ;
261 for(i=0;i<numberOfVariables;i++)
268 for(i=0;i<numberOfVariables;i++)
271 b64*
ak4[i] + b65*
ak5[i]) ;
275 for(i=0;i<numberOfVariables;i++)
277 yOut[i] =
yIn[i] + Step*(b71*DyDx[i] + b72*
ak2[i] + b73*
ak3[i] +
282 for(i=0;i<numberOfVariables;i++)
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.5
e-18 ;
305 initialPoint =
G4ThreeVector( yStart[0], yStart[1], yStart[2]);
311 if (initialPoint != finalPoint)
314 distChord = distLine;
318 distChord = (midPoint-initialPoint).mag();
341 hf1 = 5783653.0/57600000.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 ;
349 for(
int i=0; i<3; i++){
389 G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
401 tau_4 = tau_2 * tau_2;
403 bf1 = (157015080.0*tau_4 - 13107642775.0*tau_3+ 34969693132.0*tau_2- 32272833064.0*tau
404 + 11282082432.0)/11282082432.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 ;
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] ) ;
435 b81 = 6245.0/62208.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 ,
443 b91 = 632855.0/4478976.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 ;
461 for(
int i=0;i<numberOfVariables;i++)
469 for(
int i=0;i<numberOfVariables;i++)
473 b97*
ak7[i] + b98*
ak8[i] );
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 ,
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 ,
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 ,
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 ,
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 ,
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 ,
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 ,
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 ;
570 for(
int iStage=1; iStage<=9; iStage++){
574 for(
int j=0; j<=4; j++){
576 for(
int iStage=1; iStage<=9; iStage++){
577 b[iStage] += bi[iStage][j]*tauPower;
584 for(
int i=1; i<=9; i++){
587 for(
int j=0; j<=4; j++){
588 b[i] += bi[i][j]*tau;
595 for(
int i=0; i<numberOfVariables; i++){
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] );
613 for(
int i =0; i< noVars; i++)
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];
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)
void SetupInterpolation_low()
G4double * fLastInitialVector
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 * fLastFinalVector
void SetupInterpolation()
G4double DistLine(G4double yStart[], G4double yMid[], G4double yEnd[]) const
void SetupInterpolation_high()
void Interpolate(G4double tau, G4double yOut[])
G4int GetNumberOfVariables() const