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