BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSIntegratorDecapole.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2023.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "BDSDebug.hh"
20#include "BDSIntegratorDecapole.hh"
21#include "BDSMagnetStrength.hh"
22#include "BDSStep.hh"
23#include "BDSUtilities.hh"
24
25#include "globals.hh" // geant4 types / globals
26#include "G4Mag_EqRhs.hh"
27#include "G4MagIntegratorStepper.hh"
28#include "G4ThreeVector.hh"
29
30#include "CLHEP/Units/SystemOfUnits.h"
31
32#include <cmath>
33
35 G4double brho,
36 G4Mag_EqRhs* eqOfMIn):
38{
39 // B'''' = d^4By/dx^4 = Brho * (1/Brho d^4By/dx^4) = Brho * k4
40 bQuadruplePrime = brho * (*strength)["k4"] / (CLHEP::m3*CLHEP::m2);
41 zeroStrength = !BDS::IsFinite(bQuadruplePrime);
42#ifdef BDSDEBUG
43 G4cout << __METHOD_NAME__ << "B'''' = " << bQuadruplePrime << G4endl;
44#endif
45}
46
47void BDSIntegratorDecapole::AdvanceHelix(const G4double yIn[],
48 G4double h,
49 G4double yOut[],
50 G4double yErr[])
51{
52 G4ThreeVector mom = G4ThreeVector(yIn[3], yIn[4], yIn[5]);
53 G4double momMag = mom.mag();
54 G4double kappa = -eqOfM->FCof()*bQuadruplePrime / momMag;
55
56 if(std::abs(kappa)<1.e-20)
57 {
58 AdvanceDriftMag(yIn, h, yOut, yErr);
59 SetDistChord(0);
60 return;
61 }
62
63 G4ThreeVector pos = G4ThreeVector(yIn[0], yIn[1], yIn[2]);
64 BDSStep localPosMom = ConvertToLocal(pos, mom, h, false);
65 G4ThreeVector localPos = localPosMom.PreStepPoint();
66 G4ThreeVector localMom = localPosMom.PostStepPoint();
67 G4ThreeVector localMomUnit = localMom.unit();
68
69 G4double x0 = localPos.x();
70 G4double y0 = localPos.y();
71
72 // Evaluate field at the approximate midpoint of the step.
73 const G4double halfH = 0.5*h;
74 x0 = x0 + localMomUnit.x()*halfH;
75 y0 = y0 + localMomUnit.y()*halfH;
76
77 G4double xp = localMomUnit.x();
78 G4double yp = localMomUnit.y();
79 G4double zp = localMomUnit.z();
80
81 G4double xsq = std::pow(x0, 2);
82 G4double ysq = std::pow(y0, 2);
83 G4double x02My02 = (xsq - ysq);
84 G4double Bx = 4.0*x0*y0 * (x02My02);
85 G4double By = std::pow(x0,4) - 6.0 * xsq * ysq + std::pow(y0,4);
86
87 // local r'' (for curvature)
88 G4ThreeVector localA;
89 localA.setX( zp*By);
90 localA.setY(-zp*Bx);
91 localA.setZ( xp*By - yp*Bx);
92
93 localA *= kappa / 24; // 24 is actually a 4! factor.;
94
95 AdvanceChord(h,localPos,localMomUnit,localA);
96 ConvertToGlobal(localPos, localMomUnit, yOut, yErr, momMag);
97}
BDSStep ConvertToLocal(G4Step const *const step, G4bool useCurvilinear=true) const
G4double bQuadruplePrime
4th derivative of magnetic field
virtual void AdvanceHelix(const G4double yIn[], G4double h, G4double yDec[], G4double yErr[])
Calculate the new particle coordinates.
BDSIntegratorDecapole()=delete
Private default constructor to enforce use of supplied constructor.
void AdvanceDriftMag(const G4double yIn[], const G4double h, G4double yOut[], G4double yErr[]) const
Error array [6] all 0.
Common functionality for Euler integrators.
void AdvanceChord(const G4double h, G4ThreeVector &localPos, G4ThreeVector &localMom, const G4ThreeVector &localA)
G4Mag_EqRhs * eqOfM
void SetDistChord(G4double distChordIn)
Setter for distChord to private member.
void ConvertToGlobal(const G4ThreeVector &localPos, const G4ThreeVector &localMom, G4double yOut[], G4double yErr[], const G4double momScaling=1.0)
scaling of momentum in case localMom is a unit vector
Efficient storage of magnet strengths.
A simple class to represent the positions of a step.
Definition: BDSStep.hh:33
G4ThreeVector PostStepPoint() const
Accessor.
Definition: BDSStep.hh:43
G4ThreeVector PreStepPoint() const
Accessor.
Definition: BDSStep.hh:42
G4bool IsFinite(G4double value, G4double tolerance=std::numeric_limits< double >::epsilon())