Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4RK547FEq1.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 // The Butcher table of the Higham & Hall 5(4)7 method is:
27 //
28 // 0 |
29 // 2/9 | 2/9
30 // 1/3 | 1/12 1/4
31 // 1/2 | 1/8 0 3/8
32 // 3/5 | 91/500 -27/100 78/125 8/125
33 // 1 | -11/20 27/20 12/5 -36/5 5
34 // 1 | 1/12 0 27/32 -4/3 125/96 5/48
35 //----------------------------------------------------------------------------
36 // 1/12 0 27/32 -4/3 125/96 5/48 0
37 // 2/15 0 27/80 -2/15 25/48 1/24 1/10
38 
39 #include "G4RK547FEq1.hh"
40 #include "G4LineSection.hh"
41 #include "G4FieldUtils.hh"
42 
43 using namespace field_utils;
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 
54 G4RK547FEq1::G4RK547FEq1(G4EquationOfMotion* EqRhs, G4int integrationVariables)
55  : G4MagIntegratorStepper(EqRhs, integrationVariables)
56 {
57 }
58 
60  const G4double yInput[],
61  const G4double dydx[],
62  const G4double hstep,
63  G4double yOutput[],
64  G4double* dydxOutput,
65  G4double* yError) const
66 {
68  for (int i = GetNumberOfVariables(); i < GetNumberOfStateVariables(); ++i){
69  yOutput[i] = yTemp[i] = yInput[i];
70  }
71 
77 
78  const G4double
79  b21 = 2./9.,
80  b31 = 1./12., b32 = 1./4.,
81  b41 = 1./8., b42 = 0., b43 = 3./8.,
82  b51 = 91./500., b52 = -27./100., b53 = 78./125., b54 = 8./125.,
83  b61 = -11./20., b62 = 27./20., b63 = 12./5.,
84  b64 = -36./5., b65 = 5.,
85  b71 = 1./12., b72 = 0., b73 = 27./32.,
86  b74 = -4./3., b75 = 125./96., b76 = 5./48.;
87 
88  const G4double
89  dc1 = b71 - 2./15.,
90  dc2 = b72 - 0.,
91  dc3 = b73 - 27./80.,
92  dc4 = b74 + 2./15.,
93  dc5 = b75 - 25./48.,
94  dc6 = b76 - 1./24.,
95  dc7 = 0. - 1./10.;
96 
97  //RightHandSide(yInput, dydx);
98  for(int i = 0; i < GetNumberOfVariables(); ++i)
99  yTemp[i] = yInput[i] + hstep * b21 * dydx[i];
100 
101  RightHandSide(yTemp, ak2);
102  for(int i = 0; i < GetNumberOfVariables(); ++i)
103  yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
104 
105  RightHandSide(yTemp, ak3);
106  for(int i = 0;i < GetNumberOfVariables(); ++i)
107  yTemp[i] = yInput[i] + hstep * (b41 * dydx[i] + b42 * ak2[i] +
108  b43 * ak3[i]);
109 
110  RightHandSide(yTemp, ak4);
111  for(int i = 0; i < GetNumberOfVariables(); ++i)
112  yTemp[i] = yInput[i] + hstep * (b51 * dydx[i] + b52 * ak2[i] +
113  b53 * ak3[i] + b54 * ak4[i]);
114 
115  RightHandSide(yTemp, ak5);
116  for(int i = 0; i < GetNumberOfVariables(); ++i)
117  yTemp[i] = yInput[i] + hstep * (b61 * dydx[i] + b62 * ak2[i] +
118  b63 * ak3[i] + b64 * ak4[i] +
119  b65 * ak5[i]);
120 
121  RightHandSide(yTemp, ak6);
122  for(int i = 0; i < GetNumberOfVariables(); ++i)
123  yOutput[i] = yInput[i] + hstep * (b71 * dydx[i] + b72 * ak2[i] +
124  b73 * ak3[i] + b74 * ak4[i] +
125  b75 * ak5[i] + b76 * ak6[i]);
126 
127  if (dydxOutput && yError) {
128  RightHandSide(yOutput, dydxOutput);
129  for(int i = 0; i < GetNumberOfVariables(); ++i)
130  yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] + dc3 * ak3[i] +
131  dc4 * ak4[i] + dc5 * ak5[i] + dc6 * ak6[i] +
132  dc7 * dydxOutput[i]);
133  }
134 }
135 
137  const G4double yInput[],
138  const G4double dydx[],
139  G4double hstep,
140  G4double yOutput[],
141  G4double yError[])
142 {
143  copyArray(fyIn, yInput);
144  copyArray(fdydx, dydx);
145  fhstep = hstep;
146 
147  makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
148 
149  copyArray(yOutput, fyOut);
150 }
151 
153  const G4double yInput[],
154  const G4double dydx[],
155  G4double hstep,
156  G4double yOutput[],
157  G4double yError[],
158  G4double dydxOutput[])
159 {
160  copyArray(fyIn, yInput);
161  copyArray(fdydx, dydx);
162  fhstep = hstep;
163 
164  makeStep(fyIn,fdydx, fhstep, fyOut, fdydxOut, yError);
165 
166  copyArray(yOutput, fyOut);
167  copyArray(dydxOutput, fdydxOut);
168 }
169 
171 {
173  makeStep(fyIn, fdydx, fhstep / 2., yMid);
174 
175  const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
176  const G4ThreeVector mid = makeVector(yMid, Value3D::Position);
177  const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
178 
179  return G4LineSection::Distline(mid, begin, end);
180 }
G4double fdydx[G4FieldTrack::ncompSVEC]
Definition: G4RK547FEq1.hh:80
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
G4ThreeVector makeVector(const ArrayType &array, Value3D value)
G4double fyIn[G4FieldTrack::ncompSVEC]
Definition: G4RK547FEq1.hh:80
G4double fdydxOut[G4FieldTrack::ncompSVEC]
Definition: G4RK547FEq1.hh:80
void makeStep(const G4double yInput[], const G4double dydx[], const G4double hstep, G4double yOutput[], G4double *dydxOutput=nullptr, G4double *yError=nullptr) const
Definition: G4RK547FEq1.cc:59
virtual void Stepper(const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]) override
Definition: G4RK547FEq1.cc:136
G4double fhstep
Definition: G4RK547FEq1.hh:84
double G4double
Definition: G4Types.hh:76
G4int GetNumberOfStateVariables() const
void RightHandSide(const double y[], double dydx[]) const
virtual G4double DistChord() const override
Definition: G4RK547FEq1.cc:170
int G4int
Definition: G4Types.hh:78
G4double fyOut[G4FieldTrack::ncompSVEC]
Definition: G4RK547FEq1.hh:80
static const G4double ak2
G4int GetNumberOfVariables() const
G4RK547FEq1(G4EquationOfMotion *EqRhs, G4int integrationVariables=6)
Definition: G4RK547FEq1.cc:54