ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RK547FEq1.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 // G4RK547FEq1 implementation
27 //
28 // The Butcher table of the Higham & Hall 5(4)7 method is:
29 //
30 // 0 |
31 // 2/9 | 2/9
32 // 1/3 | 1/12 1/4
33 // 1/2 | 1/8 0 3/8
34 // 3/5 | 91/500 -27/100 78/125 8/125
35 // 1 | -11/20 27/20 12/5 -36/5 5
36 // 1 | 1/12 0 27/32 -4/3 125/96 5/48
37 //----------------------------------------------------------------------------
38 // 1/12 0 27/32 -4/3 125/96 5/48 0
39 // 2/15 0 27/80 -2/15 25/48 1/24 1/10
40 //
41 // Author: Dmitry Sorokin, Google Summer of Code 2017
42 // Supervision: John Apostolakis, CERN
43 // --------------------------------------------------------------------
44 
45 #include "G4RK547FEq1.hh"
46 #include "G4LineSection.hh"
47 #include "G4FieldUtils.hh"
48 
49 using namespace field_utils;
50 
51 G4RK547FEq1::G4RK547FEq1(G4EquationOfMotion* EqRhs, G4int integrationVariables)
52  : G4MagIntegratorStepper(EqRhs, integrationVariables)
53 {
54 }
55 
56 void G4RK547FEq1::makeStep( const G4double yInput[],
57  const G4double dydx[],
58  const G4double hstep,
59  G4double yOutput[],
60  G4double* dydxOutput,
61  G4double* yError ) const
62 {
63  G4double yTemp[G4FieldTrack::ncompSVEC];
64  for (G4int i=GetNumberOfVariables(); i<GetNumberOfStateVariables(); ++i)
65  {
66  yOutput[i] = yTemp[i] = yInput[i];
67  }
68 
69  G4double ak2[G4FieldTrack::ncompSVEC],
70  ak3[G4FieldTrack::ncompSVEC],
71  ak4[G4FieldTrack::ncompSVEC],
72  ak5[G4FieldTrack::ncompSVEC],
73  ak6[G4FieldTrack::ncompSVEC];
74 
75  const G4double b21 = 2./9.,
76  b31 = 1./12., b32 = 1./4.,
77  b41 = 1./8., b42 = 0., b43 = 3./8.,
78  b51 = 91./500., b52 = -27./100.,
79  b53 = 78./125., b54 = 8./125.,
80  b61 = -11./20., b62 = 27./20., b63 = 12./5.,
81  b64 = -36./5., b65 = 5.,
82  b71 = 1./12., b72 = 0., b73 = 27./32.,
83  b74 = -4./3., b75 = 125./96., b76 = 5./48.;
84 
85  const G4double dc1 = b71 - 2./15.,
86  dc2 = b72 - 0.,
87  dc3 = b73 - 27./80.,
88  dc4 = b74 + 2./15.,
89  dc5 = b75 - 25./48.,
90  dc6 = b76 - 1./24.,
91  dc7 = 0. - 1./10.;
92 
93  // RightHandSide(yInput, dydx);
94  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
95  yTemp[i] = yInput[i] + hstep * b21 * dydx[i];
96 
97  RightHandSide(yTemp, ak2);
98  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
99  yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
100 
101  RightHandSide(yTemp, ak3);
102  for(G4int i = 0;i < GetNumberOfVariables(); ++i)
103  yTemp[i] = yInput[i] + hstep * (b41 * dydx[i] + b42 * ak2[i] +
104  b43 * ak3[i]);
105 
106  RightHandSide(yTemp, ak4);
107  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
108  yTemp[i] = yInput[i] + hstep * (b51 * dydx[i] + b52 * ak2[i] +
109  b53 * ak3[i] + b54 * ak4[i]);
110 
111  RightHandSide(yTemp, ak5);
112  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
113  yTemp[i] = yInput[i] + hstep * (b61 * dydx[i] + b62 * ak2[i] +
114  b63 * ak3[i] + b64 * ak4[i] +
115  b65 * ak5[i]);
116 
117  RightHandSide(yTemp, ak6);
118  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
119  yOutput[i] = yInput[i] + hstep * (b71 * dydx[i] + b72 * ak2[i] +
120  b73 * ak3[i] + b74 * ak4[i] +
121  b75 * ak5[i] + b76 * ak6[i]);
122  if (dydxOutput && yError)
123  {
124  RightHandSide(yOutput, dydxOutput);
125  for(G4int i = 0; i < GetNumberOfVariables(); ++i)
126  yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] + dc3 * ak3[i] +
127  dc4 * ak4[i] + dc5 * ak5[i] + dc6 * ak6[i] +
128  dc7 * dydxOutput[i]);
129  }
130 }
131 
132 void G4RK547FEq1::Stepper( const G4double yInput[],
133  const G4double dydx[],
134  G4double hstep,
135  G4double yOutput[],
136  G4double yError[] )
137 {
138  copy(fyIn, yInput);
139  copy(fdydx, dydx);
140  fhstep = hstep;
141 
142  makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
143 
144  copy(yOutput, fyOut);
145 }
146 
147 void G4RK547FEq1::Stepper( const G4double yInput[],
148  const G4double dydx[],
149  G4double hstep,
150  G4double yOutput[],
151  G4double yError[],
152  G4double dydxOutput[] )
153 {
154  copy(fyIn, yInput);
155  copy(fdydx, dydx);
156  fhstep = hstep;
157 
158  makeStep(fyIn,fdydx, fhstep, fyOut, fdydxOut, yError);
159 
160  copy(yOutput, fyOut);
161  copy(dydxOutput, fdydxOut);
162 }
163 
164 G4double G4RK547FEq1::DistChord() const
165 {
166  G4double yMid[G4FieldTrack::ncompSVEC];
167  makeStep(fyIn, fdydx, fhstep / 2., yMid);
168 
169  const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
170  const G4ThreeVector mid = makeVector(yMid, Value3D::Position);
171  const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
172 
173  return G4LineSection::Distline(mid, begin, end);
174 }