Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4RKG3_Stepper.cc
이 파일의 문서화 페이지로 가기
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 //
27 // $Id: G4RKG3_Stepper.cc 68055 2013-03-13 14:43:28Z gcosmo $
28 //
29 // -------------------------------------------------------------------
30 
31 #include "G4RKG3_Stepper.hh"
32 #include "G4LineSection.hh"
33 #include "G4Mag_EqRhs.hh"
34 
36  : G4MagIntegratorStepper(EqRhs,6), hStep(0.)
37 {
38 }
39 
41 {
42 }
43 
44 void 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  G4int i;
53  G4double by15 = 1. / 15. ; // was 0.066666666 ;
54 
55  G4double yTemp[8], dydxTemp[6], yIn[8] ;
56  // Saving yInput because yInput and yOut can be aliases for same array
57  for(i=0;i<nvar;i++) yIn[i]=yInput[i];
58  yIn[6] = yInput[6];
59  yIn[7] = yInput[7];
60  G4double h = Step * 0.5;
61  hStep=Step;
62  // Do two half steps
63 
64  StepNoErr(yIn, dydx,h, yTemp,B) ;
65 
66  //Store Bfld for DistChord Calculation
67  for(i=0;i<3;i++)BfldIn[i]=B[i];
68 
69  // RightHandSide(yTemp,dydxTemp) ;
70 
71  GetEquationOfMotion()->EvaluateRhsGivenB(yTemp,B,dydxTemp) ;
72  StepNoErr(yTemp,dydxTemp,h,yOut,B);
73 
74  // Store midpoint, chord calculation
75 
76  fyMidPoint = G4ThreeVector( yTemp[0], yTemp[1], yTemp[2]);
77 
78  // Do a full Step
79 
80  h *= 2 ;
81  StepNoErr(yIn,dydx,h,yTemp,B);
82  for(i=0;i<nvar;i++)
83  {
84  yErr[i] = yOut[i] - yTemp[i] ;
85  yOut[i] += yErr[i]*by15 ; // Provides 5th order of accuracy
86  }
87 
88  //Store values for DistChord method
89 
90  fyInitial = G4ThreeVector( yIn[0], yIn[1], yIn[2]);
91  fpInitial = G4ThreeVector( yIn[3], yIn[4], yIn[5]);
92  fyFinal = G4ThreeVector( yOut[0], yOut[1], yOut[2]);
93 
94  // NormaliseTangentVector( yOut ); // Deleted
95 }
96 
97 // ---------------------------------------------------------------------------
98 
99 // Integrator for RK from G3 with evaluation of error in solution and delta
100 // geometry based on naive similarity with the case of uniform magnetic field.
101 // B1[3] is input and is the first magnetic field values
102 // B2[3] is output and is the final magnetic field values.
103 
105  const G4double*,
106  G4double,
107  G4double*,
108  G4double&,
109  G4double&,
110  const G4double*,
111  G4double* )
112 
113 {
114  G4Exception("G4RKG3_Stepper::StepWithEst()", "GeomField0001",
115  FatalException, "Method no longer used.");
116 }
117 
118 // -----------------------------------------------------------------
119 
120 
121 // Integrator RK Stepper from G3 with only two field evaluation per Step.
122 // It is used in propagation initial Step by small substeps after solution
123 // error and delta geometry considerations. B[3] is magnetic field which
124 // is passed from substep to substep.
125 
127  const G4double dydx[6],
128  G4double Step,
129  G4double tOut[8],
130  G4double B[3] ) // const
131 
132 {
133 
134  // Copy and edit the routine above, to delete alpha2, beta2, ...
135  G4double K1[7],K2[7],K3[7],K4[7] ;
136  G4double tTemp[8], yderiv[6] ;
137 
138  // Need Momentum value to give correct values to the coefficients in equation
139  // Integration on unit velocity, but tIn[3,4,5] is momentum
140  G4double mom,inverse_mom;
141  G4int i ;
142  const G4double c1=0.5,c2=0.125,c3=1./6.;
143 
144  // GetEquationOfMotion()->EvaluateRhsReturnB(tIn,dydx,B1) ;
145  // Correction for momentum not a velocity
146  // Need the protection !!! must be not zero
147  mom=std::sqrt(tIn[3]*tIn[3]+tIn[4]*tIn[4]+tIn[5]*tIn[5]);
148  inverse_mom=1./mom;
149  for(i=0;i<3;i++)
150  {
151  K1[i] = Step * dydx[i+3]*inverse_mom;
152  tTemp[i] = tIn[i] + Step*(c1*tIn[i+3]*inverse_mom + c2*K1[i]) ;
153  tTemp[i+3] = tIn[i+3] + c1*K1[i]*mom ;
154 
155  }
156 
157  GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
158 
159 
160  for(i=0;i<3;i++)
161  {
162  K2[i] = Step * yderiv[i+3]*inverse_mom;
163  tTemp[i+3] = tIn[i+3] + c1*K2[i]*mom ;
164  }
165 
166  // Given B, calculate yderiv !
167  GetEquationOfMotion()->EvaluateRhsGivenB(tTemp,B,yderiv) ;
168 
169  for(i=0;i<3;i++)
170  {
171  K3[i] = Step * yderiv[i+3]*inverse_mom;
172  tTemp[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom + c1*K3[i]) ;
173  tTemp[i+3] = tIn[i+3] + K3[i]*mom ;
174  }
175 
176 
177  // Calculates y-deriv(atives) & returns B too!
178  GetEquationOfMotion()->EvaluateRhsReturnB(tTemp,yderiv,B) ;
179 
180 
181  for(i=0;i<3;i++) // Output trajectory vector
182  {
183  K4[i] = Step * yderiv[i+3]*inverse_mom;
184  tOut[i] = tIn[i] + Step*(tIn[i+3]*inverse_mom+ (K1[i] + K2[i] + K3[i])*c3) ;
185  tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2[i] + 2*K3[i] +K4[i])*c3 ;
186  }
187  tOut[6] = tIn[6];
188  tOut[7] = tIn[7];
189  // NormaliseTangentVector( tOut );
190 
191 
192 }
193 
194 
195 // ---------------------------------------------------------------------------
196 
198  {
199  // Soon: must check whether h/R > 2 pi !!
200  // Method below is good only for < 2 pi
201  G4double distChord,distLine;
202 
203  if (fyInitial != fyFinal) {
205 
206  distChord = distLine;
207  }else{
208  distChord = (fyMidPoint-fyInitial).mag();
209  }
210 
211 
212  return distChord;
213 
214  }
215 
G4ThreeVector fyFinal
CLHEP::Hep3Vector G4ThreeVector
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
G4EquationOfMotion * GetEquationOfMotion()
void Stepper(const G4double yIn[], const G4double dydx[], G4double h, G4double yOut[], G4double yErr[])
G4double DistChord() const
G4ThreeVector fyMidPoint
G4ThreeVector fyInitial
G4ThreeVector BfldIn
double G4double
Definition: G4Types.hh:76
TCanvas * c2
Definition: plot_hist.C:75
void EvaluateRhsReturnB(const G4double y[], G4double dydx[], G4double Field[]) const
void StepNoErr(const G4double tIn[8], const G4double dydx[6], G4double Step, G4double tOut[8], G4double B[3])
G4RKG3_Stepper(G4Mag_EqRhs *EqRhs)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
int G4int
Definition: G4Types.hh:78
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])
G4ThreeVector fpInitial
virtual void EvaluateRhsGivenB(const G4double y[], const G4double B[3], G4double dydx[]) const =0
double B(double temperature)
Definition: Step.hh:41