| 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 // G4RKG3_Stepper implementation <<
27 // 26 //
28 // Created: J.Apostolakis, V.Grichine - 30.01. << 27 // $Id: G4RKG3_Stepper.cc,v 1.12 2007/04/26 12:23:55 tnikitin Exp $
>> 28 // GEANT4 tag $Name: geant4-08-03-patch-01 $
>> 29 //
29 // ------------------------------------------- 30 // -------------------------------------------------------------------
30 31
31 #include "G4RKG3_Stepper.hh" 32 #include "G4RKG3_Stepper.hh"
32 #include "G4LineSection.hh" 33 #include "G4LineSection.hh"
33 #include "G4Mag_EqRhs.hh" 34 #include "G4Mag_EqRhs.hh"
34 35
35 G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs* Eq << 36 G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs *EqRhs)
36 : G4MagIntegratorStepper(EqRhs,6) 37 : G4MagIntegratorStepper(EqRhs,6)
37 { 38 {
>> 39 // G4Exception("G4RKG3_Stepper::G4RKG3_Stepper()", "NotImplemented",
>> 40 // FatalException, "Stepper not yet available.");
38 } 41 }
39 42
40 G4RKG3_Stepper::~G4RKG3_Stepper() = default; << 43 G4RKG3_Stepper::~G4RKG3_Stepper()
>> 44 {
>> 45 }
41 46
42 void G4RKG3_Stepper::Stepper( const G4double y << 47 void G4RKG3_Stepper::Stepper( const G4double yInput[7],
43 const G4double d << 48 const G4double dydx[7],
44 G4double S << 49 G4double Step,
45 G4double y << 50 G4double yOut[7],
46 G4double y << 51 G4double yErr[])
47 { 52 {
48 G4double B[3]; 53 G4double B[3];
>> 54 // G4double yderiv[6];
>> 55 // G4double alpha2, beta2;
49 G4int nvar = 6 ; 56 G4int nvar = 6 ;
50 G4double by15 = 1. / 15. ; // was 0.06666 << 57 // G4double beTemp2, beta2=0;
51 58
52 G4double yTemp[8], dydxTemp[6], yIn[8]; << 59 G4int i;
>> 60 G4double by15 = 1. / 15. ; // was 0.066666666 ;
>> 61 G4double yTemp[7], dydxTemp[6], yIn[7] ;
>> 62 // Saving yInput because yInput and yOut can be aliases for same array
>> 63 for(i=0;i<nvar;i++) yIn[i]=yInput[i];
53 64
54 // Saving yInput because yInput and yOut ca <<
55 // <<
56 for(G4int i=0; i<nvar; ++i) <<
57 { <<
58 yIn[i]=yInput[i]; <<
59 } <<
60 yIn[6] = yInput[6]; <<
61 yIn[7] = yInput[7]; <<
62 G4double h = Step * 0.5; 65 G4double h = Step * 0.5;
63 hStep = Step; <<
64 // Do two half steps <<
65 66
66 StepNoErr(yIn, dydx,h, yTemp,B) ; << 67 // Do two half steps
67 <<
68 // Store Bfld for DistChord Calculation <<
69 // <<
70 for(auto i=0; i<3; ++i) <<
71 { <<
72 BfldIn[i] = B[i]; <<
73 } <<
74 // RightHandSide(yTemp,dydxTemp) ; <<
75 68
>> 69
>> 70 // To obtain B1 ...
>> 71 // GetEquationOfMotion()->GetFieldValue(yIn,B);
>> 72 // G4RKG3_Stepper::StepWithEst(yIn, dydx, Step, yOut,alpha2, beta2, B1, B2 );
>> 73
>> 74 StepNoErr(yIn, dydx,h, yTemp,B) ;
>> 75 // RightHandSide(yTemp,dydxTemp) ;
76 GetEquationOfMotion()->EvaluateRhsGivenB(yT 76 GetEquationOfMotion()->EvaluateRhsGivenB(yTemp,B,dydxTemp) ;
77 StepNoErr(yTemp,dydxTemp,h,yOut,B); << 77 StepNoErr(yTemp,dydxTemp,h,yOut,B); // ,beTemp2) ;
78 << 78 // beta2 += beTemp2;
>> 79 // beta2 *= 0.5;
>> 80
79 // Store midpoint, chord calculation 81 // Store midpoint, chord calculation
80 82
81 fyMidPoint = G4ThreeVector(yTemp[0], yTemp << 83 fyMidPoint = G4ThreeVector( yTemp[0], yTemp[1], yTemp[2]);
82 84
83 // Do a full Step 85 // Do a full Step
84 // << 86
85 h *= 2 ; 87 h *= 2 ;
86 StepNoErr(yIn,dydx,h,yTemp,B); << 88 StepNoErr(yIn,dydx,h,yTemp,B); // ,beTemp2) ;
87 for(G4int i=0; i<nvar; ++i) << 89 for(i=0;i<nvar;i++)
88 { 90 {
89 yErr[i] = yOut[i] - yTemp[i] ; 91 yErr[i] = yOut[i] - yTemp[i] ;
90 yOut[i] += yErr[i]*by15 ; // Pr 92 yOut[i] += yErr[i]*by15 ; // Provides 5th order of accuracy
91 } 93 }
92 94
93 // Store values for DistChord method << 95 // for(i=0;i<ncomp;i++)
94 // << 96 // {
95 fyInitial = G4ThreeVector( yIn[0], yIn[1] << 97 // fyInitial[i] = yIn[i];
96 fpInitial = G4ThreeVector( yIn[3], yIn[4] << 98 // fyFinal[i] = yOut[i];
>> 99 // }
>> 100
>> 101 fyInitial = G4ThreeVector( yIn[0], yIn[1], yIn[2]);
97 fyFinal = G4ThreeVector( yOut[0], yOut[1 102 fyFinal = G4ThreeVector( yOut[0], yOut[1], yOut[2]);
>> 103 // beta2 += beTemp2 ;
>> 104 // beta2 *= 0.5 ;
>> 105 // NormaliseTangentVector( yOut ); // Deleted
98 } 106 }
99 107
100 // ------------------------------------------- 108 // ---------------------------------------------------------------------------
101 109
102 // Integrator for RK from G3 with evaluation o 110 // Integrator for RK from G3 with evaluation of error in solution and delta
103 // geometry based on naive similarity with the 111 // geometry based on naive similarity with the case of uniform magnetic field.
104 // B1[3] is input and is the first magnetic f 112 // B1[3] is input and is the first magnetic field values
105 // B2[3] is output and is the final magnetic f 113 // B2[3] is output and is the final magnetic field values.
106 // << 114
107 void G4RKG3_Stepper::StepWithEst( const G4doub 115 void G4RKG3_Stepper::StepWithEst( const G4double*,
108 const G4doub 116 const G4double*,
109 G4doub 117 G4double,
110 G4doub 118 G4double*,
111 G4doub 119 G4double&,
112 G4doub 120 G4double&,
113 const G4doub 121 const G4double*,
114 G4doub 122 G4double* )
115 123
116 { 124 {
117 G4Exception("G4RKG3_Stepper::StepWithEst()", << 125 G4Exception("G4RKG3_Stepper::StepWithEst()", "ObsoleteMethod",
118 FatalException, "Method no longe 126 FatalException, "Method no longer used.");
119 } 127 }
120 128
121 // ------------------------------------------- 129 // -----------------------------------------------------------------
122 130
>> 131
123 // Integrator RK Stepper from G3 with only two 132 // Integrator RK Stepper from G3 with only two field evaluation per Step.
124 // It is used in propagation initial Step by s 133 // It is used in propagation initial Step by small substeps after solution
125 // error and delta geometry considerations. B[ 134 // error and delta geometry considerations. B[3] is magnetic field which
126 // is passed from substep to substep. 135 // is passed from substep to substep.
127 // << 136
128 void G4RKG3_Stepper::StepNoErr(const G4double << 137 void G4RKG3_Stepper::StepNoErr(const G4double tIn[7],
129 const G4double << 138 const G4double dydx[7],
130 G4double 139 G4double Step,
131 G4double << 140 G4double tOut[7],
132 G4double << 141 G4double B[3] ) // const
133 142
134 { << 143 {
135 144
136 // Copy and edit the routine above, to dele << 145 // Copy and edit the routine above, to delete alpha2, beta2, ...
137 // << 146 G4double K1[7],K2[7],K3[7],K4[7] ;
138 G4double K1[7], K2[7], K3[7], K4[7]; << 147 G4double tTemp[7], yderiv[6] ;
139 G4double tTemp[8]={0.0}, yderiv[6]={0.0}; << 148 // Need Momentum value to give correct values to the coefficients in equation
140 << 149 // Integration on unit velocity, but tIn[3,4,5] is momentum
141 // Need Momentum value to give correct valu << 150 G4double mom;
142 // equation. Integration on unit velocity, << 151 G4int i ;
143 <<
144 G4double mom, inverse_mom; <<
145 const G4double c1=0.5, c2=0.125, c3=1./6.; <<
146 152
>> 153 #ifdef END_CODE_G3STEPPER
>> 154 G4Exception(" G4RKG3_Stepper::StepNoErr(): method to be no longer used.");
>> 155 #else
>> 156
>> 157 // GetEquationOfMotion()->EvaluateRhsReturnB(tIn,dydx,B1) ;
147 // Correction for momentum not a velocity 158 // Correction for momentum not a velocity
148 // Need the protection !!! must be not zero << 159 mom=std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4]+tIn[5]*tIn[5]);
149 // << 160
150 mom = std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4] << 161 for(i=0;i<3;i++)
151 inverse_mom = 1./mom; << 162 {
152 for(auto i=0; i<3; ++i) << 163 K1[i] = Step * dydx[i+3]/mom;
153 { << 164 tTemp[i] = tIn[i] + Step*(0.5*tIn[i+3]/mom + 0.125*K1[i]) ;
154 K1[i] = Step * dydx[i+3]*inverse_mom; << 165 tTemp[i+3] = tIn[i+3] + 0.5*K1[i]*mom ;
155 tTemp[i] = tIn[i] + Step*(c1*tIn[i+3]*in << 166
156 tTemp[i+3] = tIn[i+3] + c1*K1[i]*mom ; <<
157 } 167 }
158 168
159 GetEquationOfMotion()->EvaluateRhsReturnB(t 169 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
>> 170
160 171
161 for(auto i=0; i<3; ++i) << 172 for(i=0;i<3;i++)
162 { 173 {
163 K2[i] = Step * yderiv[i+3]*inverse_mom; << 174 K2[i] = Step * yderiv[i+3]/mom;
164 tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ; << 175 tTemp[i+3] = tIn[i+3] + 0.5*K2[i]*mom ;
165 } 176 }
166 177
167 // Given B, calculate yderiv ! << 178 // Given B, calculate yderiv !
168 // <<
169 GetEquationOfMotion()->EvaluateRhsGivenB(tT 179 GetEquationOfMotion()->EvaluateRhsGivenB(tTemp,B,yderiv) ;
170 180
171 for(auto i=0; i<3; ++i) << 181 for(i=0;i<3;i++)
172 { 182 {
173 K3[i] = Step * yderiv[i+3]*inverse_mom; << 183 K3[i] = Step * yderiv[i+3]/mom;
174 tTemp[i] = tIn[i] + Step*(tIn[i+3]*inver << 184 tTemp[i] = tIn[i] + Step*(tIn[i+3]/mom + 0.5*K3[i]) ;
175 tTemp[i+3] = tIn[i+3] + K3[i]*mom ; 185 tTemp[i+3] = tIn[i+3] + K3[i]*mom ;
176 } 186 }
>> 187
177 188
178 // Calculates y-deriv(atives) & returns B t << 189 // Calculates y-deriv(atives) & returns B too!
179 // <<
180 GetEquationOfMotion()->EvaluateRhsReturnB(t 190 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
>> 191
181 192
182 for(auto i=0; i<3; ++i) // Output tr << 193 for(i=0;i<3;i++) // Output trajectory vector
183 { 194 {
184 K4[i] = Step * yderiv[i+3]*inverse_mom; << 195 K4[i] = Step * yderiv[i+3]/mom;
185 tOut[i] = tIn[i] + Step*(tIn[i+3]*invers << 196 tOut[i] = tIn[i] + Step*(tIn[i+3]/mom + (K1[i] + K2[i] + K3[i])/6.0) ;
186 tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2 << 197 tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2[i] + 2*K3[i] +K4[i])/6.0 ;
187 } 198 }
188 tOut[6] = tIn[6]; << 199
189 tOut[7] = tIn[7]; << 200 // NormaliseTangentVector( tOut );
>> 201
>> 202 #endif
190 } 203 }
191 <<
192 // ------------------------------------------- 204 // ---------------------------------------------------------------------------
193 <<
194 G4double G4RKG3_Stepper::DistChord() const <<
195 { <<
196 // Soon: must check whether h/R > 2 pi !! <<
197 // Method below is good only for < 2 pi <<
198 205
199 G4double distChord,distLine; << 206 G4double G4RKG3_Stepper::DistChord() const
200 << 207 {
201 if (fyInitial != fyFinal) << 208 // Soon: must check whether h/R > 2 pi !!
202 { << 209 // Method below is good only for < 2 pi
203 distLine = G4LineSection::Distline(fyMid << 210 G4double distChord,distLine;
204 distChord = distLine; << 211
205 } << 212 if (fyInitial != fyFinal) {
206 else << 213 distLine= G4LineSection::Distline(fyMidPoint,fyInitial,fyFinal );
207 { << 214
208 distChord = (fyMidPoint-fyInitial).mag() << 215 distChord = distLine;
209 } << 216 }else{
>> 217 distChord = (fyMidPoint-fyInitial).mag();
>> 218 }
210 219
211 return distChord; << 220
>> 221 return distChord;
>> 222
212 } 223 }
213 224