Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // G4HelixExplicitEuler implementation << 26 // >> 27 // $Id: G4HelixExplicitEuler.cc,v 1.8 2007-12-10 16:29:49 gunter Exp $ >> 28 // GEANT4 tag $Name: geant4-09-04-patch-02 $ >> 29 // 27 // 30 // 28 // Helix Explicit Euler: x_1 = x_0 + helix(h) 31 // Helix Explicit Euler: x_1 = x_0 + helix(h) 29 // with helix(h) being a helix piece of lengt << 32 // with helix(h) being a helix piece of length h 30 // Most simple approach for solving linear di << 33 // most simple approach for solving linear differential equations. 31 // Take the current derivative and add it to 34 // Take the current derivative and add it to the current position. 32 // 35 // 33 // Author: W.Wander <wwc@mit.edu>, 12.09.1997 << 36 // W.Wander <wwc@mit.edu> 12/09/97 34 // ------------------------------------------- 37 // ------------------------------------------------------------------- 35 38 36 #include "G4HelixExplicitEuler.hh" 39 #include "G4HelixExplicitEuler.hh" 37 #include "G4PhysicalConstants.hh" << 38 #include "G4ThreeVector.hh" 40 #include "G4ThreeVector.hh" 39 41 40 G4HelixExplicitEuler::G4HelixExplicitEuler(G4M << 41 : G4MagHelicalStepper(EqRhs) << 42 { << 43 } << 44 << 45 G4HelixExplicitEuler::~G4HelixExplicitEuler() << 46 42 47 void G4HelixExplicitEuler::Stepper( const G4do << 43 void G4HelixExplicitEuler::Stepper( const G4double yInput[7], 48 const G4do << 44 const G4double*, 49 G4do << 45 G4double Step, 50 G4do << 46 G4double yOut[7], 51 G4do << 47 G4double yErr[]) >> 48 52 { 49 { 53 // Estimation of the Stepping Angle << 50 54 // << 51 //Estimation of the Stepping Angle >> 52 55 G4ThreeVector Bfld; 53 G4ThreeVector Bfld; 56 MagFieldEvaluate(yInput, Bfld); 54 MagFieldEvaluate(yInput, Bfld); 57 55 58 const G4int nvar = 6 ; 56 const G4int nvar = 6 ; 59 G4double yTemp[8], yIn[8] ; << 57 G4int i; >> 58 G4double yTemp[7], yIn[7] ; 60 G4ThreeVector Bfld_midpoint; 59 G4ThreeVector Bfld_midpoint; 61 << 60 // Saving yInput because yInput and yOut can be aliases for same array 62 // Saving yInput because yInput and yOut can << 61 for(i=0;i<nvar;i++) yIn[i]=yInput[i]; 63 // << 64 for(G4int i=0; i<nvar; ++i) << 65 { << 66 yIn[i] = yInput[i]; << 67 } << 68 62 69 G4double h = Step * 0.5; << 63 G4double h = Step * 0.5; 70 64 71 // Do full step and two half steps << 65 // Do full step and two half steps 72 // << 66 G4double yTemp2[7]; 73 G4double yTemp2[7]; << 67 AdvanceHelix(yIn, Bfld, h, yTemp2,yTemp); 74 AdvanceHelix(yIn, Bfld, h, yTemp2,yTemp); << 68 MagFieldEvaluate(yTemp2, Bfld_midpoint) ; 75 MagFieldEvaluate(yTemp2, Bfld_midpoint) ; << 69 AdvanceHelix(yTemp2, Bfld_midpoint, h, yOut); 76 AdvanceHelix(yTemp2, Bfld_midpoint, h, yOut) << 70 77 SetAngCurve(GetAngCurve() * 2); << 71 // Error estimation >> 72 for(i=0;i<nvar;i++) { >> 73 yErr[i] = yOut[i] - yTemp[i] ; >> 74 } 78 75 79 // Error estimation << 80 // << 81 for(G4int i=0; i<nvar; ++i) << 82 { << 83 yErr[i] = yOut[i] - yTemp[i]; << 84 } << 85 } 76 } 86 77 87 G4double G4HelixExplicitEuler::DistChord() c 78 G4double G4HelixExplicitEuler::DistChord() const 88 { 79 { 89 // Implementation : must check whether h/R > 80 // Implementation : must check whether h/R > 2 pi !! 90 // If( h/R < pi) use G4LineSection::DistL 81 // If( h/R < pi) use G4LineSection::DistLine 91 // Else DistChord=R_helix 82 // Else DistChord=R_helix 92 // 83 // 93 G4double distChord; 84 G4double distChord; 94 G4double Ang_curve=GetAngCurve(); 85 G4double Ang_curve=GetAngCurve(); 95 86 96 << 87 97 if(Ang_curve<=pi) << 88 if(Ang_curve<=pi){ 98 { << 89 distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve)); 99 distChord=GetRadHelix()*(1-std::cos(0.5*An << 90 } 100 } << 91 else 101 else if(Ang_curve<twopi) << 92 if(Ang_curve<twopi){ 102 { << 93 distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve))); 103 distChord=GetRadHelix()*(1+std::cos(0.5*(t << 94 } 104 } << 95 else{ 105 else << 96 distChord=2.*GetRadHelix(); 106 { << 97 } 107 distChord=2.*GetRadHelix(); << 108 } << 109 98 110 return distChord; 99 return distChord; >> 100 111 } 101 } 112 << 102 void 113 void G4HelixExplicitEuler::DumbStepper( const << 103 G4HelixExplicitEuler::DumbStepper( const G4double yIn[], 114 << 104 G4ThreeVector Bfld, 115 << 105 G4double h, 116 << 106 G4double yOut[]) 117 { 107 { 118 AdvanceHelix(yIn, Bfld, h, yOut); << 108 >> 109 AdvanceHelix(yIn, Bfld, h, yOut); >> 110 119 } 111 } 120 112