74 int nx = 20,
int cx = 4);
82 bool backBeam,
bool backTrack):
83 Mapper(beamline, reference, backBeam, backTrack)
161 int order = field.
order();
181 static const Map ident;
182 Map translate = co + ident;
186 Series pz_trans =
sqrt((pt * pt - px * px - py * py).substitute(translate), order);
196 for(
int i = 0; i <
PSdim; i += 2) {
201 static const Vector zeroV;
202 Vector new_co = co + implicitInt4(zeroV, JgradH, length, 0.5 * length);
232 int order = field.
order();
238 if(stepsize != 0.0) {
239 int nst = (int)
ceil(length / stepsize);
242 double d_ell = length / nSlices;
258 static const Series de2 = de * de;
259 static const Series pxy2 = px * px + py * py;
260 static const Map ident;
263 double kin2 = kin * kin;
285 for(
int slice = 0; slice < nSlices; ++slice) {
287 double p1 = 1.0 + co[
PT];
288 Series de12 = p1 * p1 + 2.0 * p1 * de + de2;
289 Series H_trans = beta_inv *
sqrt(de12 + kin2, order)
291 - 2.0 * (co[
PX] * px + co[
PY] * py) - pxy2, order)
295 for(
int i = 0; i <
PSdim; i += 2) {
300 static const Vector zeroV;
301 Vector new_co = co + implicitInt4(zeroV, JgradH, d_ell, 0.5 * d_ell);
303 Vector mid_co = (co + new_co) / 2.0;
304 p1 = 1.0 + mid_co[
PT];
305 de12 = p1 * p1 + 2.0 * p1 * de + de2;
306 H_trans = beta_inv *
sqrt(de12 + kin2, order)
307 -
sqrt(de12 - mid_co[
PX] * mid_co[
PX] - mid_co[
PY] * mid_co[
PY]
308 - 2.0 * (mid_co[
PX] * px + mid_co[
PY] * py) - pxy2, order)
346 Series speed = (
c * pt) /
sqrt(pt * pt + kin * kin);
368 int order = field.
order();
374 if(stepsize != 0.0) {
375 int nst = (int)
ceil(length / stepsize);
378 double d_ell = length / nSlices;
394 static const Series de2 = de * de;
395 static const Series pxy2 = px * px + py * py;
396 static const Map ident;
400 double kin2 = kin * kin;
401 Series hx1 = (1.0 + h * x);
415 for(
int slice = 0; slice < nSlices; ++slice) {
417 double p1 = 1.0 + co[
PT];
418 Series de12 = p1 * p1 + 2.0 * p1 * de + de2;
419 Series H_trans = beta_inv *
sqrt(de12 + kin2, order)
422 - 2.0 * (co[
PX] * px + co[
PY] * py) - pxy2, order)
426 for(
int i = 0; i <
PSdim; i += 2) {
431 static const Vector zeroV;
432 Vector new_co = co + implicitInt4(zeroV, JgradH, d_ell, 0.5 * d_ell);
434 Vector mid_co = (co + new_co) / 2.0;
435 p1 = 1.0 + mid_co[
PT];
436 de12 = p1 * p1 + 2.0 * p1 * de + de2;
437 H_trans = beta_inv *
sqrt(de12 + kin2, order)
438 - (hx1 + h * mid_co[
X])
439 *
sqrt(de12 - mid_co[
PX] * mid_co[
PX] - mid_co[
PY] * mid_co[
PY]
440 - 2.0 * (mid_co[
PX] * px + mid_co[
PY] * py) - pxy2, order)
468 double Ex = scale * sep.
getEx();
469 double Ey = scale * sep.
getEy();
500 Series pz =
sqrt(pt * pt - px * px - py * py);
507 Series xt = C * itsMap[
X] + S * itsMap[
Y];
508 Series yt = C * itsMap[
Y] - S * itsMap[
X];
512 itsMap[
X] = C * xt + (S / k) * pxt;
513 itsMap[
Y] = C * yt + (S / k) * pyt;
514 itsMap[
PX] = C * pxt - (S * k) * xt;
515 itsMap[
PY] = C * pyt - (S * k) * yt;
516 itsMap[
T] += pt * (refTime / E - length / pz);
544 static const Map ident;
545 Map translate = co + ident;
548 Series E_trans =
sqrt((pt * pt + kin * kin).substitute(translate), order)
550 Series pz_trans =
sqrt((pt * pt - px * px - py * py).substitute(translate), order);
552 Series H_trans = E_trans - pz_trans;
557 for(
int i = 0; i <
PSdim; i += 2) {
562 static const Vector zeroV;
563 Vector new_co = co + implicitInt4(zeroV, JgradH, length, 0.5 * length);
578 double ca =
cos(angle);
579 double sa =
sin(angle);
580 double ta =
tan(angle);
582 int order = field.
order();
585 for(
int i = order; --i >= 1;) by = by * x + field.
normal(i);
595 itsMap[
T] -= ellovpp * p;
600 ps =
sqrt(p * p - itsMap[PX] * itsMap[PX] - itsMap[PY] * itsMap[PY], order);
601 itsMap[
X] = itsMap[
X] / (ca * (1.0 + ta * itsMap[
PX] / ps));
602 itsMap[
PX] = ca * itsMap[
PX] - sa * ps;
603 ellovpp = sa * itsMap[
X] / ps;
604 itsMap[
Y] -= ellovpp * itsMap[
PY];
605 itsMap[
T] += ellovpp * p;
614 double ca =
cos(angle);
615 double sa =
sin(angle);
616 double ta =
tan(angle);
618 int order = field.
order();
621 for(
int i = order; --i >= 1;) by = by * x + field.
normal(i);
629 itsMap[
X] = itsMap[
X] / (ca * (1.0 + ta * itsMap[
PX] / ps));
630 itsMap[
PX] = ca * itsMap[
PX] - sa * ps;
631 Series ellovpp = sa * itsMap[
X] / ps;
632 itsMap[
Y] -= ellovpp * itsMap[
PY];
633 itsMap[
T] += ellovpp * p;
638 ps =
sqrt(p * p - itsMap[PX] * itsMap[PX] - itsMap[PY] * itsMap[PY], order);
639 itsMap[
X] = itsMap[
X] / (ca * (1.0 - ta * itsMap[
PX] / ps));
640 itsMap[
PX] = ca * itsMap[
PX] + sa * ps;
641 ellovpp = sa * itsMap[
X] / ps;
642 itsMap[
Y] += ellovpp * itsMap[
PY];
643 itsMap[
T] -= ellovpp * p;
647 Vector implicitInt4(
const Vector &zin,
const VSeries &f,
double s,
double ds,
int nx,
int cx) {
659 static const double yt =
pow(2., 1 / 3.);
660 static const double ya = 1 / (2. - yt);
661 static const double yb = -yt * ya;
665 for(
int i = 0; i <
PSdim; ++i)
666 for(
int j = 0; j <
PSdim; ++j)
667 gradf[i][j] = f[i].derivative(j);
671 if(s < 0.) dsc = -dsc;
680 if(
std::abs(st + dsc) > as) dsc = s - st;
681 zt = ::implicitIntStep(zf, f, gradf, ya * dsc, nx);
682 zt = ::implicitIntStep(zt, f, gradf, yb * dsc, nx);
683 zt = ::implicitIntStep(zt, f, gradf, ya * dsc, nx);
686 std::string msg =
"Convergence not achieved within " + NumToStr<int>(cx) +
" cuts of step-size!";
693 if(ok) {zf = zt; st += dsc;}
717 static const double thresh = 1.e-8;
720 double ds2 = 0.5 * ds;
728 while(bcount < PSdim) {
730 std::string msg =
"Convergence not achieved within " + NumToStr<int>(nx) +
" iterations!";
735 Vector zt = 0.5 * (zin + z);
736 Matrix gf, idgf, idgf_inv;
737 for(
int i = 0; i <
PSdim; ++i)
738 for(
int j = 0; j <
PSdim; ++j)
739 gf[i][j] = -ds2 * gradf[i][j].evaluate(zt);
741 for(
int i = 0; i <
PSdim; ++i) idgf[i][i] += 1.;
743 idgf_inv = lu.inverse();
746 zf = idgf_inv * (zin + ds * f.
constantTerm(zt) + gf * z);
756 for(
int i = 0; i <
PSdim; ++i) {
758 {bounce[i] =
true; ++bcount;}
virtual void visitDrift(const Drift &)
Apply the algorithm to a Drift.
virtual double getExitFaceRotation() const =0
Get exit pole face rotation.
virtual void visitBeamBeam(const BeamBeam &)
Apply the algorithm to a BeamBeam.
virtual BMultipoleField & getField() override=0
Get multipole expansion of field.
virtual void visitLambertson(const Lambertson &)
Apply the algorithm to a Lambertson.
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
void applyThinSBend(const BMultipoleField &field, double scale, double h)
Thin SBend kick.
void setMinOrder(int order)
Set minimum order.
virtual BDipoleField & getField()=0
Return the corrector field.
virtual void visitMarker(const Marker &)
Apply the algorithm to a Marker.
FVps< double, 6 > VSeries
virtual double getAmplitude() const =0
Get RF amplitude.
double normal(int) const
Get component.
virtual void visitRBend(const RBend &)
Apply the algorithm to a RBend.
The field of a magnetic dipole.
virtual BMultipoleField & getField() override=0
Get multipole expansion of field.
A templated representation for matrices.
Interface for septum magnet.
virtual BMultipoleField & getField() override=0
Get multipole field.
virtual void visitSeparator(const Separator &)
Apply the algorithm to a Separator.
Interface for electrostatic separator.
Euclid3D getExitPatch() const
Get patch.
virtual double getSlices() const =0
Get number of slices.
virtual double getPhase() const =0
Get RF phase.
A templated representation for vectors.
Interface for beam position monitors.
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Interface for RF Quadrupole.
A simple arc in the XZ plane.
Tps< T > sin(const Tps< T > &x)
Sine.
virtual double getBz() const =0
Get solenoid field Bz in Teslas.
PETE_TUTree< FnCeil, typename T::PETE_Expr_t > ceil(const PETE_Expr< T > &l)
virtual void visitSBend(const SBend &)
Apply the algorithm to a SBend.
virtual double getElementLength() const override
Get design length.
virtual double getEx() const =0
Get horizontal component Ex of field in V/m.
Interface for general corrector.
virtual double getEy() const =0
Get vertical component Ey of field in V/m.
virtual void visitParallelPlate(const ParallelPlate &)
Apply the algorithm to a ParallelPlate.
Interface for beam diagnostics.
Euclid3D getEntrancePatch() const
Get patch.
virtual double getBy() const
Get vertical component.
virtual PlanarArcGeometry & getGeometry() override=0
Get SBend geometry.
virtual void visitMonitor(const Monitor &)
Apply the algorithm to a Monitor.
Tps< T > tan(const Tps< T > &x)
Tangent.
virtual double getEntryFaceRotation() const =0
Get pole entry face rotation.
Interface for drift space.
virtual void visitDiagnostic(const Diagnostic &)
Apply the algorithm to a Diagnostic.
virtual void visitBeamStripping(const BeamStripping &)
Apply the algorithm to a BeamStripping.
double getBeta() const
The relativistic beta per particle.
virtual double getStepsize() const =0
Get stepsize.
Interface for general multipole.
FVector< double, 6 > Vector
static Euclid3D YRotation(double angle)
Make rotation.
T deg(T x)
Convert radians to degrees.
Euclid3D Inverse(const Euclid3D &t)
Euclidean inverse.
virtual double getElementLength() const override
Get design length.
virtual double getStepsize() const =0
Get stepsize.
virtual double getExitFaceRotation() const =0
Get exit pole face rotation.
double getBendAngle() const
Get angle.
virtual void visitFlexibleCollimator(const FlexibleCollimator &)
Apply the algorithm to a flexible collimator.
A templated representation of a LU-decomposition.
FVector< T, N > constantTerm() const
Extract the constant part of the map.
double getQ() const
The constant charge per particle.
virtual double getElementLength() const
Get design length.
virtual double getSlices() const =0
Get number of slices.
virtual double getFrequency() const =0
Get RF frequencey.
constexpr double c
The velocity of light in m/s.
Interface for cyclotron collimator.
virtual void visitCorrector(const Corrector &)
Apply the algorithm to a Corrector.
Displacement and rotation in space.
Convergence error exception.
Abstract beam-beam interaction.
FMatrix< FTps< double, 6 >, 6, 6 > MxSeries
FTps< double, 6 > buildSBendVectorPotential(const BMultipoleField &, double h)
Construct the vector potential for an SBend.
virtual void visitRFQuadrupole(const RFQuadrupole &)
Apply the algorithm to a RFQuadrupole.
Interface for cyclotron valley.
virtual void visitDegrader(const Degrader &)
Apply the algorithm to a Degrader.
void applyThinMultipole(const BMultipoleField &field, double factor)
Thin multipole kick.
virtual double getBx() const
Get horizontal component.
FVps< T, N > ExpMap(const FTps< T, N > &H, int trunc=FTps< T, N >::EXACT)
Build the exponential series.
virtual void visitSolenoid(const Solenoid &)
Apply the algorithm to a Solenoid.
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
virtual double getElementLength() const
Get element length.
An abstract sequence of beam line components.
virtual double getEntryFaceCurvature() const =0
Get entry pole face curvature.
FTps substitute(const FMatrix< T, N, N > &M, int n) const
Substitute.
virtual double getElementLength() const
Get element length.
void applyEntranceFringe(double edge, double curve, const BMultipoleField &field, double scale)
Tps< T > sqrt(const Tps< T > &x)
Square root.
virtual double getEntryFaceRotation() const =0
Get pole entry face rotation.
FTps derivative(int var) const
Partial derivative.
The geometry for a RBend element.
double getP() const
The constant reference momentum per particle.
void applyTransform(const Euclid3D &, double refLength=0.0)
Apply transform.
The magnetic field of a multipole.
FTps< double, 6 > buildMultipoleVectorPotential(const BMultipoleField &)
Construct the vector potential for a Multipole.
double getM() const
The constant mass per particle.
virtual void visitSeptum(const Septum &)
Apply the algorithm to a Septum.
Tps< T > cos(const Tps< T > &x)
Cosine.
virtual double getExitFaceCurvature() const =0
Get exit pole face curvature.
void setMinOrder(int order)
Set minimum order.
FMatrix< double, 6, 6 > Matrix
virtual void visitProbe(const Probe &)
Apply the algorithm to a Probe.
double getCurvature() const
Get curvature.
virtual void visitCCollimator(const CCollimator &)
Apply the algorithm to a collimator.
virtual void visitMultipole(const Multipole &)
Apply the algorithm to a Multipole.
static const int EXACT
Representation of infinite precision.
virtual RBendGeometry & getGeometry() override=0
Get RBend geometry.
virtual double getBendAngle() const
Get angle.
void applyExitFringe(double edge, double curve, const BMultipoleField &field, double scale)
void setTruncOrder(int order)
Set truncation order.
virtual void visitCyclotronValley(const CyclotronValley &)
Apply the algorithm to a CyclotronValley.
Truncated power series in N variables of type T.
void applyDrift(double length)
const PartData itsReference
The reference information.
virtual void visitRFCavity(const RFCavity &)
Apply the algorithm to a RFCavity.
Vector truncated power series in n variables.
static FTps makeVariable(int var)
Make variable.
int order() const
Return order.
FVps< double, 6 > itsMap
The transfer map being built.
Interface for a Lambertson septum.