| 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 // G4DoLoMcPriRK34 implementation
27 //
28 // Created: Somnath Banerjee, Google Summer of
29 // Supervision: John Apostolakis, CERN
30 // -------------------------------------------
31
32 #include "G4DoLoMcPriRK34.hh"
33 #include "G4LineSection.hh"
34
35 // Constructor
36 //
37 G4DoLoMcPriRK34::G4DoLoMcPriRK34(G4EquationOfM
38 G4int noInte
39 G4bool prima
40 : G4MagIntegratorStepper(EqRhs, noIntegrati
41 {
42 const G4int numberOfVariables = noIntegrati
43
44 // New Chunk of memory being created for us
45
46 // aki - for storing intermediate RHS
47 //
48 ak2 = new G4double[numberOfVariables];
49 ak3 = new G4double[numberOfVariables];
50 ak4 = new G4double[numberOfVariables];
51 ak5 = new G4double[numberOfVariables];
52 ak6 = new G4double[numberOfVariables];
53
54 yTemp = new G4double[numberOfVariables] ;
55 yIn = new G4double[numberOfVariables] ;
56
57 fLastInitialVector = new G4double[numberOfV
58 fLastFinalVector = new G4double[numberOfVar
59 fLastDyDx = new G4double[numberOfVariables]
60
61 fMidVector = new G4double[numberOfVariables
62 fMidError = new G4double[numberOfVariables]
63 if( primary )
64 {
65 fAuxStepper = new G4DoLoMcPriRK34(EqRhs,
66 }
67 }
68
69 // Destructor
70 //
71 G4DoLoMcPriRK34::~G4DoLoMcPriRK34()
72 {
73 // clear all previously allocated memory fo
74
75 delete [] ak2;
76 delete [] ak3;
77 delete [] ak4;
78 delete [] ak5;
79 delete [] ak6;
80
81 delete [] yTemp;
82 delete [] yIn;
83
84 delete [] fLastInitialVector;
85 delete [] fLastFinalVector;
86 delete [] fLastDyDx;
87 delete [] fMidVector;
88 delete [] fMidError;
89
90 delete fAuxStepper;
91 }
92
93 // Stepper
94 //
95 // Passing in the value of yInput[],the first
96 // Giving back yOut and yErr arrays for output
97 //
98 void G4DoLoMcPriRK34::Stepper(const G4double y
99 const G4double D
100 G4double S
101 G4double y
102 G4double y
103 {
104 G4int i;
105
106 // The various constants defined on the bas
107 //
108 const G4double b21 = 7.0/27.0 ,
109 b31 = 7.0/72.0 ,
110 b32 = 7.0/24.0 ,
111
112 b41 = 3043.0/3528.0 ,
113 b42 = -3757.0/1176.0 ,
114 b43 = 1445.0/441.0,
115
116 b51 = 17617.0/11662.0 ,
117 b52 = -4023.0/686.0 ,
118 b53 = 9372.0/1715.0 ,
119 b54 = -66.0/595.0 ,
120
121 b61 = 29.0/238.0 ,
122 b62 = 0.0 ,
123 b63 = 216.0/385.0 ,
124 b64 = 54.0/85.0 ,
125 b65 = -7.0/22.0 ,
126
127 dc1 = 363.0/2975.0 - b61 ,
128 dc2 = 0.0 - b62 ,
129 dc3 = 981.0/1750.0 - b63,
130 dc4 = 2709.0/4250.0 - b64 ,
131 dc5 = -3.0/10.0 - b65 ,
132 dc6 = -1.0/50.0 ; //
133
134 const G4int numberOfVariables = GetNumberOf
135
136 // The number of variables to be integrated
137 //
138 yOut[7] = yTemp[7] = yIn[7];
139
140 // Saving yInput because yInput and yOut ca
141 //
142 for(i=0; i<numberOfVariables; ++i)
143 {
144 yIn[i]=yInput[i];
145 }
146
147 // RightHandSide(yIn, DyDx) ; // 1st stag
148
149 for(i=0; i<numberOfVariables; ++i)
150 {
151 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
152 }
153 RightHandSide(yTemp, ak2) ; //
154
155 for(i=0; i<numberOfVariables; ++i)
156 {
157 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] +
158 }
159 RightHandSide(yTemp, ak3) ; //
160
161 for(i=0; i<numberOfVariables; ++i)
162 {
163 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] +
164 }
165 RightHandSide(yTemp, ak4) ; //
166
167 for(i=0; i<numberOfVariables; ++i)
168 {
169 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] +
170 + b53*ak3[i] +
171 }
172 RightHandSide(yTemp, ak5) ; //
173
174 for(i=0; i<numberOfVariables; ++i)
175 {
176 yOut[i] = yIn[i] + Step*(b61*DyDx[i] +
177 + b64*ak4[i] + b
178 }
179 RightHandSide(yOut, ak6) ; //
180
181 for(i=0; i<numberOfVariables; ++i)
182 {
183
184 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i
185 + dc5*ak5[i] + dc6*ak6[i]
186
187 // Store Input and Final values, for po
188 //
189 fLastInitialVector[i] = yIn[i] ;
190 fLastFinalVector[i] = yOut[i];
191 fLastDyDx[i] = DyDx[i];
192 }
193
194 fLastStepLength = Step;
195
196 return ;
197 }
198
199 // DistChord
200 //
201 G4double G4DoLoMcPriRK34::DistChord() const
202 {
203 G4double distLine, distChord;
204 G4ThreeVector initialPoint, finalPoint, mid
205
206 // Store last initial and final points
207 // (they will be overwritten in self-Steppe
208 //
209 initialPoint = G4ThreeVector( fLastInitialV
210 fLastInitialV
211 finalPoint = G4ThreeVector( fLastFinalVec
212 fLastFinalVec
213
214 // Do half a Step using StepNoErr
215
216 fAuxStepper->Stepper( fLastInitialVector, f
217 fMidVector, fMidError
218
219 midPoint = G4ThreeVector( fMidVector[0], fM
220
221 // Use stored values of Initial and Endpoin
222 // distance of Chord
223 //
224 if (initialPoint != finalPoint)
225 {
226 distLine = G4LineSection::Distline( midP
227 distChord = distLine;
228 }
229 else
230 {
231 distChord = (midPoint-initialPoint).mag()
232 }
233 return distChord;
234 }
235
236 void G4DoLoMcPriRK34::SetupInterpolation()
237 {
238 }
239
240 void G4DoLoMcPriRK34::SetupInterpolate( const
241 const
242 const
243 {
244 // Do Nothing
245 }
246
247 void G4DoLoMcPriRK34::Interpolate( G4double ta
248 G4double yO
249 {
250 Interpolate( fLastInitialVector, fLastDyDx,
251 }
252
253 // Function to evaluate the interpolation at t
254 //
255 void G4DoLoMcPriRK34::Interpolate( const G4dou
256 const G4dou
257 const G4dou
258 G4dou
259 G4dou
260 {
261 G4double bf1, bf2, bf3, bf4, bf5, bf6;
262
263 const G4int numberOfVariables = GetNumberOf
264
265 for(G4int i=0; i<numberOfVariables; ++i)
266 {
267 yIn[i]=yInput[i];
268 }
269
270 G4double tau_2 = tau*tau, tau_3 = tau*tau_2
271
272 // Calculating the polynomials (coefficient
273 //
274 bf1 = -(162.0*tau_3 - 504.0*tau_2 + 551.0*t
275 bf2 = 0.0 ,
276 bf3 = 27.0*tau*(27.0*tau_2 - 70.0*tau + 51
277 bf4 = -27*tau*(27.0*tau_2 - 50.0*tau + 21.0
278 bf5 = 7.0*tau*(2232.0*tau_2 - 4166.0*tau +
279 bf6 = tau*(tau - 1.0)*(387.0*tau - 238.0)/1
280
281 for( G4int i=0; i<numberOfVariables; ++i)
282 {
283 yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i]
284 + bf4*ak4[i] + bf5*ak5[i] + bf6*a
285 }
286 }
287