26 #include "gsl/gsl_interp.h"
27 #include "gsl/gsl_spline.h"
44 scaleCore_m(right.scaleCore_m),
45 scaleCoreError_m(right.scaleCoreError_m),
46 phaseCore1_m(right.phaseCore1_m),
47 phaseCore2_m(right.phaseCore2_m),
48 phaseExit_m(right.phaseExit_m),
49 length_m(right.length_m),
50 startCoreField_m(right.startCoreField_m),
51 startExitField_m(right.startExitField_m),
52 mappedStartExitField_m(right.mappedStartExitField_m),
53 PeriodLength_m(right.PeriodLength_m),
54 NumCells_m(right.NumCells_m),
55 CellLength_m(right.CellLength_m),
58 autophaseVeto_m(right.autophaseVeto_m),
59 designEnergy_m(right.designEnergy_m)
67 scaleCoreError_m(0.0),
72 startCoreField_m(0.0),
73 startExitField_m(0.0),
74 mappedStartExitField_m(0.0),
80 autophaseVeto_m(false),
102 Inform msg(
"TravelingWave::addK()");
104 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
105 Vector_t tmpE_diff(0.0, 0.0, 0.0), tmpB_diff(0.0, 0.0, 0.0);
125 const double z = tmpA0(2);
177 double oldK = tempK(0);
183 K +=
Vector_t(oldK * dx, oldK * dy, 0.0);
226 double tmpcos, tmpsin;
227 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
236 Vector_t tmpE2(0.0, 0.0, 0.0), tmpB2(0.0, 0.0, 0.0);
238 const double z = tmpR(2);
279 double tmpcos, tmpsin;
280 Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
288 Vector_t tmpE2(0.0, 0.0, 0.0), tmpB2(0.0, 0.0, 0.0);
290 const double z = tmpR(2);
337 Inform msg(
"TravelingWave ", *gmsg);
338 std::stringstream errormsg;
344 double zBegin = 0.0, zEnd = 0.0, rBegin = 0.0, rEnd = 0.0;
351 errormsg <<
"FREQUENCY IN INPUT FILE DIFFERENT THAN IN FIELD MAP '" <<
filename_m +
"';\n"
355 *gmsg << errormsg_str <<
endl;
356 ERRORMSG(errormsg_str <<
"\n" << endl);
358 std::ofstream omsg(
"errormsg.txt", std::ios_base::app);
383 endField = startField - 1
e-3;
386 endField = startField - 1
e-3;
428 std::vector<double> t, E, t2, E2;
429 std::vector<std::pair<double, double> > F;
433 double phi = 0.0, tmp_phi, dphi = 0.5 *
Physics::pi / 180.;
439 if(F.size() == 0)
return 0.0;
441 N1 =
static_cast<int>(
floor(F.size() / 4.)) + 1;
442 N2 = F.size() - 2 * N1 + 1;
444 N4 =
static_cast<int>(
floor(0.5 + N2 * Mode_m));
445 Dz = F[N1 + N2].first - F[N1].first;
451 for(
int i = 1; i < N1; ++ i) {
452 E[i] = E0 + (F[i].first + F[i - 1].first) / 2. *
scale_m / mass;
455 for(
int i = N1; i < N3 - N1 + 1; ++ i) {
456 int I = (i - N1) % N2 + N1;
457 double Z = (F[I].first + F[I - 1].first) / 2. +
floor((i - N1) / N2) * Dz;
461 for(
int i = N3 - N1 + 1; i < N3; ++ i) {
462 int I = i - N3 - 1 + 2 * N1 + N2;
463 double Z = (F[I].first + F[I - 1].first) / 2. + ((
NumCells_m - 1) * Mode_m - 1) * Dz;
464 E[i] = E0 + Z *
scale_m / mass;
468 for(
int iter = 0; iter < 10; ++ iter) {
470 for(
int i = 1; i < N1; ++ i) {
471 t[i] = t[i - 1] +
getdT(i, i, E, F, mass);
472 t2[i] = t2[i - 1] +
getdT(i, i, E2, F, mass);
476 for(
int i = N1; i < N3 - N1 + 1; ++ i) {
477 int I = (i - N1) % N2 + N1;
478 int J = (i - N1 + N4) % N2 + N1;
479 t[i] = t[i - 1] +
getdT(i, I, E, F, mass);
480 t2[i] = t2[i - 1] +
getdT(i, I, E2, F, mass);
484 for(
int i = N3 - N1 + 1; i < N3; ++ i) {
485 int I = i - N3 - 1 + 2 * N1 + N2;
486 t[i] = t[i - 1] +
getdT(i, I, E, F, mass);
487 t2[i] = t2[i - 1] +
getdT(i, I, E2, F, mass);
493 tmp_phi =
atan(A / B);
497 if(q * (A *
sin(tmp_phi) + B *
cos(tmp_phi)) < 0) {
502 for(
int i = 1; i < N1; ++ i) {
505 for(
int i = N1; i < N3 - N1 + 1; ++ i) {
506 int I = (i - N1) % N2 + N1;
507 int J = (i - N1 + N4) % N2 + N1;
508 E[i] = E[i - 1] + q *
scaleCore_m * (
getdE(i, I, t, phi + phaseC1, F) +
getdE(i, J, t, phi + phaseC2, F));
510 for(
int i = N3 - N1 + 1; i < N3; ++ i) {
511 int I = i - N3 - 1 + 2 * N1 + N2;
512 E[i] = E[i - 1] + q *
scale_m *
getdE(i, I, t, phi + phaseE, F);
516 INFOMSG(
level2 <<
"estimated phase= " << tmp_phi <<
" rad = "
518 <<
"Ekin= " << E[N3 - 1] <<
" MeV" << std::setprecision(prevPrecision) <<
"\n" <<
endl);
524 for(
int i = 1; i < N1; ++ i) {
526 E2[i] = E2[i - 1] + q *
scale_m *
getdE(i, i, t, phi + dphi, F);
527 t[i] = t[i - 1] +
getdT(i, i, E, F, mass);
528 t2[i] = t2[i - 1] +
getdT(i, i, E2, F, mass);
530 E2[i] = E2[i - 1] + q *
scale_m *
getdE(i, i, t2, phi + dphi, F);
532 for(
int i = N1; i < N3 - N1 + 1; ++ i) {
533 int I = (i - N1) % N2 + N1;
534 int J = (i - N1 + N4) % N2 + N1;
535 E[i] = E[i - 1] + q *
scaleCore_m * (
getdE(i, I, t, phi + phaseC1, F) +
getdE(i, J, t, phi + phaseC2, F));
536 E2[i] = E2[i - 1] + q *
scaleCore_m * (
getdE(i, I, t, phi + phaseC1 + dphi, F) +
getdE(i, J, t, phi + phaseC2 + dphi, F));
537 t[i] = t[i - 1] +
getdT(i, I, E, F, mass);
538 t2[i] = t2[i - 1] +
getdT(i, I, E2, F, mass);
539 E[i] = E[i - 1] + q *
scaleCore_m * (
getdE(i, I, t, phi + phaseC1, F) +
getdE(i, J, t, phi + phaseC2, F));
540 E2[i] = E2[i - 1] + q *
scaleCore_m * (
getdE(i, I, t2, phi + phaseC1 + dphi, F) +
getdE(i, J, t2, phi + phaseC2 + dphi, F));
542 for(
int i = N3 - N1 + 1; i < N3; ++ i) {
543 int I = i - N3 - 1 + 2 * N1 + N2;
544 E[i] = E[i - 1] + q *
scale_m *
getdE(i, I, t, phi + phaseE, F);
545 E2[i] = E2[i - 1] + q *
scale_m *
getdE(i, I, t, phi + phaseE + dphi, F);
546 t[i] = t[i - 1] +
getdT(i, I, E, F, mass);
547 t2[i] = t2[i - 1] +
getdT(i, I, E2, F, mass);
548 E[i] = E[i - 1] + q *
scale_m *
getdE(i, I, t, phi + phaseE, F);
549 E2[i] = E2[i - 1] + q *
scale_m *
getdE(i, I, t2, phi + phaseE + dphi, F);
556 INFOMSG(
level2 <<
"estimated phase= " << tmp_phi <<
" rad = "
558 <<
"Ekin= " << E[N3 - 1] <<
" MeV" << std::setprecision(prevPrecision) <<
"\n" <<
endl);
568 std::ofstream *out) {
574 std::vector<std::pair<double, double> > F;
577 double *zvals =
new double[F.size()];
578 double *onAxisField =
new double[F.size()];
580 for(
unsigned int i = 0; i < F.size(); ++i) {
581 zvals[i] = F[i].first;
582 onAxisField[i] = F[i].second;
589 gsl_spline *onAxisInterpolants = gsl_spline_alloc(gsl_interp_cspline, F.size());
590 gsl_interp_accel *onAxisAccel = gsl_interp_accel_alloc();
591 gsl_spline_init(onAxisInterpolants, zvals, onAxisField, F.size());
594 delete[] onAxisField;
596 if (out) *out << std::setw(18) << z
599 double dz = 0.5 * p /
sqrt(1.0 + p * p) * cdt;
602 double ez =
scale_m / Ezmax *
cos(phase) * gsl_spline_eval(onAxisInterpolants, z, onAxisAccel);
603 p += ez * q * cdt / mass;
604 dz = 0.5 * p /
sqrt(1.0 + p * p) * cdt;
605 z += 0.5 * p /
sqrt(1.0 + p * p) * cdt;
618 double ez =
scaleCore_m / Ezmax *
cos(phase) * gsl_spline_eval(onAxisInterpolants, tmpz, onAxisAccel);
624 ez +=
scaleCore_m / Ezmax *
cos(phase2) * gsl_spline_eval(onAxisInterpolants, tmpz, onAxisAccel);
626 p += ez * q * cdt / mass;
627 dz = 0.5 * p /
sqrt(1.0 + p * p) * cdt;
633 if (out) *out << std::setw(18) << z
643 double ez =
scale_m / Ezmax *
cos(phase) * gsl_spline_eval(onAxisInterpolants, tmpz, onAxisAccel);
644 p += ez * q * cdt / mass;
645 dz = 0.5 * p /
sqrt(1.0 + p * p) * cdt;
651 gsl_spline_free(onAxisInterpolants);
652 gsl_interp_accel_free(onAxisAccel);
654 const double beta =
sqrt(1. - 1 / (p * p + 1.));
657 return std::pair<double, double>(p, t - tErr);
ParticleAttrib< Vector_t > P
virtual void getInfo(Inform *msg)=0
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
virtual void accept(BeamlineVisitor &) const override
Apply visitor to TravelingWave.
constexpr double e
The value of .
double getdE(const int &i, const int &I, const std::vector< double > &t, const double &phi, const std::vector< std::pair< double, double > > &F) const
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 .
virtual bool isInside(const Vector_t &r) const
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) override
virtual void addKR(int i, double t, Vector_t &K) override
double getdT(const int &i, const int &I, const std::vector< double > &E, const std::vector< std::pair< double, double > > &F, const double mass) const
virtual void getFieldDimensions(double &zBegin, double &zEnd, double &rBegin, double &rEnd) const =0
virtual const std::string & getName() const
Get element name.
static std::string typeset_msg(const std::string &msg, const std::string &title)
static Fieldmap * getFieldmap(std::string Filename, bool fast=false)
virtual double getY(int i)
constexpr double rad2deg
The conversion factor from radians to degrees.
virtual void finalise() override
double getdB(const int &i, const int &I, const std::vector< double > &t, const double &phi, const std::vector< std::pair< double, double > > &F) const
virtual void getElementDimensions(double &begin, double &end) const override
Inform & level2(Inform &inf)
virtual void goOnline(const double &kineticEnergy) override
virtual ElementBase::ElementType getType() const override
Get element type std::string.
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const =0
virtual double getGamma(int i)
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
constexpr double pi
The value of .
virtual void getDimensions(double &zBegin, double &zEnd) const override
virtual bool isInside(const Vector_t &r) const override
Fieldmap * CoreFieldmap_m
virtual double getElementLength() const override
Get design length.
virtual double getBeta(int i)
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
Vektor< double, 3 > Vector_t
virtual double getFrequency() const =0
Tps< T > sqrt(const Tps< T > &x)
Square root.
virtual double getAutoPhaseEstimate(const double &E0, const double &t0, const double &q, const double &m) override
virtual void addKT(int i, double t, Vector_t &K) override
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) override
virtual double getX0(int i)
Tps< T > cos(const Tps< T > &x)
Cosine.
constexpr double q_e
The elementary charge in As.
double getEnergy(Vector_t p, double mass)
virtual double getZ(int i)
bool isInsideTransverse(const Vector_t &r, double f=1) const
virtual void visitTravelingWave(const TravelingWave &)=0
Apply the algorithm to a RF cavity.
virtual double getX(int i)
double mappedStartExitField_m
virtual bool bends() const override
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)
double getdA(const int &i, const int &I, const std::vector< double > &t, const double &phi, const std::vector< std::pair< double, double > > &F) const
virtual bool applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B) override