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 // G4RK547FEq2 implementation 27 // 28 // The Butcher table of the Higham & Hall 5(4)7 method is: 29 // 30 // 0 | 31 // 2/13 | 2/13 32 // 2/13 | 3/52 9/52 33 // 5/9 | 12955/26244 -15925/8748 12350/6561 34 // 3/4 | -10383/52480 13923/10496 -176553/199424 505197/997120 35 // 1 | 1403/7236 -429/268 733330/309339 -7884/8911 104960/113967 36 // 1 | 181/2700 0 656903/1846800 19683/106400 34112/110565 67/800 37 //---------------------------------------------------------------------------------------------------------------------- 38 // 181/2700 0 656903/1846800 19683/106400 34112/110565 67/800 0 39 // 11377/154575 0 35378291/105729300 343359/1522850 535952/1947645 134/17175 1/12 40 // 41 // Author: Dmitry Sorokin, Google Summer of Code 2017 42 // Supervision: John Apostolakis, CERN 43 // -------------------------------------------------------------------- 44 45 #include "G4RK547FEq2.hh" 46 #include "G4LineSection.hh" 47 #include "G4FieldUtils.hh" 48 49 using namespace field_utils; 50 51 G4RK547FEq2::G4RK547FEq2(G4EquationOfMotion* EqRhs, G4int integrationVariables) 52 : G4MagIntegratorStepper(EqRhs, integrationVariables) 53 { 54 } 55 56 void G4RK547FEq2::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 = 2./13., 76 b31 = 3./52., b32 = 9./52., 77 b41 = 12955./26244., b42 = -15925./8748., 78 b43 = 12350./6561., 79 b51 = -10383./52480., b52 = 13923./10496., 80 b53 = -176553./199424., b54 = 505197./997120., 81 b61 = 1403./7236., b62 = -429./268., b63 = 733330./309339., 82 b64 = -7884./8911., b65 = 104960./113967., 83 b71 = 181./2700., b72 = 0., b73 = 656903./1846800., 84 b74 = 19683./106400., b75 = 34112./110565., 85 b76 = 67./800.; 86 87 const G4double dc1 = b71 - 11377./154575., 88 dc2 = b72 - 0., 89 dc3 = b73 - 35378291./105729300., 90 dc4 = b74 - 343359./1522850., 91 dc5 = b75 - 535952./1947645., 92 dc6 = b76 - 134./17175., 93 dc7 = 0. - 1./12.; 94 95 // RightHandSide(yInput, dydx); 96 for(G4int i = 0; i < GetNumberOfVariables(); ++i) 97 { 98 yTemp[i] = yInput[i] + hstep * b21 * dydx[i]; 99 } 100 101 RightHandSide(yTemp, ak2); 102 for(G4int i = 0; i < GetNumberOfVariables(); ++i) 103 { 104 yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]); 105 } 106 107 RightHandSide(yTemp, ak3); 108 for(G4int i = 0;i < GetNumberOfVariables(); ++i) 109 { 110 yTemp[i] = yInput[i] + hstep * (b41 * dydx[i] + b42 * ak2[i] + 111 b43 * ak3[i]); 112 } 113 114 RightHandSide(yTemp, ak4); 115 for(G4int i = 0; i < GetNumberOfVariables(); ++i) 116 { 117 yTemp[i] = yInput[i] + hstep * (b51 * dydx[i] + b52 * ak2[i] + 118 b53 * ak3[i] + b54 * ak4[i]); 119 } 120 121 RightHandSide(yTemp, ak5); 122 for(G4int i = 0; i < GetNumberOfVariables(); ++i) 123 { 124 yTemp[i] = yInput[i] + hstep * (b61 * dydx[i] + b62 * ak2[i] + 125 b63 * ak3[i] + b64 * ak4[i] + 126 b65 * ak5[i]); 127 } 128 129 RightHandSide(yTemp, ak6); 130 for(G4int i = 0; i < GetNumberOfVariables(); ++i) 131 { 132 yOutput[i] = yInput[i] + hstep * (b71 * dydx[i] + b72 * ak2[i] + 133 b73 * ak3[i] + b74 * ak4[i] + 134 b75 * ak5[i] + b76 * ak6[i]); 135 } 136 if ((dydxOutput != nullptr) && (yError != nullptr)) 137 { 138 RightHandSide(yOutput, dydxOutput); 139 for(G4int i = 0; i < GetNumberOfVariables(); ++i) 140 { 141 yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] + dc3 * ak3[i] + 142 dc4 * ak4[i] + dc5 * ak5[i] + dc6 * ak6[i] + 143 dc7 * dydxOutput[i]); 144 } 145 } 146 } 147 148 void G4RK547FEq2::Stepper( const G4double yInput[], 149 const G4double dydx[], 150 G4double hstep, 151 G4double yOutput[], 152 G4double yError[] ) 153 { 154 copy(fyIn, yInput); 155 copy(fdydx, dydx); 156 fhstep = hstep; 157 158 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError); 159 160 copy(yOutput, fyOut); 161 } 162 163 void G4RK547FEq2::Stepper( const G4double yInput[], 164 const G4double dydx[], 165 G4double hstep, 166 G4double yOutput[], 167 G4double yError[], 168 G4double dydxOutput[] ) 169 { 170 copy(fyIn, yInput); 171 copy(fdydx, dydx); 172 fhstep = hstep; 173 174 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError); 175 176 copy(yOutput, fyOut); 177 copy(dydxOutput, fdydxOut); 178 } 179 180 G4double G4RK547FEq2::DistChord() const 181 { 182 G4double yMid[G4FieldTrack::ncompSVEC]; 183 makeStep(fyIn, fdydx, fhstep / 2., yMid); 184 185 const G4ThreeVector begin = makeVector(fyIn, Value3D::Position); 186 const G4ThreeVector mid = makeVector(yMid, Value3D::Position); 187 const G4ThreeVector end = makeVector(fyOut, Value3D::Position); 188 189 return G4LineSection::Distline(mid, begin, end); 190 } 191