29 #include "gsl/gsl_interp.h"
30 #include "gsl/gsl_spline.h"
48 phase_td_m(right.phase_td_m),
49 phase_name_m(right.phase_name_m),
50 amplitude_td_m(right.amplitude_td_m),
51 amplitude_name_m(right.amplitude_name_m),
52 frequency_td_m(right.frequency_td_m),
53 frequency_name_m(right.frequency_name_m),
54 filename_m(right.filename_m),
55 scale_m(right.scale_m),
56 scaleError_m(right.scaleError_m),
57 phase_m(right.phase_m),
58 phaseError_m(right.phaseError_m),
59 frequency_m(right.frequency_m),
61 autophaseVeto_m(right.autophaseVeto_m),
62 designEnergy_m(right.designEnergy_m),
63 fieldmap_m(right.fieldmap_m),
64 startField_m(right.startField_m),
65 endField_m(right.endField_m),
69 angle_m(right.angle_m),
70 sinAngle_m(right.sinAngle_m),
71 cosAngle_m(right.cosAngle_m),
73 gapwidth_m(right.gapwidth_m),
78 num_points_m(right.num_points_m)
85 amplitude_td_m(nullptr),
86 frequency_td_m(nullptr),
94 autophaseVeto_m(false),
133 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
136 if (outOfBounds)
return true;
153 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
156 if (outOfBounds)
return true;
176 std::stringstream errormsg;
183 "The length of the field map '" +
filename_m +
"' is zero or negative");
189 errormsg <<
"FREQUENCY IN INPUT FILE DIFFERENT THAN IN FIELD MAP '" <<
filename_m <<
"';\n"
195 std::ofstream omsg(
"errormsg.txt", std::ios_base::app);
206 std::shared_ptr<AbstractTimeDependence> freq_atd,
207 std::shared_ptr<AbstractTimeDependence> ampl_atd,
208 std::shared_ptr<AbstractTimeDependence> phase_atd) {
220 "failed to open file '" +
filename_m +
"', please check if it exists");
222 *
gmsg <<
"* Read cavity voltage profile data" <<
endl;
233 "not enough data in file '" +
filename_m +
"', please check the data format");
246 *
gmsg <<
"* Cavity voltage data read successfully!" <<
endl;
326 if (
name ==
"STANDING") {
328 }
else if (
name ==
"SINGLEGAP") {
330 }
else if (
name !=
"") {
331 std::stringstream errormsg;
332 errormsg <<
getName() <<
": CAVITY TYPE " <<
name <<
" DOES NOT EXIST;";
336 std::ofstream omsg(
"errormsg.txt", std::ios_base::app);
349 return std::string(
"SINGLEGAP");
351 return std::string(
"STANDING");
368 void RFCavity::getMomentaKick(
const double normalRadius,
double momentum[],
const double t,
const double dtCorrt,
const int PID,
const double restMass,
const int chargenumber) {
372 double momentum2 = momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2];
375 double gamma =
std::sqrt(1.0 + momentum2);
376 double beta = betgam / gamma;
378 double Voltage =
spline(normalRadius, &derivate) *
scale_m * 1.0e6;
380 double Ufactor = 1.0;
386 Ufactor =
std::sin(transit_factor) / transit_factor;
391 double nphase = (frequency * (t + dtCorrt) * 1.0e-9) -
phi0_m / 180.0 *
Physics::pi ;
392 double dgam = Voltage *
std::cos(nphase) / (restMass);
395 if (tempdegree > 270.0) tempdegree -= 360.0;
399 double newmomentum2 =
std::pow(gamma, 2) - 1.0;
416 m <<
"* Cavity " <<
getName() <<
" Phase= " << tempdegree <<
" [deg] transit time factor= " << Ufactor
417 <<
" dE= " << dgam *restMass * 1.0e-6 <<
" [MeV]"
418 <<
" E_kin= " << (gamma - 1.0)*restMass * 1.0e-6 <<
" [MeV] Time dep freq = " <<
frequency_td_m->getValue(t) <<
endl;
430 "no support points!");
442 while ((ih - il) > 1) {
443 int i = (int)((il + ih) / 2.0);
466 double u = (z - x1) / dx;
469 double dy2 = -2.0 * dy;
470 double ya2 = y2a + 2.0 * y1a;
471 double dy3 = 3.0 * dy;
472 double ya3 = y2a + y1a;
473 double yb2 = dy2 + dx * ya3;
474 double yb4 = dy3 - dx * ya2;
475 splint = y1 + u * dx * y1a + u2 * yb4 + u3 * yb2;
476 *za = y1a + 2.0 * u / dx * yb4 + 3.0 * u2 / dx * yb2;
496 const double dt = 1
e-13;
508 for (
unsigned int j = 0; j < 2; ++ j) {
509 for (
unsigned int i = 0; i < 36; ++ i, phi += dphi) {
528 <<
"estimated phase= " << phimax <<
" rad = "
530 <<
"Ekin= " << Emax <<
" MeV" << std::setprecision(prevPrecision) <<
"\n" <<
endl);
537 std::vector<double> t, E, t2, E2;
538 std::vector<double> F;
539 std::vector< std::pair< double, double > > G;
540 gsl_spline *onAxisInterpolants;
541 gsl_interp_accel *onAxisAccel;
543 double phi = 0.0, tmp_phi, dphi = 0.5 *
Physics::pi / 180.;
544 double dz = 1.0, length = 0.0;
546 if (G.size() == 0)
return 0.0;
547 double begin = (G.front()).first;
548 double end = (G.back()).first;
549 std::unique_ptr<double[]> zvals(
new double[G.size()]);
550 std::unique_ptr<double[]> onAxisField(
new double[G.size()]);
552 for (
size_t j = 0; j < G.size(); ++ j) {
553 zvals[j] = G[j].first;
554 onAxisField[j] = G[j].second;
556 onAxisInterpolants = gsl_spline_alloc(gsl_interp_cspline, G.size());
557 onAxisAccel = gsl_interp_accel_alloc();
558 gsl_spline_init(onAxisInterpolants, zvals.get(), onAxisField.get(), G.size());
561 dz = length / G.size();
565 unsigned int N = (int)
floor(length / dz + 1);
570 for (
size_t j = 0; j < N; ++ j, z += dz) {
571 F[j] = gsl_spline_eval(onAxisInterpolants, z, onAxisAccel);
573 gsl_spline_free(onAxisInterpolants);
574 gsl_interp_accel_free(onAxisAccel);
582 for (
unsigned int i = 1; i < N; ++ i, z += dz) {
583 E[i] = E[i - 1] + dz *
scale_m / mass;
587 for (
int iter = 0; iter < 10; ++ iter) {
590 for (
unsigned int i = 1; i < N; ++ i) {
591 t[i] = t[i - 1] +
getdT(i, E, dz, mass);
592 t2[i] = t2[i - 1] +
getdT(i, E2, dz, mass);
598 tmp_phi =
atan(A / B);
602 if (q * (A *
sin(tmp_phi) + B *
cos(tmp_phi)) < 0) {
607 for (
unsigned int i = 1; i < N; ++ i) {
612 INFOMSG(
level2 <<
"estimated phase= " << tmp_phi <<
" rad = "
614 <<
"Ekin= " << E[N - 1] <<
" MeV" << std::setprecision(prevPrecision) <<
"\n" <<
endl);
620 for (
unsigned int i = 1; i < N; ++ i) {
625 double a = E[i], b = E2[i];
629 t[i] = t[i - 1] +
getdT(i, E, dz, mass);
630 t2[i] = t2[i - 1] +
getdT(i, E2, dz, mass);
638 double cosine_part = 0.0, sine_part = 0.0;
643 double totalEz0 =
std::cos(phi) * cosine_part -
sin(phi) * sine_part;
645 if (p0 + q * totalEz0 * (t[1] - t[0]) *
Physics::c / mass < 0) {
647 tmp_phi =
std::atan(cosine_part / sine_part);
658 <<
"estimated phase= " << tmp_phi <<
" rad = "
660 <<
"Ekin= " << E[N - 1] <<
" MeV" << std::setprecision(prevPrecision) <<
"\n" <<
endl);
670 std::ofstream *out) {
682 if (out) *out << std::setw(18) << z[2]
685 while (z(2) + dz < zend && z(2) + dz > zbegin) {
687 integrator.
push(z, p, dt);
692 if (z(2) >= zbegin && z(2) <= zend) {
695 integrator.
kick(z, p, Ef, Bf, dt);
699 integrator.
push(z, p, dt);
703 if (out) *out << std::setw(18) << z[2]
708 const double beta =
std::sqrt(1. - 1 / (
dot(p, p) + 1.));
709 const double tErr = (z(2) - zend) / (
Physics::c * beta);
711 return std::pair<double, double>(p(2), t - tErr);
727 length =
end - start;
Tps< T > cos(const Tps< T > &x)
Cosine.
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Tps< T > sin(const Tps< T > &x)
Sine.
Tps< T > sqrt(const Tps< T > &x)
Square root.
PartBunchBase< T, Dim >::ConstIterator end(PartBunchBase< T, Dim > const &bunch)
PartBunchBase< T, Dim >::ConstIterator begin(PartBunchBase< T, Dim > const &bunch)
double dot(const Vector3D &lhs, const Vector3D &rhs)
Vector dot product.
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
PETE_TBTree< FnFmod, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > fmod(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
PETE_TUTree< FnArcTan, typename T::PETE_Expr_t > atan(const PETE_Expr< T > &l)
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Inform & endl(Inform &inf)
Inform & level2(Inform &inf)
constexpr double two_pi
The value of.
constexpr double deg2rad
The conversion factor from degrees to radians.
constexpr double e
The value of.
constexpr double c
The velocity of light in m/s.
constexpr double pi
The value of.
constexpr double rad2deg
The conversion factor from radians to degrees.
T isnan(T x)
isnan function with adjusted return type
double getKineticEnergy(Vector_t p, double mass)
double getBetaGamma(double Ekin, double mass)
double getGamma(Vector_t p)
const PartData * getReference() const
double getQ() const
Access to reference data.
ParticleAttrib< Vector_t > P
virtual void visitRFCavity(const RFCavity &)=0
Apply the algorithm to a RF cavity.
Interface for a single beam element.
PartBunchBase< double, 3 > * RefPartBunch_m
virtual const std::string & getName() const
Get element name.
virtual double getElementLength() const
Get design length.
virtual void setElementLength(double length)
Set design length.
bool isInsideTransverse(const Vector_t &r) const
virtual std::string getComponentType() const override
void setFrequencyModel(std::shared_ptr< AbstractTimeDependence > time_dep)
virtual bool isInside(const Vector_t &r) const override
virtual double getPhasem() const
virtual bool bends() const override
virtual double getRmax() const
void setPerpenDistance(double pdis)
void getMomentaKick(const double normalRadius, double momentum[], const double t, const double dtCorrt, const int PID, const double restMass, const int chargenumber)
used in OPAL-cycl
virtual double getAzimuth() const
double getdE(const int &i, const std::vector< double > &t, const double &dz, const double &phi, const double &frequency, const std::vector< double > &F) const
virtual void accept(BeamlineVisitor &) const override
Apply visitor to RFCavity.
double getdB(const int &i, const std::vector< double > &t, const double &dz, const double &frequency, const std::vector< double > &F) const
std::unique_ptr< double[]> DvDr_m
virtual void finalise() override
void setRmin(double rmin)
void setPhaseModel(std::shared_ptr< AbstractTimeDependence > time_dep)
virtual std::pair< double, double > trackOnAxisParticle(const double &p0, const double &t0, const double &dt, const double &q, const double &mass, std::ofstream *out=NULL)
virtual void getElementDimensions(double &begin, double &end) const override
virtual double getCosAzimuth() const
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
std::shared_ptr< AbstractTimeDependence > frequency_td_m
std::unique_ptr< double[]> RNormal_m
virtual void setPhasem(double phase)
virtual void setComponentType(std::string name) override
virtual double getSinAzimuth() const
void setPhi0(double phi0)
virtual bool applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B) override
double spline(double z, double *za)
void setAzimuth(double angle)
std::unique_ptr< double[]> VrNormal_m
virtual void goOnline(const double &kineticEnergy) override
virtual double getCycFrequency() const
double getdA(const int &i, const std::vector< double > &t, const double &dz, const double &frequency, const std::vector< double > &F) const
virtual double getGapWidth() const
virtual double getAutoPhaseEstimate(const double &E0, const double &t0, const double &q, const double &m)
virtual void getDimensions(double &zBegin, double &zEnd) const override
virtual double getPerpenDistance() const
virtual ElementBase::ElementType getType() const override
Get element type std::string.
void setGapWidth(double gapwidth)
virtual double getElementLength() const override
Get design length.
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) override
virtual void goOffline() override
virtual double getPhi0() const
virtual double getAutoPhaseEstimateFallback(double E0, double t0, double q, double m)
std::string frequency_name_m
void setRmax(double rmax)
void setAmplitudeModel(std::shared_ptr< AbstractTimeDependence > time_dep)
virtual double getRmin() const
double getdT(const int &i, const std::vector< double > &E, const double &dz, const double mass) const
virtual void getInfo(Inform *msg)=0
virtual void getFieldDimensions(double &zBegin, double &zEnd) const =0
static Fieldmap * getFieldmap(std::string Filename, bool fast=false)
virtual bool isInside(const Vector_t &) const
virtual void getOnaxisEz(std::vector< std::pair< double, double > > &onaxis)
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const =0
virtual double getFrequency() const =0
static std::string typeset_msg(const std::string &msg, const std::string &title)
void kick(const Vector_t &R, Vector_t &P, const Vector_t &Ef, const Vector_t &Bf, const double &dt) const
void push(Vector_t &R, const Vector_t &P, const double &dt) const
Vektor< double, 3 > Vector_t