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