| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4BogackiShampine23 implementation
27 //
28 // Bogacki-Shampine - 4 - 3(2) non-FSAL imple
29 //
30 // Implementation of the method proposed in t
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. 32
34 //
35 // The Bogacki shampine method has the followi
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
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
55 G4BogackiShampine23::G4BogackiShampine23(G4Equ
56 G4int
57 : G4MagIntegratorStepper(EqRhs, integrationV
58 {
59 SetIntegrationOrder(3);
60 SetFSAL(true);
61 }
62
63 void G4BogackiShampine23::makeStep(const G4dou
64 const G4dou
65 const G4dou
66 G4double yO
67 G4double* d
68 G4double* y
69 {
70
71 G4double yTemp[G4FieldTrack::ncompSVEC];
72 for(G4int i = GetNumberOfVariables(); i < Ge
73 {
74 yOutput[i] = yTemp[i] = yInput[i];
75 }
76
77 G4double ak2[G4FieldTrack::ncompSVEC],
78 ak3[G4FieldTrack::ncompSVEC];
79
80 const G4double b21 = 0.5 ,
81 b31 = 0., b32 = 3.0 / 4.0,
82 b41 = 2.0 / 9.0, b42 = 1.0 /
83
84 const G4double dc1 = b41 - 7.0 / 24.0, dc2
85 dc3 = b43 - 1.0 / 3.0, dc4
86
87 // RightHandSide(yInput, dydx);
88 for(G4int i = 0; i < GetNumberOfVariables();
89 {
90 yTemp[i] = yInput[i] + b21 * hstep * dydx[
91 }
92
93 RightHandSide(yTemp, ak2);
94 for(G4int i = 0; i < GetNumberOfVariables();
95 {
96 yTemp[i] = yInput[i] + hstep * (b31 * dydx
97 }
98
99 RightHandSide(yTemp, ak3);
100 for(G4int i = 0; i < GetNumberOfVariables();
101 {
102 yOutput[i] = yInput[i] + hstep * (b41*dydx
103 }
104
105 if ((dydxOutput != nullptr) && (yError != nu
106 {
107 RightHandSide(yOutput, dydxOutput);
108 for(G4int i = 0; i < GetNumberOfVariables(
109 {
110 yError[i] = hstep * (dc1 * dydx[i] + dc2
111 dc3 * ak3[i] + dc4
112 }
113 }
114 }
115
116 void G4BogackiShampine23::Stepper(const G4doub
117 const G4doub
118 G4double hst
119 G4double yOu
120 G4double yEr
121 {
122 copy(fyIn, yInput);
123 copy(fdydx, dydx);
124 fhstep = hstep;
125
126 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOu
127
128 copy(yOutput, fyOut);
129 }
130
131 void G4BogackiShampine23::Stepper(const G4doub
132 const G4doub
133 G4double hst
134 G4double yOu
135 G4double yEr
136 G4double dyd
137 {
138 copy(fyIn, yInput);
139 copy(fdydx, dydx);
140 fhstep = hstep;
141
142 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOu
143
144 copy(yOutput, fyOut);
145 copy(dydxOutput, fdydxOut);
146 }
147
148 G4double G4BogackiShampine23::DistChord() cons
149 {
150 G4double yMid[G4FieldTrack::ncompSVEC];
151 makeStep(fyIn, fdydx, fhstep / 2., yMid);
152
153 const G4ThreeVector begin = makeVector(fyIn,
154 const G4ThreeVector mid = makeVector(yMid, V
155 const G4ThreeVector end = makeVector(fyOut,
156
157 return G4LineSection::Distline(mid, begin, e
158 }
159