77 G4int noIntegrationVariables,
80 fLastStepLength(-1.0),
82 fPreparedInterpolation(false)
85 const G4int numberOfVariables = noIntegrationVariables;
101 for (
int i = 0; i < 6; i++) {
144 for (
int i = 0; i < 6; i++) {
169 for(
G4int i=0; i < numberOfVariables; i++ ){
170 dyDxLast[i] =
ak9[i];
190 b31 = 2.0/27.0 , b32 = 4.0/27.0,
192 b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
194 b51 = 68.0/297.0, b52 = -4.0/11.0,
195 b53 = 42.0/143.0, b54 = 1960.0/3861.0,
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,
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,
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;
224 dc1 = b81 - 2479.0 / 34992.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;
239 for(i=0;i<numberOfVariables;i++)
247 for(i=0;i<numberOfVariables;i++)
249 yTemp[i] =
yIn[i] + b21*Step*DyDx[i] ;
253 for(i=0;i<numberOfVariables;i++)
255 yTemp[i] =
yIn[i] + Step*(b31*DyDx[i] + b32*
ak2[i]) ;
259 for(i=0;i<numberOfVariables;i++)
261 yTemp[i] =
yIn[i] + Step*(b41*DyDx[i] + b42*
ak2[i] + b43*
ak3[i]) ;
265 for(i=0;i<numberOfVariables;i++)
272 for(i=0;i<numberOfVariables;i++)
275 b64*
ak4[i] + b65*
ak5[i]) ;
279 for(i=0;i<numberOfVariables;i++)
286 for(i=0;i<numberOfVariables;i++)
288 yOut[i] =
yIn[i] + Step*(b81*DyDx[i] + b82*
ak2[i] + b83*
ak3[i] +
294 for(i=0;i<numberOfVariables;i++)
296 yErr[i] = Step*(dc1*DyDx[i] + dc2*
ak2[i] + dc3*
ak3[i] + dc4*
ak4[i] +
350 if (initialPoint != finalPoint)
353 distChord = distLine;
357 distChord = (midPoint-initialPoint).mag();
369 a93 = 10256301.0/35409920.0 ,
370 a94 = 2307361.0/17971200.0 ,
371 a95 = -387.0/102400.0 ,
373 a97 = -7267.0/215040.0 ,
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 ,
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 ,
395 a1110 = -1403317093.0/11371610250.0 ;
406 for(
int i=0; i<numberOfVariables; i++){
409 a97*
ak7[i] + a98*
ak8[i] );
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] );
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] +
431 int nwant = numberOfVariables;
437 for (
int l = 0; l < nwant; l++) {
439 p[5][l] =
bi[5][6]*
ak5[l] +
448 bi[5][5]*ak5[l]) + ((
bi[4][5]*
ak4[l] +
450 bi[11][5]*
ak11[l]) +
bi[1][5]*dydx[l]);
454 bi[5][4]*ak5[l]) + ((
bi[10][4]*
ak10[l] +
456 bi[11][4]*
ak11[l]) +
bi[1][4]*dydx[l]);
458 p[2][l] =
bi[5][3]*ak5[l] +
bi[6][3]*
ak6[l] +
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]) +
472 for (
int i = 0; i < 6; i++) {
473 for (
int l = 0; l < nwant; l++) {
483 for(
int i=1; i<= 11; i++)
486 for(
int i=1; i<=6; i++)
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 ,
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 ,
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 ,
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 ,
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 ,
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 ,
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 ,
532 bi[9][6] = 55.0 / 3.0 ,
535 bi[9][3] = 199.0 / 3.0 ,
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 ,
555 G4int numberOfVariables = this->GetNumberOfVariables();
556 assert( fPreparedConstants);
558 G4Exception(
"G4BogackiShampine45::InterpolateHigh()",
"GeomField0001",
567 G4int nwant = numberOfVariables;
568 const G4int norder= 6;
571 for (l = 0; l < nwant; l++) {
572 yOut[l] =
p[norder-1][l] * tau;
574 for (k = norder - 2; k >= 1; k--) {
575 for (l = 0; l < nwant; l++) {
576 yOut[l] = ( yOut[l] +
p[k][l] ) * tau;
579 for (l = 0; l < nwant; l++) {
580 yOut[l] = ( yOut[l] + Step * ak8[l] ) * tau + yIn[l];
590 for(
int iStage=1; iStage<=11; iStage++){
593 for(
int j=6; j>=1; j--){
594 b[iStage] += bi[iStage][j] * tau;
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] );
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void SetupInterpolationHigh()
CLHEP::Hep3Vector G4ThreeVector
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
static bool fPreparedConstants
G4double * fLastInitialVector
void InterpolateHigh(G4double tau, G4double yOut[]) const
G4double DistChord() const
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
bool fPreparedInterpolation
G4BogackiShampine45 * fAuxStepper
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)
G4double * fLastFinalVector
static const G4double ak2
G4int GetNumberOfVariables() const