ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RK547FEq2.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RK547FEq2.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 // G4RK547FEq2 implementation
27 //
28 // The Butcher table of the Higham & Hall 5(4)7 method is:
29 //
30 // 0 |
31 // 2/13 | 2/13
32 // 2/13 | 3/52 9/52
33 // 5/9 | 12955/26244 -15925/8748 12350/6561
34 // 3/4 | -10383/52480 13923/10496 -176553/199424 505197/997120
35 // 1 | 1403/7236 -429/268 733330/309339 -7884/8911 104960/113967
36 // 1 | 181/2700 0 656903/1846800 19683/106400 34112/110565 67/800
37 //----------------------------------------------------------------------------------------------------------------------
38 // 181/2700 0 656903/1846800 19683/106400 34112/110565 67/800 0
39 // 11377/154575 0 35378291/105729300 343359/1522850 535952/1947645 134/17175 1/12
40 //
41 // Author: Dmitry Sorokin, Google Summer of Code 2017
42 // Supervision: John Apostolakis, CERN
43 // --------------------------------------------------------------------
44 
45 #include "G4RK547FEq2.hh"
46 #include "G4LineSection.hh"
47 #include "G4FieldUtils.hh"
48 
49 using namespace field_utils;
50 
51 G4RK547FEq2::G4RK547FEq2(G4EquationOfMotion* EqRhs, G4int integrationVariables)
52  : G4MagIntegratorStepper(EqRhs, integrationVariables)
53 {
54 }
55 
56 void G4RK547FEq2::makeStep( const G4double yInput[],
57  const G4double dydx[],
58  const G4double hstep,
59  G4double yOutput[],
60  G4double* dydxOutput,
61  G4double* yError ) const
62 {
63  G4double yTemp[G4FieldTrack::ncompSVEC];
64  for (G4int i=GetNumberOfVariables(); i<GetNumberOfStateVariables(); ++i)
65  {
66  yOutput[i] = yTemp[i] = yInput[i];
67  }
68 
69  G4double ak2[G4FieldTrack::ncompSVEC],
70  ak3[G4FieldTrack::ncompSVEC],
71  ak4[G4FieldTrack::ncompSVEC],
72  ak5[G4FieldTrack::ncompSVEC],
73  ak6[G4FieldTrack::ncompSVEC];
74 
75  const G4double b21 = 2./13.,
76  b31 = 3./52., b32 = 9./52.,
77  b41 = 12955./26244., b42 = -15925./8748.,
78  b43 = 12350./6561.,
79  b51 = -10383./52480., b52 = 13923./10496.,
80  b53 = -176553./199424., b54 = 505197./997120.,
81  b61 = 1403./7236., b62 = -429./268., b63 = 733330./309339.,
82  b64 = -7884./8911., b65 = 104960./113967.,
83  b71 = 181./2700., b72 = 0., b73 = 656903./1846800.,
84  b74 = 19683./106400., b75 = 34112./110565.,
85  b76 = 67./800.;
86 
87  const G4double dc1 = b71 - 11377./154575.,
88  dc2 = b72 - 0.,
89  dc3 = b73 - 35378291./105729300.,
90  dc4 = b74 - 343359./1522850.,
91  dc5 = b75 - 535952./1947645.,
92  dc6 = b76 - 134./17175.,
93  dc7 = 0. - 1./12.;
94 
95  // RightHandSide(yInput, dydx);
96  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
97  yTemp[i] = yInput[i] + hstep * b21 * dydx[i];
98 
99  RightHandSide(yTemp, ak2);
100  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
101  yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
102 
103  RightHandSide(yTemp, ak3);
104  for(G4int i = 0;i < GetNumberOfVariables(); ++i)
105  yTemp[i] = yInput[i] + hstep * (b41 * dydx[i] + b42 * ak2[i] +
106  b43 * ak3[i]);
107 
108  RightHandSide(yTemp, ak4);
109  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
110  yTemp[i] = yInput[i] + hstep * (b51 * dydx[i] + b52 * ak2[i] +
111  b53 * ak3[i] + b54 * ak4[i]);
112 
113  RightHandSide(yTemp, ak5);
114  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
115  yTemp[i] = yInput[i] + hstep * (b61 * dydx[i] + b62 * ak2[i] +
116  b63 * ak3[i] + b64 * ak4[i] +
117  b65 * ak5[i]);
118 
119  RightHandSide(yTemp, ak6);
120  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
121  yOutput[i] = yInput[i] + hstep * (b71 * dydx[i] + b72 * ak2[i] +
122  b73 * ak3[i] + b74 * ak4[i] +
123  b75 * ak5[i] + b76 * ak6[i]);
124  if (dydxOutput && yError)
125  {
126  RightHandSide(yOutput, dydxOutput);
127  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
128  yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] + dc3 * ak3[i] +
129  dc4 * ak4[i] + dc5 * ak5[i] + dc6 * ak6[i] +
130  dc7 * dydxOutput[i]);
131  }
132 }
133 
134 void G4RK547FEq2::Stepper( const G4double yInput[],
135  const G4double dydx[],
136  G4double hstep,
137  G4double yOutput[],
138  G4double yError[] )
139 {
140  copy(fyIn, yInput);
141  copy(fdydx, dydx);
142  fhstep = hstep;
143 
144  makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
145 
146  copy(yOutput, fyOut);
147 }
148 
149 void G4RK547FEq2::Stepper( const G4double yInput[],
150  const G4double dydx[],
151  G4double hstep,
152  G4double yOutput[],
153  G4double yError[],
154  G4double dydxOutput[] )
155 {
156  copy(fyIn, yInput);
157  copy(fdydx, dydx);
158  fhstep = hstep;
159 
160  makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
161 
162  copy(yOutput, fyOut);
163  copy(dydxOutput, fdydxOut);
164 }
165 
166 G4double G4RK547FEq2::DistChord() const
167 {
168  G4double yMid[G4FieldTrack::ncompSVEC];
169  makeStep(fyIn, fdydx, fhstep / 2., yMid);
170 
171  const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
172  const G4ThreeVector mid = makeVector(yMid, Value3D::Position);
173  const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
174 
175  return G4LineSection::Distline(mid, begin, end);
176 }