Geant4 Cross Reference |
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 // G4ConstRK4 27 // 28 // class description: 29 // 30 // G4ConstRK4 performs the integration of one step with error calculation 31 // in constant magnetic field. The integration method is the same as in 32 // ClassicalRK4. The field value is assumed constant for the step. 33 // This field evaluation is called only once per step. 34 // G4ConstRK4 can be used only for magnetic fields. 35 36 // Created: J.Apostolakis, T.Nikitina - 18.09.2008 37 // ------------------------------------------------------------------- 38 #ifndef G4CONSTRK4_HH 39 #define G4CONSTRK4_HH 40 41 #include "G4MagErrorStepper.hh" 42 #include "G4EquationOfMotion.hh" 43 #include "G4Mag_EqRhs.hh" 44 45 class G4ConstRK4 : public G4MagErrorStepper 46 { 47 public: 48 49 G4ConstRK4(G4Mag_EqRhs* EquationMotion, G4int numberOfStateVariables=8); 50 ~G4ConstRK4() override; 51 52 G4ConstRK4(const G4ConstRK4&) = delete; 53 G4ConstRK4& operator=(const G4ConstRK4&) = delete; 54 // Copy constructor and assignment operator not allowed 55 56 void Stepper( const G4double y[], 57 const G4double dydx[], 58 G4double h, 59 G4double yout[], 60 G4double yerr[] ) override; 61 void DumbStepper( const G4double yIn[], 62 const G4double dydx[], 63 G4double h, 64 G4double yOut[] ) override ; 65 G4double DistChord() const override; 66 67 inline void RightHandSideConst(const G4double y[], 68 G4double dydx[] ) const; 69 70 inline void GetConstField(const G4double y[], G4double Field[]); 71 72 G4int IntegratorOrder() const override { return 4; } 73 74 private: 75 76 G4ThreeVector fInitialPoint, fMidPoint, fFinalPoint; 77 // Data stored in order to find the chord 78 G4double *dydxm, *dydxt, *yt; // scratch space - not state 79 G4double *yInitial, *yMiddle, *dydxMid, *yOneStep; 80 G4Mag_EqRhs* fEq = nullptr; 81 G4double Field[3]; 82 }; 83 84 // Inline methods 85 86 inline void G4ConstRK4::RightHandSideConst(const G4double y[], 87 G4double dydx[] ) const 88 { 89 90 G4double momentum_mag_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5]; 91 G4double inv_momentum_magnitude = 1.0 / std::sqrt( momentum_mag_square ); 92 93 G4double cof = fEq->FCof()*inv_momentum_magnitude; 94 95 dydx[0] = y[3]*inv_momentum_magnitude; // (d/ds)x = Vx/V 96 dydx[1] = y[4]*inv_momentum_magnitude; // (d/ds)y = Vy/V 97 dydx[2] = y[5]*inv_momentum_magnitude; // (d/ds)z = Vz/V 98 99 dydx[3] = cof*(y[4]*Field[2] - y[5]*Field[1]) ; // Ax = a*(Vy*Bz - Vz*By) 100 dydx[4] = cof*(y[5]*Field[0] - y[3]*Field[2]) ; // Ay = a*(Vz*Bx - Vx*Bz) 101 dydx[5] = cof*(y[3]*Field[1] - y[4]*Field[0]) ; // Az = a*(Vx*By - Vy*Bx) 102 } 103 104 inline void G4ConstRK4::GetConstField(const G4double y[], G4double B[]) 105 { 106 G4double PositionAndTime[4]; 107 108 PositionAndTime[0] = y[0]; 109 PositionAndTime[1] = y[1]; 110 PositionAndTime[2] = y[2]; 111 // Global Time 112 PositionAndTime[3] = y[7]; 113 fEq -> GetFieldValue(PositionAndTime, B); 114 } 115 116 #endif 117