| 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 // G4RK547FEq3 implementation 26 // G4RK547FEq3 implementation
27 // 27 //
28 // The Butcher table of the Higham & Hall 5(4) 28 // The Butcher table of the Higham & Hall 5(4)7 method is:
29 // 29 //
30 // 0 | 30 // 0 |
31 // 11/45 | 11/45 31 // 11/45 | 11/45
32 // 11/30 | 11/120 11/40 32 // 11/30 | 11/120 11/40
33 // 55/56 | 106865/87808 -408375/87808 33 // 55/56 | 106865/87808 -408375/87808 193875/43904
34 // 9/10 | 79503/121000 -1053/440 34 // 9/10 | 79503/121000 -1053/440 147753/56870 27048/710875
35 // 1 | 89303/78045 -2025/473 35 // 1 | 89303/78045 -2025/473 994650/244541 -2547216/28122215 475/2967
36 // 1 | 1247/10890 0 36 // 1 | 1247/10890 0 57375/108053 -1229312/1962015 125/207 43/114
37 //-------------------------------------------- 37 //---------------------------------------------------------------------------------------------------------------------
38 // 1247/10890 0 38 // 1247/10890 0 57375/108053 -1229312/1962015 125/207 43/114 0
39 // 21487/185130 0 39 // 21487/185130 0 963225/1836901 -39864832/33354255 2575/3519 4472/4845 -1/10
40 // 40 //
41 // Author: Dmitry Sorokin, Google Summer of Co 41 // Author: Dmitry Sorokin, Google Summer of Code 2017
42 // Supervision: John Apostolakis, CERN 42 // Supervision: John Apostolakis, CERN
43 // ------------------------------------------- 43 // --------------------------------------------------------------------
44 44
45 #include "G4RK547FEq3.hh" 45 #include "G4RK547FEq3.hh"
46 #include "G4LineSection.hh" 46 #include "G4LineSection.hh"
47 #include "G4FieldUtils.hh" 47 #include "G4FieldUtils.hh"
48 48
49 using namespace field_utils; 49 using namespace field_utils;
50 50
51 G4RK547FEq3::G4RK547FEq3(G4EquationOfMotion* E 51 G4RK547FEq3::G4RK547FEq3(G4EquationOfMotion* EqRhs, G4int integrationVariables)
52 : G4MagIntegratorStepper(EqRhs, integrationV 52 : G4MagIntegratorStepper(EqRhs, integrationVariables)
53 { 53 {
54 } 54 }
55 55
56 void G4RK547FEq3::makeStep( const G4double yIn 56 void G4RK547FEq3::makeStep( const G4double yInput[],
57 const G4double dyd 57 const G4double dydx[],
58 const G4double hst 58 const G4double hstep,
59 G4double yOu 59 G4double yOutput[],
60 G4double* dy 60 G4double* dydxOutput,
61 G4double* yE 61 G4double* yError ) const
62 { 62 {
63 G4double yTemp[G4FieldTrack::ncompSVEC]; 63 G4double yTemp[G4FieldTrack::ncompSVEC];
64 for (G4int i=GetNumberOfVariables(); i<Get 64 for (G4int i=GetNumberOfVariables(); i<GetNumberOfStateVariables(); ++i)
65 { 65 {
66 yOutput[i] = yTemp[i] = yInput[i]; 66 yOutput[i] = yTemp[i] = yInput[i];
67 } 67 }
68 68
69 G4double ak2[G4FieldTrack::ncompSVEC], 69 G4double ak2[G4FieldTrack::ncompSVEC],
70 ak3[G4FieldTrack::ncompSVEC], 70 ak3[G4FieldTrack::ncompSVEC],
71 ak4[G4FieldTrack::ncompSVEC], 71 ak4[G4FieldTrack::ncompSVEC],
72 ak5[G4FieldTrack::ncompSVEC], 72 ak5[G4FieldTrack::ncompSVEC],
73 ak6[G4FieldTrack::ncompSVEC]; 73 ak6[G4FieldTrack::ncompSVEC];
74 74
75 const G4double b21 = 11./45., 75 const G4double b21 = 11./45.,
76 b31 = 11./120., b32 = 11./4 76 b31 = 11./120., b32 = 11./40.,
77 b41 = 106865./87808., b42 = 77 b41 = 106865./87808., b42 = -408375./87808.,
78 b43 = 193875./43904., 78 b43 = 193875./43904.,
79 b51 = 79503./121000., b52 = 79 b51 = 79503./121000., b52 = -1053./440.,
80 b53 = 147753./56870., b 80 b53 = 147753./56870., b54 = 27048./710875.,
81 b61 = 89303./78045., b62 = 81 b61 = 89303./78045., b62 = -2025./473.,
82 b63 = 994650./244541., 82 b63 = 994650./244541., b64 = -2547216./28122215.,
83 b65 = 475./2967., 83 b65 = 475./2967.,
84 b71 = 1247./10890., b72 = 0 84 b71 = 1247./10890., b72 = 0., b73 = 57375./108053.,
85 b74 = -1229312./1962015 85 b74 = -1229312./1962015., b75 = 125./207.,
86 b76 = 43./114.; 86 b76 = 43./114.;
87 87
88 const G4double dc1 = b71 - 21487./185130., 88 const G4double dc1 = b71 - 21487./185130.,
89 dc2 = b72 - 0., 89 dc2 = b72 - 0.,
90 dc3 = b73 - 963225./1836901 90 dc3 = b73 - 963225./1836901.,
91 dc4 = b74 + 39864832./33354 91 dc4 = b74 + 39864832./33354255.,
92 dc5 = b75 - 2575./3519., 92 dc5 = b75 - 2575./3519.,
93 dc6 = b76 - 4472./4845., 93 dc6 = b76 - 4472./4845.,
94 dc7 = 0. + 1./10.; 94 dc7 = 0. + 1./10.;
95 95
96 // RightHandSide(yInput, dydx); 96 // RightHandSide(yInput, dydx);
97 for(G4int i = 0; i < GetNumberOfVariables( 97 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
98 { << 98 yTemp[i] = yInput[i] + hstep * b21 * dydx[i];
99 yTemp[i] = yInput[i] + hstep * b21 * dyd <<
100 } <<
101 99
102 RightHandSide(yTemp, ak2); 100 RightHandSide(yTemp, ak2);
103 for(G4int i = 0; i < GetNumberOfVariables( 101 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
104 { << 102 yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
105 yTemp[i] = yInput[i] + hstep * (b31 * dy <<
106 } <<
107 103
108 RightHandSide(yTemp, ak3); 104 RightHandSide(yTemp, ak3);
109 for(G4int i = 0;i < GetNumberOfVariables() 105 for(G4int i = 0;i < GetNumberOfVariables(); ++i)
110 { << 106 yTemp[i] = yInput[i] + hstep * (b41 * dydx[i] + b42 * ak2[i] +
111 yTemp[i] = yInput[i] + hstep * (b41 * dy << 107 b43 * ak3[i]);
112 b43 * ak <<
113 } <<
114 108
115 RightHandSide(yTemp, ak4); 109 RightHandSide(yTemp, ak4);
116 for(G4int i = 0; i < GetNumberOfVariables( 110 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
117 { << 111 yTemp[i] = yInput[i] + hstep * (b51 * dydx[i] + b52 * ak2[i] +
118 yTemp[i] = yInput[i] + hstep * (b51 * dy << 112 b53 * ak3[i] + b54 * ak4[i]);
119 b53 * ak <<
120 } <<
121 113
122 RightHandSide(yTemp, ak5); 114 RightHandSide(yTemp, ak5);
123 for(G4int i = 0; i < GetNumberOfVariables( 115 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
124 { << 116 yTemp[i] = yInput[i] + hstep * (b61 * dydx[i] + b62 * ak2[i] +
125 yTemp[i] = yInput[i] + hstep * (b61 * dy << 117 b63 * ak3[i] + b64 * ak4[i] +
126 b63 * ak << 118 b65 * ak5[i]);
127 b65 * ak <<
128 } <<
129 119
130 RightHandSide(yTemp, ak6); 120 RightHandSide(yTemp, ak6);
131 for(G4int i = 0; i < GetNumberOfVariables( 121 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
132 { << 122 yOutput[i] = yInput[i] + hstep * (b71 * dydx[i] + b72 * ak2[i] +
133 yOutput[i] = yInput[i] + hstep * (b71 * << 123 b73 * ak3[i] + b74 * ak4[i] +
134 b73 * << 124 b75 * ak5[i] + b76 * ak6[i]);
135 b75 * << 125 if (dydxOutput && yError)
136 } <<
137 if ((dydxOutput != nullptr) && (yError != <<
138 { 126 {
139 RightHandSide(yOutput, dydxOutput); 127 RightHandSide(yOutput, dydxOutput);
140 for(G4int i = 0; i < GetNumberOfVariab 128 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
141 { << 129 yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] + dc3 * ak3[i] +
142 yError[i] = hstep * (dc1 * dydx[i] + << 130 dc4 * ak4[i] + dc5 * ak5[i] + dc6 * ak6[i] +
143 dc4 * ak4[i] + << 131 dc7 * dydxOutput[i]);
144 dc7 * dydxOutpu <<
145 } <<
146 } 132 }
147 } 133 }
148 134
149 void G4RK547FEq3::Stepper( const G4double yInp 135 void G4RK547FEq3::Stepper( const G4double yInput[],
150 const G4double dydx 136 const G4double dydx[],
151 G4double hste 137 G4double hstep,
152 G4double yOut 138 G4double yOutput[],
153 G4double yErr 139 G4double yError[] )
154 { 140 {
155 copy(fyIn, yInput); 141 copy(fyIn, yInput);
156 copy(fdydx, dydx); 142 copy(fdydx, dydx);
157 fhstep = hstep; 143 fhstep = hstep;
158 144
159 makeStep(fyIn, fdydx, fhstep, fyOut, fdydx 145 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
160 146
161 copy(yOutput, fyOut); 147 copy(yOutput, fyOut);
162 } 148 }
163 149
164 void G4RK547FEq3::Stepper( const G4double yInp 150 void G4RK547FEq3::Stepper( const G4double yInput[],
165 const G4double dydx 151 const G4double dydx[],
166 G4double hste 152 G4double hstep,
167 G4double yOut 153 G4double yOutput[],
168 G4double yErr 154 G4double yError[],
169 G4double dydx 155 G4double dydxOutput[] )
170 { 156 {
171 copy(fyIn, yInput); 157 copy(fyIn, yInput);
172 copy(fdydx, dydx); 158 copy(fdydx, dydx);
173 fhstep = hstep; 159 fhstep = hstep;
174 160
175 makeStep(fyIn, fdydx, fhstep, fyOut, fdydx 161 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
176 162
177 copy(yOutput, fyOut); 163 copy(yOutput, fyOut);
178 copy(dydxOutput, fdydxOut); 164 copy(dydxOutput, fdydxOut);
179 } 165 }
180 166
181 G4double G4RK547FEq3::DistChord() const 167 G4double G4RK547FEq3::DistChord() const
182 { 168 {
183 G4double yMid[G4FieldTrack::ncompSVEC]; 169 G4double yMid[G4FieldTrack::ncompSVEC];
184 makeStep(fyIn, fdydx, fhstep / 2., yMid); 170 makeStep(fyIn, fdydx, fhstep / 2., yMid);
185 171
186 const G4ThreeVector begin = makeVector(fyI 172 const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
187 const G4ThreeVector mid = makeVector(yMid, 173 const G4ThreeVector mid = makeVector(yMid, Value3D::Position);
188 const G4ThreeVector end = makeVector(fyOut 174 const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
189 175
190 return G4LineSection::Distline(mid, begin, 176 return G4LineSection::Distline(mid, begin, end);
191 } 177 }
192 178