60 G4int noIntegrationVariables,
64 const G4int numberOfVariables = noIntegrationVariables;
69 ak2 =
new G4double[numberOfVariables];
70 ak3 =
new G4double[numberOfVariables];
71 ak4 =
new G4double[numberOfVariables];
72 ak5 =
new G4double[numberOfVariables];
73 ak6 =
new G4double[numberOfVariables];
74 ak7 =
new G4double[numberOfVariables];
76 yTemp =
new G4double[numberOfVariables] ;
77 yIn =
new G4double[numberOfVariables] ;
79 pseudoDydx_for_DistChord =
new G4double[numberOfVariables];
81 fLastInitialVector =
new G4double[numberOfVariables] ;
82 fLastFinalVector =
new G4double[numberOfVariables] ;
83 fLastDyDx =
new G4double[numberOfVariables];
85 fMidVector =
new G4double[numberOfVariables];
86 fMidError =
new G4double[numberOfVariables];
108 delete[] fLastInitialVector;
109 delete[] fLastFinalVector;
116 delete[] pseudoDydx_for_DistChord;
140 b31 = 3.0/40.0, b32 = 9.0/40.0 ,
142 b41 = 44.0/45.0, b42 = -56.0/15.0, b43 = 32.0/9.0,
144 b51 = 19372.0/6561.0, b52 = -25360.0/2187.0, b53 = 64448.0/6561.0,
147 b61 = 9017.0/3168.0 , b62 = -355.0/33.0,
148 b63 = 46732.0/5247.0 , b64 = 49.0/176.0 ,
149 b65 = -5103.0/18656.0 ,
151 b71 = 35.0/384.0, b72 = 0.,
152 b73 = 500.0/1113.0, b74 = 125.0/192.0,
153 b75 = -2187.0/6784.0, b76 = 11.0/84.0,
160 dc1 = b71 - 5179.0/57600.0,
162 dc3 = b73 - 7571.0/16695.0,
163 dc4 = b74 - 393.0/640.0,
164 dc5 = b75 + 92097.0/339200.0,
165 dc6 = b76 - 187.0/2100.0,
173 yOut[7] = yTemp[7] = yIn[7];
176 for(i=0;i<numberOfVariables;i++)
187 for(i=0;i<numberOfVariables;i++)
189 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
193 for(i=0;i<numberOfVariables;i++)
195 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
199 for(i=0;i<numberOfVariables;i++)
201 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
205 for(i=0;i<numberOfVariables;i++)
207 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
212 for(i=0;i<numberOfVariables;i++)
214 yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
215 b64*ak4[i] + b65*ak5[i]) ;
219 for(i=0;i<numberOfVariables;i++)
221 yOut[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
222 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
226 for(i=0;i<numberOfVariables;i++)
229 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
230 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ;
234 fLastInitialVector[i] = yIn[i] ;
235 fLastFinalVector[i] = yOut[i];
236 fLastDyDx[i] = DyDx[i];
237 nextDydx[i] = ak7[i];
242 fLastStepLength = Step;
258 fLastInitialVector[1], fLastInitialVector[2]);
260 fLastFinalVector[1], fLastFinalVector[2]);
264 fAuxStepper->
Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
265 fMidVector, fMidError, pseudoDydx_for_DistChord );
267 midPoint =
G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
273 if (initialPoint != finalPoint)
276 distChord = distLine;
280 distChord = (midPoint-initialPoint).mag();
293 bf1, bf2, bf3, bf4, bf5, bf6, bf7;
303 for(
int i=0;i<numberOfVariables;i++)
313 bf1 = (157015080.0*tau_4 - 13107642775.0*tau_3+ 34969693132.0*tau_2- 32272833064.0*tau
314 + 11282082432.0)/11282082432.0,
316 bf3 = - 100.0*tau*(15701508.0*tau_3 - 914128567.0*tau_2 + 2074956840.0*tau
317 - 1323431896.0)/32700410799.0,
318 bf4 = 25.0*tau*(94209048.0*tau_3- 1518414297.0*tau_2+ 2460397220.0*tau - 889289856.0)/5641041216.0 ,
319 bf5 = -2187.0*tau*(52338360.0*tau_3 - 451824525.0*tau_2 + 687873124.0*tau - 259006536.0)/199316789632.0 ,
320 bf6 = 11.0*tau*(106151040.0*tau_3- 661884105.0*tau_2 + 946554244.0*tau - 361440756.0)/2467955532.0 ,
321 bf7 = tau*(1.0 - tau)*(8293050.0*tau_2 - 82437520.0*tau + 44764047.0)/ 29380423.0 ;
324 for(
int i=0; i<numberOfVariables; i++){
325 yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i] + bf4*ak4[i]
326 + bf5*ak5[i] + bf6*ak6[i] + bf7*ak7[i] ) ;
339 b81 = 6245.0/62208.0 ,
341 b83 = 8875.0/103032.0 ,
342 b84 = -125.0/1728.0 ,
343 b85 = 801.0/13568.0 ,
344 b86 = -13519.0/368064.0 ,
345 b87 = 11105.0/368064.0 ,
347 b91 = 632855.0/4478976.0 ,
349 b93 = 4146875.0/6491016.0 ,
350 b94 = 5490625.0/14183424.0 ,
351 b95 = -15975.0/108544.0 ,
352 b96 = 8295925.0/220286304.0 ,
353 b97 = -1779595.0/62938944.0 ,
354 b98 = -805.0/4104.0 ;
359 for(
int i=0;i<numberOfVariables;i++)
366 ak8 =
new G4double[numberOfVariables];
367 ak9 =
new G4double[numberOfVariables];
370 for(
int i=0;i<numberOfVariables;i++)
372 yTemp[i] = yIn[i] + Step*( b81*dydx[i] + b82*ak2[i] + b83*ak3[i] +
373 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
378 for(
int i=0;i<numberOfVariables;i++)
380 yTemp[i] = yIn[i] + Step * ( b91*dydx[i] + b92*ak2[i] + b93*ak3[i] +
381 b94*ak4[i] + b95*ak5[i] + b96*ak6[i] +
382 b97*ak7[i] + b98*ak8[i] );
403 bi[1][1] = -38039.0/7040.0 ,
404 bi[1][2] = 125923.0/10560.0 ,
405 bi[1][3] = -19683.0/1760.0 ,
406 bi[1][4] = 3303.0/880.0 ,
419 bi[3][1] = -12500.0/4081.0 ,
420 bi[3][2] = 205000.0/12243.0 ,
421 bi[3][3] = -90000.0/4081.0 ,
422 bi[3][4] = 36000.0/4081.0 ,
427 bi[4][1] = -3125.0/704.0 ,
428 bi[4][2] = 25625.0/1056.0 ,
429 bi[4][3] = -5625.0/176.0 ,
430 bi[4][4] = 1125.0/88.0 ,
435 bi[5][1] = 164025.0/74624.0 ,
436 bi[5][2] = -448335.0/37312.0 ,
437 bi[5][3] = 295245.0/18656.0 ,
438 bi[5][4] = -59049.0/9328.0 ,
443 bi[6][1] = -25.0/28.0 ,
444 bi[6][2] = 205.0/42.0 ,
445 bi[6][3] = -45.0/7.0 ,
446 bi[6][4] = 18.0/7.0 ,
451 bi[7][1] = -2.0/11.0 ,
452 bi[7][2] = 73.0/55.0 ,
453 bi[7][3] = -171.0/55.0 ,
454 bi[7][4] = 108.0/55.0 ,
459 bi[8][1] = 189.0/22.0 ,
460 bi[8][2] = -1593.0/55.0 ,
461 bi[8][3] = 3537.0/110.0 ,
462 bi[8][4] = -648.0/55.0 ,
467 bi[9][1] = 351.0/110.0 ,
468 bi[9][2] = -999.0/55.0 ,
469 bi[9][3] = 2943.0/110.0 ,
470 bi[9][4] = -648.0/55.0 ;
475 for(
G4int i = 0; i< numberOfVariables; i++)
481 for(
int i=1; i<=9; i++){
484 for(
int j=0; j<=4; j++){
485 b[i] += bi[i][j]*tau;
490 for(
int i=0; i<numberOfVariables; i++){
491 yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
492 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
493 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] );
CLHEP::Hep3Vector G4ThreeVector
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[], G4double nextDydx[])
G4double DistChord() const
void RightHandSide(const double y[], double dydx[])
void SetupInterpolate(const G4double yInput[], const G4double dydx[], const G4double Step)
G4FSALDormandPrince745(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
void Interpolate(const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
~G4FSALDormandPrince745()
G4int GetNumberOfVariables() const
void interpolate(const G4double yInput[], const G4double dydx[], G4double yOut[], G4double Step, G4double tau)