BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
BDSIntegratorOctupole.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 "BDSIntegratorOctupole.hh"
21#include "BDSMagnetStrength.hh"
22#include "BDSStep.hh"
23#include "BDSUtilities.hh"
24
25#include "G4Mag_EqRhs.hh"
26#include "G4MagIntegratorStepper.hh"
27#include "G4ThreeVector.hh"
28
29#include "CLHEP/Units/SystemOfUnits.h"
30
31#include <cmath>
32
33BDSIntegratorOctupole::BDSIntegratorOctupole(BDSMagnetStrength const* strength,
34 G4double brho,
35 G4Mag_EqRhs* eqOfMIn):
37{
38 // B''' = d^3By/dx^3 = Brho * (1/Brho d^3By/dx^3) = Brho * k3
39 bTriplePrime = brho * (*strength)["k3"] / (CLHEP::m3*CLHEP::m);
40
41 zeroStrength = !BDS::IsFinite(bTriplePrime);
42#ifdef BDSDEBUG
43 G4cout << __METHOD_NAME__ << "B''' = " << bTriplePrime << G4endl;
44#endif
45}
46
47void BDSIntegratorOctupole::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()*bTriplePrime / momMag;
55
56 if(std::abs(kappa) < 1e-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 y3fac = y0*(y0*y0 - 3*x0*x0);
82 G4double x3fac = x0*(x0*x0 - 3*y0*y0);
83
84 // local r'' (for curvature)
85 G4ThreeVector localA;
86 localA.setX(-zp*x3fac);
87 localA.setY(-zp*y3fac);
88 localA.setZ( xp*x3fac + yp*y3fac);
89
90 localA *= kappa / 6; // 6 is actually a 3! factor.;
91
92 AdvanceChord(h,localPos,localMomUnit,localA);
93 ConvertToGlobal(localPos, localMomUnit, yOut, yErr, momMag);
94}
BDSStep ConvertToLocal(G4Step const *const step, G4bool useCurvilinear=true) const
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
virtual void AdvanceHelix(const G4double yIn[], G4double h, G4double yOut[], G4double yErr[])
Calculate the new particle coordinates.
G4double bTriplePrime
3rd differential of field.
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())