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),
66 length_m(right.length_m),
70 angle_m(right.angle_m),
71 sinAngle_m(right.sinAngle_m),
72 cosAngle_m(right.cosAngle_m),
74 gapwidth_m(right.gapwidth_m),
79 num_points_m(right.num_points_m)
88 amplitude_td_m(nullptr),
89 frequency_td_m(nullptr),
97 autophaseVeto_m(false),
152 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
182 double kx = 0.0, ky = 0.0;
186 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
191 kx += -cf *
scale_m * tmpE(0) * cwtf;
192 ky += -cf *
scale_m * tmpE(1) * cwtf;
201 K +=
Vector_t(KR(1) * dx + kx, KR(1) * dy + ky, 0.0);
217 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
219 if (tmpR(2) >= 0.0 &&
222 if (outOfBounds)
return true;
238 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
240 if (tmpR(2) >= 0.0 &&
243 if (outOfBounds)
return true;
262 double rBegin = 0.0, rEnd = 0.0;
263 Inform msg(
"RFCavity ", *gmsg);
264 std::stringstream errormsg;
274 errormsg <<
"FREQUENCY IN INPUT FILE DIFFERENT THAN IN FIELD MAP '" <<
filename_m <<
"';\n"
280 std::ofstream omsg(
"errormsg.txt", std::ios_base::app);
289 endField = startField - 1
e-3;
295 std::shared_ptr<AbstractTimeDependence> freq_atd,
296 std::shared_ptr<AbstractTimeDependence> ampl_atd,
297 std::shared_ptr<AbstractTimeDependence> phase_atd) {
310 "failed to open file '" +
filename_m +
"', please check if it exists");
312 *gmsg <<
"* Read cavity voltage profile data" <<
endl;
323 "not enough data in file '" +
filename_m +
"', please check the data format");
325 in >>
RNormal_m[i] >> VrNormal_m[i] >> DvDr_m[i];
336 *gmsg <<
"* Cavity voltage data read successfully!" <<
endl;
417 if(name ==
"STANDING") {
419 }
else if(name ==
"SINGLEGAP") {
421 }
else if(name !=
"") {
422 std::stringstream errormsg;
423 errormsg <<
getName() <<
": CAVITY TYPE " << name <<
" DOES NOT EXIST;";
427 ofstream omsg(
"errormsg.txt", ios_base::app);
428 omsg << errormsg_str <<
endl;
440 return std::string(
"SINGLEGAP");
442 return std::string(
"STANDING");
459 void RFCavity::getMomentaKick(
const double normalRadius,
double momentum[],
const double t,
const double dtCorrt,
const int PID,
const double restMass,
const int chargenumber) {
466 double momentum2 = momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2];
467 double betgam =
sqrt(momentum2);
469 double gamma =
sqrt(1.0 + momentum2);
470 double beta = betgam / gamma;
474 double transit_factor = 0.0;
475 double Ufactor = 1.0;
480 transit_factor = 0.5 * frequency *
gapwidth_m * 1.0e-3 / (
c * beta);
481 Ufactor =
sin(transit_factor) / transit_factor;
486 double nphase = (frequency * (t + dtCorrt) * 1.0
e-9) -
phi0_m / 180.0 *
pi ;
487 double dgam = Voltage *
cos(nphase) / (restMass);
489 double tempdegree =
fmod(nphase * 360.0 /
two_pi, 360.0);
490 if(tempdegree > 270.0) tempdegree -= 360.0;
494 double newmomentum2 =
pow(gamma, 2) - 1.0;
497 double ptheta =
sqrt(newmomentum2 -
pow(pr, 2));
499 double py = pr * sinAngle_m + ptheta *
cosAngle_m;
501 double rotate = -derivate * (
scale_m * 1.0e6) / ((
rmax_m -
rmin_m) / 1000.0) *
sin(nphase) / (frequency *
two_pi) / (betgam * restMass /
c / chargenumber);
504 momentum[0] =
cos(rotate) * px +
sin(rotate) * py;
505 momentum[1] = -
sin(rotate) * px +
cos(rotate) * py;
511 m <<
"* Cavity " <<
getName() <<
" Phase= " << tempdegree <<
" [deg] transit time factor= " << Ufactor
512 <<
" dE= " << dgam *restMass * 1.0e-6 <<
" [MeV]"
513 <<
" E_kin= " << (gamma - 1.0)*restMass * 1.0e-6 <<
" [MeV] Time dep freq = " <<
frequency_td_m->getValue(t) <<
endl;
525 "no support points!");
537 while((ih - il) > 1) {
538 int i = (int)((il + ih) / 2.0);
561 double u = (z - x1) / dx;
564 double dy2 = -2.0 * dy;
565 double ya2 = y2a + 2.0 * y1a;
566 double dy3 = 3.0 * dy;
567 double ya3 = y2a + y1a;
568 double yb2 = dy2 + dx * ya3;
569 double yb4 = dy3 - dx * ya2;
570 splint = y1 + u * dx * y1a + u2 * yb4 + u3 * yb2;
571 *za = y1a + 2.0 * u / dx * yb4 + 3.0 * u2 / dx * yb2;
591 const double dt = 1
e-13;
595 double dphi =
pi / 18;
604 for (
unsigned int j = 0; j < 2; ++ j) {
605 for (
unsigned int i = 0; i < 36; ++ i, phi += dphi) {
620 phimax =
fmod(phimax, Physics::two_pi);
624 <<
"estimated phase= " << phimax <<
" rad = "
626 <<
"Ekin= " << Emax <<
" MeV" << setprecision(prevPrecision) <<
"\n" <<
endl);
633 vector<double> t, E, t2, E2;
634 std::vector< double > F;
635 std::vector< std::pair< double, double > > G;
636 gsl_spline *onAxisInterpolants;
637 gsl_interp_accel *onAxisAccel;
641 double phi = 0.0, tmp_phi, dphi = 0.5 *
Physics::pi / 180.;
642 double dz = 1.0, length = 0.0;
644 double begin = (G.front()).first;
645 double end = (G.back()).first;
646 std::unique_ptr<double[]> zvals(
new double[G.size()]);
647 std::unique_ptr<double[]> onAxisField(
new double[G.size()]);
649 for(
size_t j = 0; j < G.size(); ++ j) {
650 zvals[j] = G[j].first;
651 onAxisField[j] = G[j].second;
653 onAxisInterpolants = gsl_spline_alloc(gsl_interp_cspline, G.size());
654 onAxisAccel = gsl_interp_accel_alloc();
655 gsl_spline_init(onAxisInterpolants, zvals.get(), onAxisField.get(), G.size());
657 length = end - begin;
658 dz = length / G.size();
662 N = (int)
floor(length / dz + 1);
667 for(
size_t j = 0; j < N; ++ j, z += dz) {
668 F[j] = gsl_spline_eval(onAxisInterpolants, z, onAxisAccel);
670 gsl_spline_free(onAxisInterpolants);
671 gsl_interp_accel_free(onAxisAccel);
679 for(
unsigned int i = 1; i < N; ++ i, z += dz) {
680 E[i] = E[i - 1] + dz *
scale_m / mass;
684 for(
int iter = 0; iter < 10; ++ iter) {
686 for(
unsigned int i = 1; i < N; ++ i) {
687 t[i] = t[i - 1] +
getdT(i, E, dz, mass);
688 t2[i] = t2[i - 1] +
getdT(i, E2, dz, mass);
694 tmp_phi =
atan(A / B);
698 if(q * (A *
sin(tmp_phi) + B *
cos(tmp_phi)) < 0) {
703 for(
unsigned int i = 1; i < N; ++ i) {
708 INFOMSG(
level2 <<
"estimated phase= " << tmp_phi <<
" rad = "
710 <<
"Ekin= " << E[N - 1] <<
" MeV" << setprecision(prevPrecision) <<
"\n" <<
endl);
716 for(
unsigned int i = 1; i < N; ++ i) {
721 double a = E[i], b = E2[i];
725 t[i] = t[i - 1] +
getdT(i, E, dz, mass);
726 t2[i] = t2[i - 1] +
getdT(i, E2, dz, mass);
734 double cosine_part = 0.0, sine_part = 0.0;
735 double p0 =
sqrt((E0 / mass + 1) * (E0 / mass + 1) - 1);
739 double totalEz0 =
cos(phi) * cosine_part -
sin(phi) * sine_part;
741 if(p0 + q * totalEz0 * (t[1] - t[0]) *
Physics::c / mass < 0) {
743 tmp_phi =
atan(cosine_part / sine_part);
754 <<
"estimated phase= " << tmp_phi <<
" rad = "
756 <<
"Ekin= " << E[N - 1] <<
" MeV" << setprecision(prevPrecision) <<
"\n" <<
endl);
766 std::ofstream *out) {
775 double dz = 0.5 * p(2) /
sqrt(1.0 +
dot(p, p)) * cdt;
778 if (out) *out << std::setw(18) << z[2]
781 while(z(2) + dz < zend && z(2) + dz > zbegin) {
783 integrator.push(z, p, dt);
786 if(z(2) >= zbegin && z(2) <= zend) {
791 integrator.kick(z, p, Ef, Bf, dt);
793 dz = 0.5 * p(2) /
sqrt(1.0 +
dot(p, p)) * cdt;
795 integrator.push(z, p, dt);
799 if (out) *out << std::setw(18) << z[2]
804 const double beta =
sqrt(1. - 1 / (
dot(p, p) + 1.));
805 const double tErr = (z(2) - zend) / (
Physics::c * beta);
807 return pair<double, double>(p(2), t - tErr);
821 double start, end, tmp;
823 length = end - start;
ParticleAttrib< Vector_t > P
virtual void getInfo(Inform *msg)=0
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
constexpr double e
The value of .
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
virtual void addKR(int i, double t, Vector_t &K) override
void setRmin(double rmin)
virtual ElementBase::ElementType getType() const override
Get element type std::string.
double getP(double E, double mass)
virtual void goOffline() override
virtual void getOnaxisEz(std::vector< std::pair< double, double > > &onaxis)
Tps< T > sin(const Tps< T > &x)
Sine.
constexpr double two_pi
The value of .
double spline(double z, double *za)
virtual bool isInside(const Vector_t &r) const
double getdB(const int &i, const std::vector< double > &t, const double &dz, const double &frequency, const std::vector< double > &F) const
virtual void setComponentType(std::string name) override
virtual double getElementLength() const override
Get design length.
std::string toUpper(const std::string &str)
void setRmax(double rmax)
virtual void getFieldDimensions(double &zBegin, double &zEnd, double &rBegin, double &rEnd) const =0
virtual void visitRFCavity(const RFCavity &)=0
Apply the algorithm to a RF cavity.
const PartData * getReference() const
virtual void getElementDimensions(double &begin, double &end) const override
virtual const std::string & getName() const
Get element name.
virtual void finalise() override
static std::string typeset_msg(const std::string &msg, const std::string &title)
std::string frequency_name_m
static Fieldmap * getFieldmap(std::string Filename, bool fast=false)
virtual double getY(int i)
void setAmplitudeModel(std::shared_ptr< AbstractTimeDependence > time_dep)
virtual double getRmin() const
constexpr double rad2deg
The conversion factor from radians to degrees.
double dot(const Vector3D &lhs, const Vector3D &rhs)
Vector dot product.
double getdT(const int &i, const std::vector< double > &E, const double &dz, const double mass) 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
std::unique_ptr< double[]> RNormal_m
std::unique_ptr< double[]> VrNormal_m
void setFrequencyModel(std::shared_ptr< AbstractTimeDependence > time_dep)
constexpr double m_e
The electron rest mass in GeV.
Inform & level2(Inform &inf)
virtual void getDimensions(double &zBegin, double &zEnd) const override
virtual double getAzimuth() const
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const =0
void setPhaseModel(std::shared_ptr< AbstractTimeDependence > time_dep)
virtual double getGamma(int i)
void setPhi0(double phi0)
virtual bool bends() const override
std::unique_ptr< double[]> DvDr_m
virtual double getPhi0() const
virtual double getGapWidth() const
void setAzimuth(double angle)
constexpr double pi
The value of .
void setElType(ElemType elt)
set the element type as enumeration needed in the envelope tracker
virtual double getPerpenDistance() const
virtual double getElementLength() const
Get design length.
virtual std::string getComponentType() const override
virtual double getBeta(int i)
virtual double getPhasem() const
PETE_TUTree< FnArcTan, typename T::PETE_Expr_t > atan(const PETE_Expr< T > &l)
constexpr double c
The velocity of light in m/s.
virtual double getY0(int i)
PartBunchBase< double, 3 > * RefPartBunch_m
void splint(double xa[], double ya[], double y2a[], int n, double x, double *y)
void setGapWidth(double gapwidth)
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Vektor< double, 3 > Vector_t
virtual double getFrequency() const =0
double getdA(const int &i, const std::vector< double > &t, const double &dz, const double &frequency, const std::vector< double > &F) const
PETE_TBTree< FnFmod, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > fmod(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
void setPerpenDistance(double pdis)
Tps< T > sqrt(const Tps< T > &x)
Square root.
virtual void addKT(int i, double t, Vector_t &K) override
virtual bool applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B) override
virtual void setPhasem(double phase)
virtual double getAutoPhaseEstimateFallback(double E0, double t0, double q, double m)
virtual double getCycFrequency() const
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 getRmax() const
virtual double getX0(int i)
virtual double getSinAzimuth() const
Tps< T > cos(const Tps< T > &x)
Cosine.
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) override
constexpr double q_e
The elementary charge in As.
double getQ() const
Access to reference data.
virtual void accept(BeamlineVisitor &) const override
Apply visitor to RFCavity.
T isnan(T x)
isnan function with adjusted return type
virtual bool isInside(const Vector_t &r) const override
double getEnergy(Vector_t p, double mass)
std::shared_ptr< AbstractTimeDependence > frequency_td_m
virtual double getZ(int i)
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 double getCosAzimuth() const
bool isInsideTransverse(const Vector_t &r, double f=1) const
Interface for a single beam element.
virtual double getX(int i)
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const =0
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Inform & endl(Inform &inf)
virtual double getAutoPhaseEstimate(const double &E0, const double &t0, const double &q, const double &m)
virtual void goOnline(const double &kineticEnergy) override