Geant4-11
G4HelixExplicitEuler.cc
Go to the documentation of this file.
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// G4HelixExplicitEuler implementation
27//
28// Helix Explicit Euler: x_1 = x_0 + helix(h)
29// with helix(h) being a helix piece of length h.
30// Most simple approach for solving linear differential equations.
31// Take the current derivative and add it to the current position.
32//
33// Author: W.Wander <wwc@mit.edu>, 12.09.1997
34// -------------------------------------------------------------------
35
38#include "G4ThreeVector.hh"
39
41 : G4MagHelicalStepper(EqRhs)
42{
43}
44
46{
47}
48
50 const G4double*,
51 G4double Step,
52 G4double yOut[7],
53 G4double yErr[] )
54{
55 // Estimation of the Stepping Angle
56 //
57 G4ThreeVector Bfld;
58 MagFieldEvaluate(yInput, Bfld);
59
60 const G4int nvar = 6 ;
61 G4double yTemp[8], yIn[8] ;
62 G4ThreeVector Bfld_midpoint;
63
64 // Saving yInput because yInput and yOut can be aliases for same array
65 //
66 for(G4int i=0; i<nvar; ++i)
67 {
68 yIn[i] = yInput[i];
69 }
70
71 G4double h = Step * 0.5;
72
73 // Do full step and two half steps
74 //
75 G4double yTemp2[7];
76 AdvanceHelix(yIn, Bfld, h, yTemp2,yTemp);
77 MagFieldEvaluate(yTemp2, Bfld_midpoint) ;
78 AdvanceHelix(yTemp2, Bfld_midpoint, h, yOut);
80
81 // Error estimation
82 //
83 for(G4int i=0; i<nvar; ++i)
84 {
85 yErr[i] = yOut[i] - yTemp[i];
86 }
87}
88
90{
91 // Implementation : must check whether h/R > 2 pi !!
92 // If( h/R < pi) use G4LineSection::DistLine
93 // Else DistChord=R_helix
94 //
95 G4double distChord;
96 G4double Ang_curve=GetAngCurve();
97
98
99 if(Ang_curve<=pi)
100 {
101 distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve));
102 }
103 else if(Ang_curve<twopi)
104 {
105 distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve)));
106 }
107 else
108 {
109 distChord=2.*GetRadHelix();
110 }
111
112 return distChord;
113}
114
116 G4ThreeVector Bfld,
117 G4double h,
118 G4double yOut[] )
119{
120 AdvanceHelix(yIn, Bfld, h, yOut);
121}
static constexpr double twopi
Definition: G4SIunits.hh:56
static constexpr double pi
Definition: G4SIunits.hh:55
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
G4double DistChord() const
G4HelixExplicitEuler(G4Mag_EqRhs *EqRhs)
void Stepper(const G4double y[], const G4double *, G4double h, G4double yout[], G4double yerr[])
void DumbStepper(const G4double y[], G4ThreeVector Bfld, G4double h, G4double yout[])
void AdvanceHelix(const G4double yIn[], G4ThreeVector Bfld, G4double h, G4double yHelix[], G4double yHelix2[]=0)
G4double GetRadHelix() const
void MagFieldEvaluate(const G4double y[], G4ThreeVector &Bfield)
void SetAngCurve(const G4double Ang)
G4double GetAngCurve() const