Geant4-11
G4RKG3_Stepper.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4RKG3_Stepper implementation
27//
28// Created: J.Apostolakis, V.Grichine - 30.01.1997
29// -------------------------------------------------------------------
30
31#include "G4RKG3_Stepper.hh"
32#include "G4LineSection.hh"
33#include "G4Mag_EqRhs.hh"
34
36 : G4MagIntegratorStepper(EqRhs,6)
37{
38}
39
41{
42}
43
44void G4RKG3_Stepper::Stepper( const G4double yInput[8],
45 const G4double dydx[6],
46 G4double Step,
47 G4double yOut[8],
48 G4double yErr[] )
49{
50 G4double B[3];
51 G4int nvar = 6 ;
52 G4double by15 = 1. / 15. ; // was 0.066666666 ;
53
54 G4double yTemp[8], dydxTemp[6], yIn[8];
55
56 // Saving yInput because yInput and yOut can be aliases for same array
57 //
58 for(G4int i=0; i<nvar; ++i)
59 {
60 yIn[i]=yInput[i];
61 }
62 yIn[6] = yInput[6];
63 yIn[7] = yInput[7];
64 G4double h = Step * 0.5;
65 hStep = Step;
66 // Do two half steps
67
68 StepNoErr(yIn, dydx,h, yTemp,B) ;
69
70 // Store Bfld for DistChord Calculation
71 //
72 for(auto i=0; i<3; ++i)
73 {
74 BfldIn[i] = B[i];
75 }
76 // RightHandSide(yTemp,dydxTemp) ;
77
78 GetEquationOfMotion()->EvaluateRhsGivenB(yTemp,B,dydxTemp) ;
79 StepNoErr(yTemp,dydxTemp,h,yOut,B);
80
81 // Store midpoint, chord calculation
82
83 fyMidPoint = G4ThreeVector(yTemp[0], yTemp[1], yTemp[2]);
84
85 // Do a full Step
86 //
87 h *= 2 ;
88 StepNoErr(yIn,dydx,h,yTemp,B);
89 for(G4int i=0; i<nvar; ++i)
90 {
91 yErr[i] = yOut[i] - yTemp[i] ;
92 yOut[i] += yErr[i]*by15 ; // Provides 5th order of accuracy
93 }
94
95 // Store values for DistChord method
96 //
97 fyInitial = G4ThreeVector( yIn[0], yIn[1], yIn[2]);
98 fpInitial = G4ThreeVector( yIn[3], yIn[4], yIn[5]);
99 fyFinal = G4ThreeVector( yOut[0], yOut[1], yOut[2]);
100}
101
102// ---------------------------------------------------------------------------
103
104// Integrator for RK from G3 with evaluation of error in solution and delta
105// geometry based on naive similarity with the case of uniform magnetic field.
106// B1[3] is input and is the first magnetic field values
107// B2[3] is output and is the final magnetic field values.
108//
110 const G4double*,
111 G4double,
112 G4double*,
113 G4double&,
114 G4double&,
115 const G4double*,
116 G4double* )
117
118{
119 G4Exception("G4RKG3_Stepper::StepWithEst()", "GeomField0001",
120 FatalException, "Method no longer used.");
121}
122
123// -----------------------------------------------------------------
124
125// Integrator RK Stepper from G3 with only two field evaluation per Step.
126// It is used in propagation initial Step by small substeps after solution
127// error and delta geometry considerations. B[3] is magnetic field which
128// is passed from substep to substep.
129//
131 const G4double dydx[6],
132 G4double Step,
133 G4double tOut[8],
134 G4double B[3] )
135
136{
137
138 // Copy and edit the routine above, to delete alpha2, beta2, ...
139 //
140 G4double K1[7], K2[7], K3[7], K4[7];
141 G4double tTemp[8]={0.0}, yderiv[6]={0.0};
142
143 // Need Momentum value to give correct values to the coefficients in
144 // equation. Integration on unit velocity, but tIn[3,4,5] is momentum
145
146 G4double mom, inverse_mom;
147 const G4double c1=0.5, c2=0.125, c3=1./6.;
148
149 // Correction for momentum not a velocity
150 // Need the protection !!! must be not zero
151 //
152 mom = std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4]+tIn[5]*tIn[5]);
153 inverse_mom = 1./mom;
154 for(auto i=0; i<3; ++i)
155 {
156 K1[i] = Step * dydx[i+3]*inverse_mom;
157 tTemp[i] = tIn[i] + Step*(c1*tIn[i+3]*inverse_mom + c2*K1[i]) ;
158 tTemp[i+3] = tIn[i+3] + c1*K1[i]*mom ;
159 }
160
161 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
162
163 for(auto i=0; i<3; ++i)
164 {
165 K2[i] = Step * yderiv[i+3]*inverse_mom;
166 tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ;
167 }
168
169 // Given B, calculate yderiv !
170 //
171 GetEquationOfMotion()->EvaluateRhsGivenB(tTemp,B,yderiv) ;
172
173 for(auto i=0; i<3; ++i)
174 {
175 K3[i] = Step * yderiv[i+3]*inverse_mom;
176 tTemp[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom + c1*K3[i]) ;
177 tTemp[i+3] = tIn[i+3] + K3[i]*mom ;
178 }
179
180 // Calculates y-deriv(atives) & returns B too!
181 //
182 GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
183
184 for(auto i=0; i<3; ++i) // Output trajectory vector
185 {
186 K4[i] = Step * yderiv[i+3]*inverse_mom;
187 tOut[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom+ (K1[i]+K2[i]+K3[i])*c3) ;
188 tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2[i] + 2*K3[i] +K4[i])*c3 ;
189 }
190 tOut[6] = tIn[6];
191 tOut[7] = tIn[7];
192}
193
194// ---------------------------------------------------------------------------
195
197{
198 // Soon: must check whether h/R > 2 pi !!
199 // Method below is good only for < 2 pi
200
201 G4double distChord,distLine;
202
203 if (fyInitial != fyFinal)
204 {
206 distChord = distLine;
207 }
208 else
209 {
210 distChord = (fyMidPoint-fyInitial).mag();
211 }
212
213 return distChord;
214}
G4double B(G4double temperature)
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
CLHEP::Hep3Vector G4ThreeVector
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
virtual void EvaluateRhsGivenB(const G4double y[], const G4double B[3], G4double dydx[]) const =0
void EvaluateRhsReturnB(const G4double y[], G4double dydx[], G4double Field[]) const
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
G4EquationOfMotion * GetEquationOfMotion()
G4ThreeVector fyInitial
void StepWithEst(const G4double tIn[8], const G4double dydx[6], G4double Step, G4double tOut[8], G4double &alpha2, G4double &beta2, const G4double B1[3], G4double B2[3])
void Stepper(const G4double yIn[], const G4double dydx[], G4double h, G4double yOut[], G4double yErr[])
G4ThreeVector fyMidPoint
G4ThreeVector fyFinal
G4double DistChord() const
G4ThreeVector fpInitial
void StepNoErr(const G4double tIn[8], const G4double dydx[6], G4double Step, G4double tOut[8], G4double B[3])
G4ThreeVector BfldIn
G4RKG3_Stepper(G4Mag_EqRhs *EqRhs)