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 // 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