diff --git a/src/Algorithms/CavityAutophaser.cpp b/src/Algorithms/CavityAutophaser.cpp
index 2fa5ade912b68fa87dc666a13d99370eac1af316..307c9c8bb2cb6cb52d0c6e4f5b5823930a292a52 100644
--- a/src/Algorithms/CavityAutophaser.cpp
+++ b/src/Algorithms/CavityAutophaser.cpp
@@ -35,21 +35,28 @@ double CavityAutophaser::getPhaseAtMaxEnergy(const Vector_t &R,
}
initialP_m = Vector_t(0, 0, euclidean_norm(P));
- double tErr = (initialR_m(2) - R(2)) * sqrt(dot(P,P) + 1.0) / (P(2) * Physics::c);
- double initialEnergy = Util::getEnergy(P, itsReference_m.getM()) * 1e-6;
- double newPhase = 0.0, AstraPhase = 0.0;
- double initialPhase = guessCavityPhase(t + tErr);
- double optimizedPhase = 0.0;
- double finalEnergy = 0.0;
- RFCavity *element = static_cast<RFCavity *>(itsCavity_m.get());
- double amplitude = element->getAmplitudem();
- double designEnergy = element->getDesignEnergy();
- double originalPhase = element->getPhasem();
- bool apVeto = element->getAutophaseVeto();
- double basePhase = std::fmod(element->getFrequencym() * (t + tErr), Physics::two_pi);
+ RFCavity *element = static_cast<RFCavity *>(itsCavity_m.get());
+ bool apVeto = element->getAutophaseVeto();
+ double originalPhase = element->getPhasem();
+ double tErr = (initialR_m(2) - R(2)) * sqrt(dot(P,P) + 1.0) / (P(2) * Physics::c);
+ double optimizedPhase = 0.0;
+ double finalEnergy = 0.0;
+ double newPhase = 0.0;
+ double amplitude = element->getAmplitudem();
+ double basePhase = std::fmod(element->getFrequencym() * (t + tErr), Physics::two_pi);
if (!apVeto) {
+ double initialEnergy = Util::getEnergy(P, itsReference_m.getM()) * 1e-6;
+ double AstraPhase = 0.0;
+ double initialPhase = guessCavityPhase(t + tErr);
+ double designEnergy = element->getDesignEnergy();
+
+ if (amplitude < 0.0) {
+ amplitude = -amplitude;
+ element->setAmplitudem(amplitude);
+ }
+
if (amplitude == 0.0 && designEnergy <= 0.0) {
throw OpalException("CavityAutophaser::getPhaseAtMaxEnergy()",
"neither amplitude or design energy given to cavity " + element->getName());
diff --git a/src/Classic/AbsBeamline/RFCavity.cpp b/src/Classic/AbsBeamline/RFCavity.cpp
index c2489e8791d060b2ba2a8131a1e7c67065e0f094..c3d262ddc20516831ecd8269245d4251d6ab65d4 100644
--- a/src/Classic/AbsBeamline/RFCavity.cpp
+++ b/src/Classic/AbsBeamline/RFCavity.cpp
@@ -604,6 +604,35 @@ ElementBase::ElementType RFCavity::getType() const {
return RFCAVITY;
}
+double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, double mass) {
+ const double dt = 1e-13;
+ const double dphi = pi / 18;
+ const double p0 = Util::getP(E0, mass);
+
+ double phi = 0.0;
+ setPhasem(phi);
+ std::pair<double, double> ret = trackOnAxisParticle(E0 / mass, t0, dt, q, mass);
+ double phimax = ret.first;
+ double Emax = ret.second;
+ phi += dphi;
+
+ for (unsigned int i = 1; i < 36; ++ i, phi += dphi) {
+ setPhasem(phi);
+ ret = trackOnAxisParticle(p0, t0, dt, q, mass);
+ if (ret.second > Emax) {
+ Emax = ret.second;
+ phimax = ret.second;
+ }
+ }
+
+ const int prevPrecision = Ippl::Info->precision(8);
+ INFOMSG(level2 << "estimated phase= " << phimax << " rad = "
+ << phimax * Physics::rad2deg << " deg \n"
+ << "Ekin= " << Emax << " MeV" << setprecision(prevPrecision) << endl);
+
+ return phimax;
+}
+
double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const double &q, const double &mass) {
vector<double> t, E, t2, E2;
std::vector< double > F;
@@ -695,8 +724,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
E2[i] += q * scale_m * getdE(i, t2, dz, phi + dphi, frequency_m, F);
double a = E[i], b = E2[i];
if (std::isnan(a) || std::isnan(b)) {
- throw GeneralClassicException("RFCavity::getAutoPhaseEstimate",
- "solution is nan");
+ return getAutoPhaseEstimateFallback(E0, t0, q, mass);
}
t[i] = t[i - 1] + getdT(i, E, dz, mass);
t2[i] = t2[i - 1] + getdT(i, E2, dz, mass);
diff --git a/src/Classic/AbsBeamline/RFCavity.h b/src/Classic/AbsBeamline/RFCavity.h
index 75c729f47c41eac085151f0e3ea9cc6fe6e52135..e4682c819c3a994ac910b5d3aa57bf8142daf7a3 100644
--- a/src/Classic/AbsBeamline/RFCavity.h
+++ b/src/Classic/AbsBeamline/RFCavity.h
@@ -97,12 +97,13 @@ public:
virtual bool getAutophaseVeto() const;
virtual double getAutoPhaseEstimate(const double & E0, const double & t0, const double & q, const double & m);
+ virtual double getAutoPhaseEstimateFallback(double E0, double t0, double q, double m);
virtual std::pair<double, double> trackOnAxisParticle(const double & p0,
- const double & t0,
- const double & dt,
- const double & q,
- const double & mass);
+ const double & t0,
+ const double & dt,
+ const double & q,
+ const double & mass);
virtual void addKR(int i, double t, Vector_t &K);
@@ -503,4 +504,4 @@ std::string RFCavity::getFrequencyModelName() {
return frequency_name_m;
}
-#endif // CLASSIC_RFCavity_HH
+#endif // CLASSIC_RFCavity_HH
\ No newline at end of file
diff --git a/src/Classic/Utilities/Util.h b/src/Classic/Utilities/Util.h
index 5ea8395ac627deec94823522b3c8859931c8e61a..58fcaeb3cf64fea52d72cdaecf329b6ee0dc59fd 100644
--- a/src/Classic/Utilities/Util.h
+++ b/src/Classic/Utilities/Util.h
@@ -26,6 +26,12 @@ namespace Util {
return (getGamma(p) - 1.0) * mass;
}
+ inline
+ double getP(double E, double mass) {
+ double gamma = E / mass + 1;
+ return sqrt(std::pow(gamma, 2.0) - 1.0);
+ }
+
inline
std::string getTimeString(double time, unsigned int precision = 3) {
std::string timeUnit(" [ps]");
@@ -127,12 +133,12 @@ namespace Util {
std::string chargeUnit(" [fC]");
charge *= 1e15;
-
+
if (std::abs(charge) > 1000.0) {
charge /= 1000.0;
chargeUnit = std::string(" [pC]");
}
-
+
if (std::abs(charge) > 1000.0) {
charge /= 1000.0;
chargeUnit = std::string(" [nC]");