| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4BogackiShampine23 implementation << 26 // Bogacki-Shampine - 4 - 3(2) non-FSAL implementation by Somnath Banerjee
>> 27 // Supervision / code review: John Apostolakis
27 // 28 //
28 // Bogacki-Shampine - 4 - 3(2) non-FSAL imple << 29 // Somnath's work was sponsored by Google as part of the Google Summer of
>> 30 // Code 2015, as part of the CERN / SFT organisation.p
>> 31 // ===================================================================
29 // 32 //
30 // Implementation of the method proposed in t 33 // Implementation of the method proposed in the publication
31 // "A 3(2) pair of Runge - Kutta formulas" << 34 // “A 3(2) pair of Runge - Kutta formulas,”
32 // by P. Bogacki and L. F. Shampine, << 35 // by P. Bogacki and L. F. Shampine,
33 // Appl. Math. Lett., vol. 2, no. 4, pp. 32 << 36 // Appl. Math. Lett., vol. 2, no. 4, pp. 321–325, Jan. 1989.
34 // 37 //
35 // The Bogacki shampine method has the followi << 38 // First version: 20 May 2015
36 // 39 //
37 // 0 | << 40 // History
38 // 1/2|1/2 << 41 // -----------------------------
39 // 3/4|0 3/4 << 42 // Created by Somnath Banerjee on 20 May 2015
40 // 1 |2/9 1/3 4/9 << 43 ///////////////////////////////////////////////////////////////////////////////
41 // ------------------- << 44
42 // |2/9 1/3 4/9 0 << 45
43 // |7/24 1/4 1/3 1/8 << 46 /*
44 // << 47
45 // Created: Somnath Banerjee, Google Summer of << 48 This contains the stepper function of the G4BogackiShampine23 class
46 // Supervision: John Apostolakis, CERN << 49
47 // ------------------------------------------- << 50 The Bogacki shampine method has the following Butcher's tableau
>> 51
>> 52 0 |
>> 53 1/2|1/2
>> 54 3/4|0 3/4
>> 55 1 |2/9 1/3 4/9
>> 56 -------------------
>> 57 |2/9 1/3 4/9 0
>> 58 |7/24 1/4 1/3 1/8
>> 59
>> 60 */
48 61
49 #include "G4BogackiShampine23.hh" 62 #include "G4BogackiShampine23.hh"
50 #include "G4LineSection.hh" 63 #include "G4LineSection.hh"
51 #include "G4FieldUtils.hh" 64 #include "G4FieldUtils.hh"
52 65
53 using namespace field_utils; 66 using namespace field_utils;
54 67
55 G4BogackiShampine23::G4BogackiShampine23(G4Equ 68 G4BogackiShampine23::G4BogackiShampine23(G4EquationOfMotion* EqRhs,
56 G4int << 69 G4int integrationVariables):
57 : G4MagIntegratorStepper(EqRhs, integrationV << 70 G4MagIntegratorStepper(EqRhs, integrationVariables)
58 { 71 {
>> 72
59 SetIntegrationOrder(3); 73 SetIntegrationOrder(3);
60 SetFSAL(true); 74 SetFSAL(true);
61 } 75 }
62 76
63 void G4BogackiShampine23::makeStep(const G4dou 77 void G4BogackiShampine23::makeStep(const G4double yInput[],
64 const G4dou 78 const G4double dydx[],
65 const G4dou 79 const G4double hstep,
66 G4double yO 80 G4double yOutput[],
67 G4double* d 81 G4double* dydxOutput,
68 G4double* y 82 G4double* yError) const
69 { 83 {
70 84
71 G4double yTemp[G4FieldTrack::ncompSVEC]; 85 G4double yTemp[G4FieldTrack::ncompSVEC];
72 for(G4int i = GetNumberOfVariables(); i < Ge << 86 for(G4int i = GetNumberOfVariables(); i < GetNumberOfStateVariables(); ++i)
73 { <<
74 yOutput[i] = yTemp[i] = yInput[i]; 87 yOutput[i] = yTemp[i] = yInput[i];
75 } << 88
76 89
77 G4double ak2[G4FieldTrack::ncompSVEC], 90 G4double ak2[G4FieldTrack::ncompSVEC],
78 ak3[G4FieldTrack::ncompSVEC]; 91 ak3[G4FieldTrack::ncompSVEC];
79 92
80 const G4double b21 = 0.5 , << 93 const G4double b21 = 0.5 ,
81 b31 = 0., b32 = 3.0 / 4.0, << 94 b31 = 0., b32 = 3.0 / 4.0,
82 b41 = 2.0 / 9.0, b42 = 1.0 / << 95 b41 = 2.0 / 9.0, b42 = 1.0 / 3.0, b43 = 4.0 / 9.0;
83 96
84 const G4double dc1 = b41 - 7.0 / 24.0, dc2 << 97 const G4double dc1 = b41 - 7.0 / 24.0, dc2 = b42 - 1.0 / 4.0,
85 dc3 = b43 - 1.0 / 3.0, dc4 << 98 dc3 = b43 - 1.0 / 3.0, dc4 = - 1.0 / 8.0;
86 99
87 // RightHandSide(yInput, dydx); << 100 // RightHandSide(yInput, dydx);
88 for(G4int i = 0; i < GetNumberOfVariables(); << 101 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
89 { << 102 yTemp[i] = yInput[i] + b21 * hstep * dydx[i];
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 103
105 if ((dydxOutput != nullptr) && (yError != nu << 104 RightHandSide(yTemp, ak2);
106 { <<
107 RightHandSide(yOutput, dydxOutput); <<
108 for(G4int i = 0; i < GetNumberOfVariables( 105 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
109 { << 106 yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
110 yError[i] = hstep * (dc1 * dydx[i] + dc2 << 107
111 dc3 * ak3[i] + dc4 << 108 RightHandSide(yTemp, ak3);
112 } << 109 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
>> 110 yOutput[i] = yInput[i] + hstep * (b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
>> 111
>> 112 if (dydxOutput && yError) {
>> 113 RightHandSide(yOutput, dydxOutput);
>> 114 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
>> 115 yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] +
>> 116 dc3 * ak3[i] + dc4 * dydxOutput[i]);
>> 117
113 } 118 }
114 } 119 }
115 120
116 void G4BogackiShampine23::Stepper(const G4doub 121 void G4BogackiShampine23::Stepper(const G4double yInput[],
117 const G4doub 122 const G4double dydx[],
118 G4double hst 123 G4double hstep,
119 G4double yOu 124 G4double yOutput[],
120 G4double yEr 125 G4double yError[])
121 { 126 {
122 copy(fyIn, yInput); 127 copy(fyIn, yInput);
123 copy(fdydx, dydx); 128 copy(fdydx, dydx);
124 fhstep = hstep; 129 fhstep = hstep;
125 130
126 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOu 131 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
127 132
128 copy(yOutput, fyOut); 133 copy(yOutput, fyOut);
129 } 134 }
130 135
131 void G4BogackiShampine23::Stepper(const G4doub 136 void G4BogackiShampine23::Stepper(const G4double yInput[],
132 const G4doub 137 const G4double dydx[],
133 G4double hst 138 G4double hstep,
134 G4double yOu 139 G4double yOutput[],
135 G4double yEr 140 G4double yError[],
136 G4double dyd 141 G4double dydxOutput[])
137 { 142 {
138 copy(fyIn, yInput); 143 copy(fyIn, yInput);
139 copy(fdydx, dydx); 144 copy(fdydx, dydx);
140 fhstep = hstep; 145 fhstep = hstep;
141 146
142 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOu 147 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
143 148
144 copy(yOutput, fyOut); 149 copy(yOutput, fyOut);
145 copy(dydxOutput, fdydxOut); 150 copy(dydxOutput, fdydxOut);
146 } 151 }
147 152
148 G4double G4BogackiShampine23::DistChord() cons 153 G4double G4BogackiShampine23::DistChord() const
149 { 154 {
150 G4double yMid[G4FieldTrack::ncompSVEC]; 155 G4double yMid[G4FieldTrack::ncompSVEC];
151 makeStep(fyIn, fdydx, fhstep / 2., yMid); 156 makeStep(fyIn, fdydx, fhstep / 2., yMid);
152 157
153 const G4ThreeVector begin = makeVector(fyIn, 158 const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
154 const G4ThreeVector mid = makeVector(yMid, V 159 const G4ThreeVector mid = makeVector(yMid, Value3D::Position);
155 const G4ThreeVector end = makeVector(fyOut, 160 const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
156 161
157 return G4LineSection::Distline(mid, begin, e 162 return G4LineSection::Distline(mid, begin, end);
158 } 163 }
159 164