ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BogackiShampine23.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BogackiShampine23.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 // G4BogackiShampine23 implementation
27 //
28 // Bogacki-Shampine - 4 - 3(2) non-FSAL implementation
29 //
30 // Implementation of the method proposed in the publication
31 // "A 3(2) pair of Runge - Kutta formulas"
32 // by P. Bogacki and L. F. Shampine,
33 // Appl. Math. Lett., vol. 2, no. 4, pp. 321-325, Jan. 1989.
34 //
35 // The Bogacki shampine method has the following Butcher's tableau
36 //
37 // 0 |
38 // 1/2|1/2
39 // 3/4|0 3/4
40 // 1 |2/9 1/3 4/9
41 // -------------------
42 // |2/9 1/3 4/9 0
43 // |7/24 1/4 1/3 1/8
44 //
45 // Created: Somnath Banerjee, Google Summer of Code 2015, 20 May 2015
46 // Supervision: John Apostolakis, CERN
47 // --------------------------------------------------------------------
48 
49 #include "G4BogackiShampine23.hh"
50 #include "G4LineSection.hh"
51 #include "G4FieldUtils.hh"
52 
53 using namespace field_utils;
54 
56  G4int integrationVariables)
57  : G4MagIntegratorStepper(EqRhs, integrationVariables)
58 {
60  SetFSAL(true);
61 }
62 
64  const G4double dydx[],
65  const G4double hstep,
66  G4double yOutput[],
67  G4double* dydxOutput,
68  G4double* yError) const
69 {
70 
73  {
74  yOutput[i] = yTemp[i] = yInput[i];
75  }
76 
79 
80  const G4double b21 = 0.5 ,
81  b31 = 0., b32 = 3.0 / 4.0,
82  b41 = 2.0 / 9.0, b42 = 1.0 / 3.0, b43 = 4.0 / 9.0;
83 
84  const G4double dc1 = b41 - 7.0 / 24.0, dc2 = b42 - 1.0 / 4.0,
85  dc3 = b43 - 1.0 / 3.0, dc4 = - 1.0 / 8.0;
86 
87  // RightHandSide(yInput, dydx);
88  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
89  {
90  yTemp[i] = yInput[i] + b21 * hstep * dydx[i];
91  }
92 
93  RightHandSide(yTemp, ak2);
94  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
95  {
96  yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
97  }
98 
99  RightHandSide(yTemp, ak3);
100  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
101  {
102  yOutput[i] = yInput[i] + hstep * (b41*dydx[i] + b42*ak2[i] + b43*ak3[i]);
103  }
104 
105  if (dydxOutput && yError)
106  {
107  RightHandSide(yOutput, dydxOutput);
108  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
109  {
110  yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] +
111  dc3 * ak3[i] + dc4 * dydxOutput[i]);
112  }
113  }
114 }
115 
117  const G4double dydx[],
118  G4double hstep,
119  G4double yOutput[],
120  G4double yError[])
121 {
122  copy(fyIn, yInput);
123  copy(fdydx, dydx);
124  fhstep = hstep;
125 
126  makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
127 
128  copy(yOutput, fyOut);
129 }
130 
132  const G4double dydx[],
133  G4double hstep,
134  G4double yOutput[],
135  G4double yError[],
136  G4double dydxOutput[])
137 {
138  copy(fyIn, yInput);
139  copy(fdydx, dydx);
140  fhstep = hstep;
141 
142  makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
143 
144  copy(yOutput, fyOut);
145  copy(dydxOutput, fdydxOut);
146 }
147 
149 {
151  makeStep(fyIn, fdydx, fhstep / 2., yMid);
152 
153  const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
154  const G4ThreeVector mid = makeVector(yMid, Value3D::Position);
155  const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
156 
157  return G4LineSection::Distline(mid, begin, end);
158 }