Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4RK547FEq3.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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     {
 99       yTemp[i] = yInput[i] + hstep * b21 * dydx[i];
100     }
101 
102     RightHandSide(yTemp, ak2);
103     for(G4int i = 0; i < GetNumberOfVariables(); ++i)
104     {
105       yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
106     }
107 
108     RightHandSide(yTemp, ak3);
109     for(G4int i = 0;i < GetNumberOfVariables(); ++i)
110     {
111       yTemp[i] = yInput[i] + hstep * (b41 * dydx[i] + b42 * ak2[i] +
112                                       b43 * ak3[i]);
113     }
114 
115     RightHandSide(yTemp, ak4);
116     for(G4int i = 0; i < GetNumberOfVariables(); ++i)
117     {
118       yTemp[i] = yInput[i] + hstep * (b51 * dydx[i] + b52 * ak2[i] +
119                                       b53 * ak3[i] + b54 * ak4[i]);
120     }
121 
122     RightHandSide(yTemp, ak5);
123     for(G4int i = 0; i < GetNumberOfVariables(); ++i)
124     {
125       yTemp[i] = yInput[i] + hstep * (b61 * dydx[i] + b62 * ak2[i] +
126                                       b63 * ak3[i] + b64 * ak4[i] +
127                                       b65 * ak5[i]);
128     }
129 
130     RightHandSide(yTemp, ak6);
131     for(G4int i = 0; i < GetNumberOfVariables(); ++i)
132     {
133       yOutput[i] = yInput[i] + hstep * (b71 * dydx[i] + b72 * ak2[i] +
134                                         b73 * ak3[i] + b74 * ak4[i] +
135                                         b75 * ak5[i] + b76 * ak6[i]);
136     }
137     if ((dydxOutput != nullptr) && (yError != nullptr))
138     {
139         RightHandSide(yOutput, dydxOutput);
140         for(G4int i = 0; i < GetNumberOfVariables(); ++i)
141         {
142           yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] + dc3 * ak3[i] +
143                                dc4 * ak4[i] + dc5 * ak5[i] + dc6 * ak6[i] +
144                                dc7 * dydxOutput[i]);
145         }
146     }
147 }
148 
149 void G4RK547FEq3::Stepper( const G4double yInput[],
150                            const G4double dydx[],
151                                  G4double hstep,
152                                  G4double yOutput[],
153                                  G4double yError[] )
154 {
155     copy(fyIn, yInput);
156     copy(fdydx, dydx);
157     fhstep = hstep;
158 
159     makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
160 
161     copy(yOutput, fyOut);
162 }
163 
164 void G4RK547FEq3::Stepper( const G4double yInput[],
165                            const G4double dydx[],
166                                  G4double hstep,
167                                  G4double yOutput[],
168                                  G4double yError[],
169                                  G4double dydxOutput[] )
170 {
171     copy(fyIn, yInput);
172     copy(fdydx, dydx);
173     fhstep = hstep;
174 
175     makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
176 
177     copy(yOutput, fyOut);
178     copy(dydxOutput, fdydxOut);
179 }
180 
181 G4double G4RK547FEq3::DistChord() const
182 {
183     G4double yMid[G4FieldTrack::ncompSVEC];
184     makeStep(fyIn, fdydx, fhstep / 2., yMid);
185 
186     const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
187     const G4ThreeVector mid = makeVector(yMid, Value3D::Position);
188     const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
189 
190     return G4LineSection::Distline(mid, begin, end);
191 }
192