| 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 // G4RK547FEq2 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 // 2/13 | 2/13 29 // 2/13 | 2/13
32 // 2/13 | 3/52 9/52 << 30 // 2/13 | 3/52 9/52
33 // 5/9 | 12955/26244 -15925/8748 1 31 // 5/9 | 12955/26244 -15925/8748 12350/6561
34 // 3/4 | -10383/52480 13923/10496 - 32 // 3/4 | -10383/52480 13923/10496 -176553/199424 505197/997120
35 // 1 | 1403/7236 -429/268 7 33 // 1 | 1403/7236 -429/268 733330/309339 -7884/8911 104960/113967
36 // 1 | 181/2700 0 6 34 // 1 | 181/2700 0 656903/1846800 19683/106400 34112/110565 67/800
37 //-------------------------------------------- 35 //----------------------------------------------------------------------------------------------------------------------
38 // 181/2700 0 6 36 // 181/2700 0 656903/1846800 19683/106400 34112/110565 67/800 0
39 // 11377/154575 0 3 37 // 11377/154575 0 35378291/105729300 343359/1522850 535952/1947645 134/17175 1/12
40 // <<
41 // Author: Dmitry Sorokin, Google Summer of Co <<
42 // Supervision: John Apostolakis, CERN <<
43 // ------------------------------------------- <<
44 38
45 #include "G4RK547FEq2.hh" 39 #include "G4RK547FEq2.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 namespace {
>> 46
>> 47 void copyArray(G4double dst[], const G4double src[])
>> 48 {
>> 49 memcpy(dst, src, sizeof(G4double) * G4FieldTrack::ncompSVEC);
>> 50 }
>> 51
>> 52 } // namespace
>> 53
51 G4RK547FEq2::G4RK547FEq2(G4EquationOfMotion* E 54 G4RK547FEq2::G4RK547FEq2(G4EquationOfMotion* EqRhs, G4int integrationVariables)
52 : G4MagIntegratorStepper(EqRhs, integrationV << 55 : G4MagIntegratorStepper(EqRhs, integrationVariables)
53 { 56 {
54 } 57 }
55 58
56 void G4RK547FEq2::makeStep( const G4double yIn << 59 void G4RK547FEq2::makeStep(
57 const G4double dyd << 60 const G4double yInput[],
58 const G4double hst << 61 const G4double dydx[],
59 G4double yOu << 62 const G4double hstep,
60 G4double* dy << 63 G4double yOutput[],
61 G4double* yE << 64 G4double* dydxOutput,
>> 65 G4double* yError) const
62 { 66 {
63 G4double yTemp[G4FieldTrack::ncompSVEC]; 67 G4double yTemp[G4FieldTrack::ncompSVEC];
64 for (G4int i=GetNumberOfVariables(); i<Get << 68 for (int i = GetNumberOfVariables(); i < GetNumberOfStateVariables(); ++i){
65 { <<
66 yOutput[i] = yTemp[i] = yInput[i]; 69 yOutput[i] = yTemp[i] = yInput[i];
67 } 70 }
68 71
69 G4double ak2[G4FieldTrack::ncompSVEC], 72 G4double ak2[G4FieldTrack::ncompSVEC],
70 ak3[G4FieldTrack::ncompSVEC], 73 ak3[G4FieldTrack::ncompSVEC],
71 ak4[G4FieldTrack::ncompSVEC], 74 ak4[G4FieldTrack::ncompSVEC],
72 ak5[G4FieldTrack::ncompSVEC], 75 ak5[G4FieldTrack::ncompSVEC],
73 ak6[G4FieldTrack::ncompSVEC]; 76 ak6[G4FieldTrack::ncompSVEC];
74 77
75 const G4double b21 = 2./13., << 78 const G4double
76 b31 = 3./52., b32 = 9./52., << 79 b21 = 2./13.,
77 b41 = 12955./26244., b42 = << 80 b31 = 3./52., b32 = 9./52.,
78 b43 = 12350./6561., << 81 b41 = 12955./26244., b42 = -15925./8748., b43 = 12350./6561.,
79 b51 = -10383./52480., b52 = << 82 b51 = -10383./52480., b52 = 13923./10496., b53 = -176553./199424.,
80 b53 = -176553./199424., << 83 b54 = 505197./997120.,
81 b61 = 1403./7236., b62 = -4 << 84 b61 = 1403./7236., b62 = -429./268., b63 = 733330./309339.,
82 b64 = -7884./8911., b65 << 85 b64 = -7884./8911., b65 = 104960./113967.,
83 b71 = 181./2700., b72 = 0., << 86 b71 = 181./2700., b72 = 0., b73 = 656903./1846800.,
84 b74 = 19683./106400., b << 87 b74 = 19683./106400., b75 = 34112./110565., b76 = 67./800.;
85 b76 = 67./800.; << 88
86 << 89 const G4double
87 const G4double dc1 = b71 - 11377./154575., << 90 dc1 = b71 - 11377./154575.,
88 dc2 = b72 - 0., << 91 dc2 = b72 - 0.,
89 dc3 = b73 - 35378291./10572 << 92 dc3 = b73 - 35378291./105729300.,
90 dc4 = b74 - 343359./1522850 << 93 dc4 = b74 - 343359./1522850.,
91 dc5 = b75 - 535952./1947645 << 94 dc5 = b75 - 535952./1947645.,
92 dc6 = b76 - 134./17175., << 95 dc6 = b76 - 134./17175.,
93 dc7 = 0. - 1./12.; << 96 dc7 = 0. - 1./12.;
94 << 97
95 // RightHandSide(yInput, dydx); << 98 //RightHandSide(yInput, dydx);
96 for(G4int i = 0; i < GetNumberOfVariables( << 99 for(int i = 0; i < GetNumberOfVariables(); ++i)
97 { << 100 yTemp[i] = yInput[i] + hstep * b21 * dydx[i];
98 yTemp[i] = yInput[i] + hstep * b21 * dyd <<
99 } <<
100 101
101 RightHandSide(yTemp, ak2); 102 RightHandSide(yTemp, ak2);
102 for(G4int i = 0; i < GetNumberOfVariables( << 103 for(int i = 0; i < GetNumberOfVariables(); ++i)
103 { << 104 yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
104 yTemp[i] = yInput[i] + hstep * (b31 * dy <<
105 } <<
106 105
107 RightHandSide(yTemp, ak3); 106 RightHandSide(yTemp, ak3);
108 for(G4int i = 0;i < GetNumberOfVariables() << 107 for(int i = 0;i < GetNumberOfVariables(); ++i)
109 { << 108 yTemp[i] = yInput[i] + hstep * (b41 * dydx[i] + b42 * ak2[i] +
110 yTemp[i] = yInput[i] + hstep * (b41 * dy << 109 b43 * ak3[i]);
111 b43 * ak <<
112 } <<
113 110
114 RightHandSide(yTemp, ak4); 111 RightHandSide(yTemp, ak4);
115 for(G4int i = 0; i < GetNumberOfVariables( << 112 for(int i = 0; i < GetNumberOfVariables(); ++i)
116 { << 113 yTemp[i] = yInput[i] + hstep * (b51 * dydx[i] + b52 * ak2[i] +
117 yTemp[i] = yInput[i] + hstep * (b51 * dy << 114 b53 * ak3[i] + b54 * ak4[i]);
118 b53 * ak <<
119 } <<
120 115
121 RightHandSide(yTemp, ak5); 116 RightHandSide(yTemp, ak5);
122 for(G4int i = 0; i < GetNumberOfVariables( << 117 for(int i = 0; i < GetNumberOfVariables(); ++i)
123 { << 118 yTemp[i] = yInput[i] + hstep * (b61 * dydx[i] + b62 * ak2[i] +
124 yTemp[i] = yInput[i] + hstep * (b61 * dy << 119 b63 * ak3[i] + b64 * ak4[i] +
125 b63 * ak << 120 b65 * ak5[i]);
126 b65 * ak <<
127 } <<
128 121
129 RightHandSide(yTemp, ak6); 122 RightHandSide(yTemp, ak6);
130 for(G4int i = 0; i < GetNumberOfVariables( << 123 for(int i = 0; i < GetNumberOfVariables(); ++i)
131 { << 124 yOutput[i] = yInput[i] + hstep * (b71 * dydx[i] + b72 * ak2[i] +
132 yOutput[i] = yInput[i] + hstep * (b71 * << 125 b73 * ak3[i] + b74 * ak4[i] +
133 b73 * << 126 b75 * ak5[i] + b76 * ak6[i]);
134 b75 * << 127
135 } << 128 if (dydxOutput && yError) {
136 if ((dydxOutput != nullptr) && (yError != <<
137 { <<
138 RightHandSide(yOutput, dydxOutput); 129 RightHandSide(yOutput, dydxOutput);
139 for(G4int i = 0; i < GetNumberOfVariab << 130 for(int i = 0; i < GetNumberOfVariables(); ++i)
140 { << 131 yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] + dc3 * ak3[i] +
141 yError[i] = hstep * (dc1 * dydx[i] + << 132 dc4 * ak4[i] + dc5 * ak5[i] + dc6 * ak6[i] +
142 dc4 * ak4[i] + << 133 dc7 * dydxOutput[i]);
143 dc7 * dydxOutpu <<
144 } <<
145 } 134 }
146 } 135 }
147 136
148 void G4RK547FEq2::Stepper( const G4double yInp << 137 void G4RK547FEq2::Stepper(
149 const G4double dydx << 138 const G4double yInput[],
150 G4double hste << 139 const G4double dydx[],
151 G4double yOut << 140 G4double hstep,
152 G4double yErr << 141 G4double yOutput[],
>> 142 G4double yError[])
153 { 143 {
154 copy(fyIn, yInput); << 144 copyArray(fyIn, yInput);
155 copy(fdydx, dydx); << 145 copyArray(fdydx, dydx);
156 fhstep = hstep; 146 fhstep = hstep;
157 147
158 makeStep(fyIn, fdydx, fhstep, fyOut, fdydx 148 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
159 149
160 copy(yOutput, fyOut); << 150 copyArray(yOutput, fyOut);
161 } 151 }
162 152
163 void G4RK547FEq2::Stepper( const G4double yInp << 153 void G4RK547FEq2::Stepper(
164 const G4double dydx << 154 const G4double yInput[],
165 G4double hste << 155 const G4double dydx[],
166 G4double yOut << 156 G4double hstep,
167 G4double yErr << 157 G4double yOutput[],
168 G4double dydx << 158 G4double yError[],
>> 159 G4double dydxOutput[])
169 { 160 {
170 copy(fyIn, yInput); << 161 copyArray(fyIn, yInput);
171 copy(fdydx, dydx); << 162 copyArray(fdydx, dydx);
172 fhstep = hstep; 163 fhstep = hstep;
173 164
174 makeStep(fyIn, fdydx, fhstep, fyOut, fdydx 165 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
175 166
176 copy(yOutput, fyOut); << 167 copyArray(yOutput, fyOut);
177 copy(dydxOutput, fdydxOut); << 168 copyArray(dydxOutput, fdydxOut);
178 } 169 }
179 170
180 G4double G4RK547FEq2::DistChord() const 171 G4double G4RK547FEq2::DistChord() const
181 { 172 {
182 G4double yMid[G4FieldTrack::ncompSVEC]; 173 G4double yMid[G4FieldTrack::ncompSVEC];
183 makeStep(fyIn, fdydx, fhstep / 2., yMid); 174 makeStep(fyIn, fdydx, fhstep / 2., yMid);
184 175
185 const G4ThreeVector begin = makeVector(fyI 176 const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
186 const G4ThreeVector mid = makeVector(yMid, 177 const G4ThreeVector mid = makeVector(yMid, Value3D::Position);
187 const G4ThreeVector end = makeVector(fyOut 178 const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
188 179
189 return G4LineSection::Distline(mid, begin, 180 return G4LineSection::Distline(mid, begin, end);
190 } 181 }
191 182