Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4BogackiShampine23.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 // G4BogackiShampine23 implementation
 27 //
 28 //  Bogacki-Shampine - 4 - 3(2) non-FSAL implementation
 29 //
 30 //  Implementation of the method proposed in the publication
 31 //   "A 3(2) pair of Runge - Kutta formulas"
 32 //    by P. Bogacki and L. F. Shampine,
 33 //    Appl. Math. Lett., vol. 2, no. 4, pp. 321-325, Jan. 1989.
 34 // 
 35 // The Bogacki shampine method has the following Butcher's tableau
 36 //
 37 // 0  |
 38 // 1/2|1/2
 39 // 3/4|0        3/4
 40 // 1  |2/9      1/3     4/9
 41 // -------------------
 42 //    |2/9      1/3     4/9    0
 43 //    |7/24 1/4 1/3 1/8
 44 //
 45 // Created: Somnath Banerjee, Google Summer of Code 2015, 20 May 2015
 46 // Supervision: John Apostolakis, CERN
 47 // --------------------------------------------------------------------
 48 
 49 #include "G4BogackiShampine23.hh"
 50 #include "G4LineSection.hh"
 51 #include "G4FieldUtils.hh"
 52 
 53 using namespace field_utils;
 54 
 55 G4BogackiShampine23::G4BogackiShampine23(G4EquationOfMotion* EqRhs,
 56                                          G4int integrationVariables)
 57   : G4MagIntegratorStepper(EqRhs, integrationVariables)
 58 {
 59   SetIntegrationOrder(3);
 60   SetFSAL(true);
 61 }
 62 
 63 void G4BogackiShampine23::makeStep(const G4double yInput[],
 64                                    const G4double dydx[],
 65                                    const G4double hstep,
 66                                    G4double yOutput[],
 67                                    G4double* dydxOutput,
 68                                    G4double* yError) const
 69 {
 70 
 71   G4double yTemp[G4FieldTrack::ncompSVEC];
 72   for(G4int i = GetNumberOfVariables(); i < GetNumberOfStateVariables(); ++i)
 73   { 
 74     yOutput[i] = yTemp[i] = yInput[i];
 75   }
 76 
 77   G4double ak2[G4FieldTrack::ncompSVEC],
 78            ak3[G4FieldTrack::ncompSVEC];
 79 
 80   const G4double b21 = 0.5 ,
 81                  b31 = 0., b32 = 3.0 / 4.0,
 82                  b41 = 2.0 / 9.0, b42 = 1.0 / 3.0, b43 = 4.0 / 9.0;
 83 
 84   const G4double dc1 = b41 - 7.0 / 24.0,  dc2 = b42 - 1.0 / 4.0,
 85                  dc3 = b43 - 1.0 / 3.0,   dc4 = - 1.0 / 8.0;
 86  
 87   // RightHandSide(yInput, dydx);
 88   for(G4int i = 0; i < GetNumberOfVariables(); ++i)
 89   {
 90     yTemp[i] = yInput[i] + b21 * hstep * dydx[i];
 91   }
 92     
 93   RightHandSide(yTemp, ak2);
 94   for(G4int i = 0; i < GetNumberOfVariables(); ++i)
 95   {
 96     yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
 97   }
 98 
 99   RightHandSide(yTemp, ak3);
100   for(G4int i = 0; i < GetNumberOfVariables(); ++i)
101   {
102     yOutput[i] = yInput[i] + hstep * (b41*dydx[i] + b42*ak2[i] + b43*ak3[i]);
103   }
104     
105   if ((dydxOutput != nullptr) && (yError != nullptr))
106   {
107     RightHandSide(yOutput, dydxOutput);
108     for(G4int i = 0; i < GetNumberOfVariables(); ++i)
109     {
110       yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] + 
111                            dc3 * ak3[i] + dc4 * dydxOutput[i]);
112     }
113   }
114 }
115 
116 void G4BogackiShampine23::Stepper(const G4double yInput[],
117                                   const G4double dydx[],
118                                   G4double hstep,
119                                   G4double yOutput[],
120                                   G4double yError[])
121 {
122   copy(fyIn, yInput);
123   copy(fdydx, dydx);
124   fhstep = hstep;
125 
126   makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
127 
128   copy(yOutput, fyOut);
129 }
130 
131 void G4BogackiShampine23::Stepper(const G4double yInput[],
132                                   const G4double dydx[],
133                                   G4double hstep,
134                                   G4double yOutput[],
135                                   G4double yError[],
136                                   G4double dydxOutput[])
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   copy(dydxOutput, fdydxOut);
146 }
147 
148 G4double G4BogackiShampine23::DistChord() const
149 {
150   G4double yMid[G4FieldTrack::ncompSVEC];
151   makeStep(fyIn, fdydx, fhstep / 2., yMid);
152 
153   const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
154   const G4ThreeVector mid = makeVector(yMid, Value3D::Position);
155   const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
156 
157   return G4LineSection::Distline(mid, begin, end);
158 }
159