diff --git a/.gitignore b/.gitignore
index 6dbeb02b5764e4217656694881b4f01af02cd48f..1b64b6e89ba109e0ab247936a6be7ab51145a13c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,6 +3,7 @@
*.o
*.a
*.aux
+optimizer/Tests/*.exe
build
CMakeCache.txt
CMakeFiles
diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp
index 0e1943b378eda1ec35aeaf58c640f57dbbd808af..25244ea458835feccf5c338267e83e9a6f88d33c 100644
--- a/src/Algorithms/ParallelCyclotronTracker.cpp
+++ b/src/Algorithms/ParallelCyclotronTracker.cpp
@@ -53,6 +53,7 @@
#include "AbsBeamline/CyclotronValley.h"
#include "AbsBeamline/Stripper.h"
#include "AbsBeamline/VariableRFCavity.h"
+#include "AbsBeamline/VariableRFCavityFringeField.h"
#include "AbstractObjects/Element.h"
@@ -829,6 +830,16 @@ void ParallelCyclotronTracker::visitVariableRFCavity(const VariableRFCavity &cav
"Need to define a RINGDEFINITION to use VariableRFCavity element");
}
+void ParallelCyclotronTracker::visitVariableRFCavityFringeField
+ (const VariableRFCavityFringeField &cav) {
+ *gmsg << "Adding Variable RF Cavity with Fringe Field" << endl;
+ if (opalRing_m != NULL)
+ opalRing_m->appendElement(cav);
+ else
+ throw OpalException("ParallelCyclotronTracker::visitVariableRFCavityFringeField",
+ "Need to define a RINGDEFINITION to use VariableRFCavity element");
+}
+
/**
*
*
diff --git a/src/Algorithms/ParallelCyclotronTracker.h b/src/Algorithms/ParallelCyclotronTracker.h
index 3ecd7dc72a5788b95cd94416ee3f218b1d9199d3..4555dd546bd898df9ca2cf9e3e6644df2ed1e6c2 100644
--- a/src/Algorithms/ParallelCyclotronTracker.h
+++ b/src/Algorithms/ParallelCyclotronTracker.h
@@ -33,6 +33,7 @@ class PlanarArcGeometry;
class Ring;
class SBend3D;
class VariableRFCavity;
+class VariableRFCavityFringeField;
class Offset;
// Class ParallelCyclotronTracker
@@ -177,6 +178,10 @@ public:
/// Apply the algorithm to a VariabelRFCavity.
virtual void visitVariableRFCavity(const VariableRFCavity &cav);
+ /// Apply the algorithm to a VariabelRFCavity.
+ virtual void visitVariableRFCavityFringeField
+ (const VariableRFCavityFringeField &cav);
+
/// Apply the algorithm to the top-level beamline.
// overwrite the execute-methode from DefaultVisitor
virtual void execute();
diff --git a/src/BasicActions/DumpEMFields.cpp b/src/BasicActions/DumpEMFields.cpp
index cded6540ab1b3c4a437a5d9a6350f67c53284dfd..9ea7d346e060d6e7ddd677b8feb563d18ef4167f 100644
--- a/src/BasicActions/DumpEMFields.cpp
+++ b/src/BasicActions/DumpEMFields.cpp
@@ -41,7 +41,7 @@ std::unordered_set<DumpEMFields*> DumpEMFields::dumpsSet_m;
std::string DumpEMFields::dumpemfields_docstring =
std::string("The \"DUMPEMFIELDS\" statement dumps a field map to a user-defined")+
std::string(" field file, for checking that fields are generated correctly.")+
-std::string(" The fields are written out on a grid in Cartesian space and time.");
+std::string(" The fields are written out on a grid in space and time.");
DumpEMFields::DumpEMFields() :
Action(20, "DUMPEMFIELDS", dumpemfields_docstring.c_str()),
@@ -196,7 +196,8 @@ void DumpEMFields::writeFields(Component* field) {
}
void DumpEMFields::checkInt(double real, std::string name, double tolerance) {
- if (fabs(floor(real) - real) > tolerance) {
+ real += tolerance; // prevent rounding error
+ if (fabs(floor(real) - real) > 2*tolerance) {
throw OpalException("DumpEMFields::checkInt",
"Value for "+name+
" should be an integer but a real value was found");
@@ -210,30 +211,29 @@ void DumpEMFields::checkInt(double real, std::string name, double tolerance) {
void DumpEMFields::writeHeader(std::ofstream& fout) const {
fout << grid_m->end().toInteger() << "\n";
if (coordinates_m == CYLINDRICAL) {
- fout << 1 << " r [mm]\n";
- fout << 2 << " phi [degree]\n";
+ fout << 1 << " r [mm]\n";
+ fout << 2 << " phi [degree]\n";
} else {
- fout << 1 << " x [mm]\n";
- fout << 2 << " y [mm]\n";
+ fout << 1 << " x [mm]\n";
+ fout << 2 << " y [mm]\n";
}
- fout << 3 << " z [mm]\n";
- fout << 4 << " t [ns]\n";
+ fout << 3 << " z [mm]\n";
+ fout << 4 << " t [ns]\n";
if (coordinates_m == CYLINDRICAL) {
- fout << 5 << " Br [kGauss]\n";
- fout << 6 << " Bphi [kGauss]\n";
+ fout << 5 << " Br [kGauss]\n";
+ fout << 6 << " Bphi [kGauss]\n";
+ fout << 7 << " Bz [kGauss]\n";
+ fout << 8 << " Er [MV/m]\n";
+ fout << 9 << " Ephi [MV/m]\n";
+ fout << 10 << " Ez [MV/m]\n";
} else {
- fout << 5 << " Bx [kGauss]\n";
- fout << 6 << " By [kGauss]\n";
+ fout << 5 << " Bx [kGauss]\n";
+ fout << 6 << " By [kGauss]\n";
+ fout << 7 << " Bz [kGauss]\n";
+ fout << 8 << " Ex [MV/m]\n";
+ fout << 9 << " Ey [MV/m]\n";
+ fout << 10 << " Ez [MV/m]\n";
}
- fout << 7 << " Bz [kGauss]\n";
- if (coordinates_m == CYLINDRICAL) {
- fout << 5 << " Er [MV/m]\n";
- fout << 6 << " Ephi [MV/m]\n";
- } else {
- fout << 5 << " Ex [MV/m]\n";
- fout << 6 << " Ey [MV/m]\n";
- }
- fout << 10 << " Ez [MV/m]\n";
fout << 0 << std::endl;
}
diff --git a/src/Classic/AbsBeamline/BeamlineVisitor.h b/src/Classic/AbsBeamline/BeamlineVisitor.h
index 15c19cc95d8348f87ff1d03da2d089320e4c5113..1e1899c77ef93312fa32224d18963a32435d297c 100644
--- a/src/Classic/AbsBeamline/BeamlineVisitor.h
+++ b/src/Classic/AbsBeamline/BeamlineVisitor.h
@@ -54,6 +54,7 @@ class RBend;
class RBend3D;
class RFCavity;
class VariableRFCavity;
+class VariableRFCavityFringeField;
class TravelingWave;
class RFQuadrupole;
class SBend;
@@ -171,6 +172,10 @@ public:
/// Apply the algorithm to a variable RF cavity.
virtual void visitVariableRFCavity(const VariableRFCavity &) = 0;
+ /// Apply the algorithm to a variable RF cavity with Fringe Field.
+ virtual void visitVariableRFCavityFringeField
+ (const VariableRFCavityFringeField &) = 0;
+
/// Apply the algorithm to a RF cavity.
virtual void visitTravelingWave(const TravelingWave &) = 0;
diff --git a/src/Classic/AbsBeamline/CMakeLists.txt b/src/Classic/AbsBeamline/CMakeLists.txt
index 2e8dab69e68e481423d58cff3ffb92046f3e00c0..d5dfedfef8ac7dbfa68c231ef6652f7739c72922 100644
--- a/src/Classic/AbsBeamline/CMakeLists.txt
+++ b/src/Classic/AbsBeamline/CMakeLists.txt
@@ -43,6 +43,7 @@ set (_SRCS
Stripper.cpp
TravelingWave.cpp
VariableRFCavity.cpp
+ VariableRFCavityFringeField.cpp
)
include_directories (
@@ -96,6 +97,7 @@ set (HDRS
Stripper.h
TravelingWave.h
VariableRFCavity.h
+ VariableRFCavityFringeField.h
)
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline")
\ No newline at end of file
diff --git a/src/Classic/AbsBeamline/EndFieldModel/EndFieldModel.h b/src/Classic/AbsBeamline/EndFieldModel/EndFieldModel.h
index c2e0580669834bfee6f6b1ba3c6af94c6ca3e644..f276bd9971c7f2924704a7f656fbec12ddcc1eca 100644
--- a/src/Classic/AbsBeamline/EndFieldModel/EndFieldModel.h
+++ b/src/Classic/AbsBeamline/EndFieldModel/EndFieldModel.h
@@ -32,7 +32,7 @@
#include <vector>
namespace endfieldmodel {
-
+
class EndFieldModel {
public:
virtual ~EndFieldModel() {;}
diff --git a/src/Classic/AbsBeamline/EndFieldModel/Tanh.cpp b/src/Classic/AbsBeamline/EndFieldModel/Tanh.cpp
index 43ff3146535fae866ecfdb0a396eef76ff1ecff8..3aa3fed0c4ae5eff147fbb884688b12207ee51bd 100644
--- a/src/Classic/AbsBeamline/EndFieldModel/Tanh.cpp
+++ b/src/Classic/AbsBeamline/EndFieldModel/Tanh.cpp
@@ -104,5 +104,11 @@ std::vector< std::vector<int> > Tanh::getTanhDiffIndices(size_t n) {
return _tdi[n];
}
+
+std::ostream& Tanh::print(std::ostream& out) const {
+ out << "Tanh model with centre length: " << _x0
+ << " end length: " << _lambda;
+ return out;
+}
}
diff --git a/src/Classic/AbsBeamline/EndFieldModel/Tanh.h b/src/Classic/AbsBeamline/EndFieldModel/Tanh.h
index b3418eb46b632f3df48db46082824b6493bb61b9..7d441f2d32b45cc5d40881dab66d7a579660a795 100644
--- a/src/Classic/AbsBeamline/EndFieldModel/Tanh.h
+++ b/src/Classic/AbsBeamline/EndFieldModel/Tanh.h
@@ -101,7 +101,7 @@ class Tanh : public EndFieldModel {
inline void setX0(double x0) {_x0 = x0;}
/** Bug - does nothing */
- std::ostream& print(std::ostream& out) const {return out;}
+ std::ostream& print(std::ostream& out) const;
private:
double _x0, _lambda;
diff --git a/src/Classic/AbsBeamline/SpecificElementVisitor.h b/src/Classic/AbsBeamline/SpecificElementVisitor.h
index 62af615761fd52a80b135300bed27462fca86e71..6c2fe198c48f16accc372245e5bac3b8d04233aa 100644
--- a/src/Classic/AbsBeamline/SpecificElementVisitor.h
+++ b/src/Classic/AbsBeamline/SpecificElementVisitor.h
@@ -27,6 +27,7 @@
#include "AbsBeamline/Ring.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/VariableRFCavity.h"
+#include "AbsBeamline/VariableRFCavityFringeField.h"
#include "AbsBeamline/TravelingWave.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
@@ -139,6 +140,10 @@ public:
/// Apply the algorithm to a RF cavity.
virtual void visitVariableRFCavity(const VariableRFCavity &vcav);
+ /// Apply the algorithm to a RF cavity.
+ virtual void visitVariableRFCavityFringeField
+ (const VariableRFCavityFringeField &vcav);
+
/// Apply the algorithm to a RF cavity.
virtual void visitRFCavity(const RFCavity &);
@@ -348,10 +353,18 @@ void SpecificElementVisitor<ELEM>::visitRBend3D(const RBend3D &element) {
}
template<class ELEM>
-void SpecificElementVisitor<ELEM>::visitVariableRFCavity(const VariableRFCavity &element) {
+void SpecificElementVisitor<ELEM>::visitVariableRFCavity
+ (const VariableRFCavity &element) {
CastsTrait<ELEM, VariableRFCavity>::apply(allElementsOfTypeE, element);
}
+template<class ELEM>
+void SpecificElementVisitor<ELEM>::visitVariableRFCavityFringeField
+ (const VariableRFCavityFringeField &element) {
+ CastsTrait<ELEM, VariableRFCavityFringeField>::apply(allElementsOfTypeE,
+ element);
+}
+
template<class ELEM>
void SpecificElementVisitor<ELEM>::visitRFCavity(const RFCavity &element) {
CastsTrait<ELEM, RFCavity>::apply(allElementsOfTypeE, element);
diff --git a/src/Classic/AbsBeamline/VariableRFCavity.cpp b/src/Classic/AbsBeamline/VariableRFCavity.cpp
index 4cc23ebd3bd379f30e488f9f7a46b8f3beecf40a..41a681e4a6db669ca8a26a997e24a30ad4950332 100644
--- a/src/Classic/AbsBeamline/VariableRFCavity.cpp
+++ b/src/Classic/AbsBeamline/VariableRFCavity.cpp
@@ -55,13 +55,13 @@ VariableRFCavity& VariableRFCavity::operator=(const VariableRFCavity& rhs) {
setPhaseModel(NULL);
setAmplitudeModel(NULL);
setFrequencyModel(NULL);
- if (rhs._phase_td != NULL)
- setPhaseModel(std::shared_ptr<AbstractTimeDependence>(rhs._phase_td->clone()));
- if (rhs._amplitude_td != NULL) {
- setAmplitudeModel(std::shared_ptr<AbstractTimeDependence>(rhs._amplitude_td->clone()));
+ if (rhs.phaseTD_m != NULL)
+ setPhaseModel(std::shared_ptr<AbstractTimeDependence>(rhs.phaseTD_m->clone()));
+ if (rhs.amplitudeTD_m != NULL) {
+ setAmplitudeModel(std::shared_ptr<AbstractTimeDependence>(rhs.amplitudeTD_m->clone()));
}
- if (rhs._frequency_td != NULL)
- setFrequencyModel(std::shared_ptr<AbstractTimeDependence>(rhs._frequency_td->clone()));
+ if (rhs.frequencyTD_m != NULL)
+ setFrequencyModel(std::shared_ptr<AbstractTimeDependence>(rhs.frequencyTD_m->clone()));
phaseName_m = rhs.phaseName_m;
amplitudeName_m = rhs.amplitudeName_m;
frequencyName_m = rhs.frequencyName_m;
@@ -86,27 +86,27 @@ void VariableRFCavity::initNull() {
}
std::shared_ptr<AbstractTimeDependence> VariableRFCavity::getAmplitudeModel() const {
- return _amplitude_td;
+ return amplitudeTD_m;
}
std::shared_ptr<AbstractTimeDependence> VariableRFCavity::getPhaseModel() const {
- return _phase_td;
+ return phaseTD_m;
}
std::shared_ptr<AbstractTimeDependence> VariableRFCavity::getFrequencyModel() const {
- return _frequency_td;
+ return frequencyTD_m;
}
void VariableRFCavity::setAmplitudeModel(std::shared_ptr<AbstractTimeDependence> amplitude_td) {
- _amplitude_td = amplitude_td;
+ amplitudeTD_m = amplitude_td;
}
void VariableRFCavity::setPhaseModel(std::shared_ptr<AbstractTimeDependence> phase_td) {
- _phase_td = phase_td;
+ phaseTD_m = phase_td;
}
void VariableRFCavity::setFrequencyModel(std::shared_ptr<AbstractTimeDependence> frequency_td) {
- _frequency_td = frequency_td;
+ frequencyTD_m = frequency_td;
}
StraightGeometry &VariableRFCavity::getGeometry() {
@@ -146,9 +146,9 @@ bool VariableRFCavity::apply(const Vector_t &R, const Vector_t &P,
return true;
}
- double E0 = _amplitude_td->getValue(t);
- double f = _frequency_td->getValue(t) * 1.0E-3; // need GHz on the element we have MHz
- double phi = _phase_td->getValue(t);
+ double E0 = amplitudeTD_m->getValue(t);
+ double f = frequencyTD_m->getValue(t) * 1.0E-3; // need GHz on the element we have MHz
+ double phi = phaseTD_m->getValue(t);
E = Vector_t(0., 0., E0*sin(Physics::two_pi * f * t + phi));
return false;
}
@@ -173,6 +173,11 @@ ElementBase* VariableRFCavity::clone() const {
}
void VariableRFCavity::accept(BeamlineVisitor& visitor) const {
+ initialise();
+ visitor.visitVariableRFCavity(*this);
+}
+
+void VariableRFCavity::initialise() const {
VariableRFCavity* cavity = const_cast<VariableRFCavity*>(this);
std::shared_ptr<AbstractTimeDependence> phaseTD =
AbstractTimeDependence::getTimeDependence(phaseName_m);
@@ -183,7 +188,6 @@ void VariableRFCavity::accept(BeamlineVisitor& visitor) const {
std::shared_ptr<AbstractTimeDependence> amplitudeTD =
AbstractTimeDependence::getTimeDependence(amplitudeName_m);
cavity->setAmplitudeModel(std::shared_ptr<AbstractTimeDependence>(amplitudeTD->clone()));
- visitor.visitVariableRFCavity(*this);
if (halfHeight_m < 1e-9 || halfWidth_m < 1e-9)
throw GeneralClassicException("VariableRFCavity::accept",
diff --git a/src/Classic/AbsBeamline/VariableRFCavity.h b/src/Classic/AbsBeamline/VariableRFCavity.h
index 7a319151277909b3d501ff439aa4ce5b01ac1fcd..ead071aa6d29007624adb84326fc579e5fa66bdd 100644
--- a/src/Classic/AbsBeamline/VariableRFCavity.h
+++ b/src/Classic/AbsBeamline/VariableRFCavity.h
@@ -39,7 +39,7 @@ class Fieldmap;
*
* Generates a field like
* E = E0*a(t)*sin{f(t)*t-q(t)}
- * B = B0*a(t)*cos{f(t)*t-q(t)}
+ * B = 0
* where E0, B0 are user defined fields, a(t), f(t), q(t) are time
* dependent amplitude, frequency, phase respectively; it is assumed that these
* quantities vary sufficiently slowly that Maxwell is satisfied.
@@ -208,9 +208,11 @@ class VariableRFCavity: public Component {
/// Not implemented
virtual const EMField &getField() const;
protected:
- std::shared_ptr<AbstractTimeDependence> _phase_td;
- std::shared_ptr<AbstractTimeDependence> _amplitude_td;
- std::shared_ptr<AbstractTimeDependence> _frequency_td;
+ void initNull();
+ void initialise() const;
+ std::shared_ptr<AbstractTimeDependence> phaseTD_m;
+ std::shared_ptr<AbstractTimeDependence> amplitudeTD_m;
+ std::shared_ptr<AbstractTimeDependence> frequencyTD_m;
std::string phaseName_m;
std::string amplitudeName_m;
std::string frequencyName_m;
@@ -221,21 +223,19 @@ class VariableRFCavity: public Component {
StraightGeometry geometry;
static const double lengthUnit_m;
- private:
- void initNull();
-
+private:
};
double VariableRFCavity::getAmplitude(double time) const {
- return _amplitude_td->getValue(time);
+ return amplitudeTD_m->getValue(time);
}
double VariableRFCavity::getPhase(double time) const {
- return _phase_td->getValue(time);
+ return phaseTD_m->getValue(time);
}
double VariableRFCavity::getFrequency(double time) const {
- return _frequency_td->getValue(time);
+ return frequencyTD_m->getValue(time);
}
double VariableRFCavity::getHeight() const {
diff --git a/src/Classic/AbsBeamline/VariableRFCavityFringeField.cpp b/src/Classic/AbsBeamline/VariableRFCavityFringeField.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c5c6712a10a5d8c7c0affa2a0663ba03a4789a94
--- /dev/null
+++ b/src/Classic/AbsBeamline/VariableRFCavityFringeField.cpp
@@ -0,0 +1,233 @@
+/*
+ * Copyright (c) 2014, Chris Rogers
+ * All rights reserved.
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ * 3. Neither the name of STFC nor the names of its contributors may be used to
+ * endorse or promote products derived from this software without specific
+ * prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "Physics/Physics.h"
+#include "Algorithms/PartBunchBase.h"
+#include "AbsBeamline/BeamlineVisitor.h"
+
+#include "AbsBeamline/EndFieldModel/EndFieldModel.h"
+#include "AbsBeamline/VariableRFCavityFringeField.h"
+
+VariableRFCavityFringeField::VariableRFCavityFringeField(const std::string &name) : VariableRFCavity(name) {
+ initNull(); // initialise everything to NULL
+}
+
+VariableRFCavityFringeField::VariableRFCavityFringeField() : VariableRFCavity() {
+ initNull(); // initialise everything to NULL
+}
+
+VariableRFCavityFringeField::VariableRFCavityFringeField(const VariableRFCavityFringeField& var) : VariableRFCavity() {
+ initNull(); // initialise everything to NULL
+ *this = var;
+}
+
+VariableRFCavityFringeField& VariableRFCavityFringeField::operator=(const VariableRFCavityFringeField& rhs) {
+
+ if (&rhs == this) {
+ return *this;
+ }
+ VariableRFCavity::operator=(rhs);
+ setEndField(rhs.endField_m);
+ zCentre_m = rhs.zCentre_m;
+ f_m = rhs.f_m;
+ g_m = rhs.f_m;
+ h_m = rhs.f_m;
+ return *this;
+}
+
+VariableRFCavityFringeField::~VariableRFCavityFringeField() {
+}
+
+ElementBase* VariableRFCavityFringeField::clone() const {
+ return new VariableRFCavityFringeField(*this);
+}
+
+void VariableRFCavityFringeField::initNull() {
+ VariableRFCavity::initNull();
+ endField_m = std::shared_ptr<endfieldmodel::EndFieldModel>();
+ zCentre_m = 0;
+}
+
+void VariableRFCavityFringeField::accept(BeamlineVisitor& visitor) const {
+ VariableRFCavity::initialise();
+ VariableRFCavityFringeField* cavity =
+ const_cast<VariableRFCavityFringeField*>(this);
+ cavity->initialiseCoefficients();
+ visitor.visitVariableRFCavity(*this);
+}
+
+bool VariableRFCavityFringeField::apply(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B) {
+ if (R[2] > _length || R[2] < 0.) {
+ return true;
+ }
+ if (std::abs(R[0]) > halfWidth_m || std::abs(R[1]) > halfHeight_m) {
+ return true;
+ }
+ double z_pos = R[2]-zCentre_m;
+ double E0 = amplitudeTD_m->getValue(t);
+ double omega = Physics::two_pi*frequencyTD_m->getValue(t) * 1.0E-3; // need GHz on the element we have MHz
+ double phi = phaseTD_m->getValue(t);
+ double E_sin_t = E0*sin(omega * t + phi);
+ double B_cos_t = E0*cos(omega * t + phi); // 1/c^2 factor in the h_n coefficients
+
+ std::vector<double> y_power(maxOrder_m+1, 0.);
+ y_power[0] = 1.;
+ for (size_t i = 1; i < y_power.size(); ++i) {
+ y_power[i] = y_power[i-1]*R[1];
+ }
+
+ // d^i f0 dz^i
+ std::vector<double> endField(maxOrder_m/2+2, 0.);
+ for (size_t i = 0; i < endField.size(); ++i) {
+ endField[i] = endField_m->function(z_pos, i);
+ }
+
+ // omega^i
+ std::vector<double> omegaPower(maxOrder_m+1, 1.);
+ for (size_t i = 1; i < omegaPower.size(); ++i) {
+ omegaPower[i] = omegaPower[i-1]*omega;
+ }
+
+ E = Vector_t(0., 0., 0.);
+ B = Vector_t(0., 0., 0.);
+ // even power of y
+ for (size_t n = 0; n <= maxOrder_m ; n += 2) { // power of y
+ double fCoeff = 0.;
+ size_t index = n/2;
+ for (size_t i = 0; i < f_m[index].size(); i += 2) { // derivative of f
+ fCoeff += f_m[index][i]*endField[i]*omegaPower[n-i];
+ }
+ E[2] += E_sin_t*y_power[n]*fCoeff;
+ }
+ // odd power of y
+ for (size_t n = 1; n <= maxOrder_m; n += 2) {
+ double gCoeff = 0.;
+ double hCoeff = 0.;
+ size_t index = (n-1)/2;
+ //std::cerr << "n: " << n << std::endl;
+ for (size_t j = 0; j < g_m[index].size(); ++j) {
+ gCoeff += g_m[index][j]*endField[j]*omegaPower[n-j];
+ }
+ for (size_t j = 0; j < h_m[index].size(); ++j) {
+ hCoeff += h_m[index][j]*endField[j]*omegaPower[n-j];
+ //std::cerr << "j: " << j << " " << hCoeff << " ";
+ }
+ //std::cerr << std::endl;
+ E[1] += E_sin_t*y_power[n]*gCoeff;
+ B[0] += B_cos_t*y_power[n]*hCoeff;
+ //std::cerr << "APPLY B " << n << " " << B[0] << " " << hCoeff << std::endl;
+ }
+ B *= 1e4; //B natural units are kT; convert to kGauss
+ return false;
+}
+
+bool VariableRFCavityFringeField::apply(const size_t &i, const double &t,
+ Vector_t &E, Vector_t &B) {
+ return apply(RefPartBunch_m->R[i], RefPartBunch_m->P[i], t, E, B);
+}
+
+bool VariableRFCavityFringeField::applyToReferenceParticle(const Vector_t &R,
+ const Vector_t &P,
+ const double &t,
+ Vector_t &E,
+ Vector_t &B) {
+ return apply(R, P, t, E, B);
+}
+
+void VariableRFCavityFringeField::initialise(PartBunchBase<double, 3> *bunch,
+ double &startField,
+ double &endField) {
+ VariableRFCavity::initialise(bunch, startField, endField);
+ initialiseCoefficients();
+}
+
+void printVector(std::ostream& out, std::vector< std::vector<double> > vec) {
+ for (size_t i = 0; i < vec.size(); ++i) {
+ out << std::setw(3) << i;
+ for (size_t j = 0; j < vec[i].size(); ++j) {
+ out << " " << std::setw(14) << vec[i][j];
+ }
+ out << std::endl;
+ }
+}
+
+void VariableRFCavityFringeField::initialiseCoefficients() {
+ f_m = std::vector< std::vector<double> >();
+ g_m = std::vector< std::vector<double> >();
+ h_m = std::vector< std::vector<double> >();
+ f_m.push_back(std::vector<double>(1, 1.));
+ double c_l = Physics::c*1e-6; // speed of light in mm/ns
+ // generate f_{n+2} term
+ // note frequency term has to be added at apply(time) as we have
+ // time-dependent frequency
+ for(size_t n = 0; n+2 <= maxOrder_m; n += 2) {
+ // n denotes the subscript on f_n
+ // n+2 is the subscript on g_{n+2} and terms proportional to y^{n+2}
+ std::vector<double> f_n = f_m.back(); // f_n
+ std::vector<double> f_np2 = std::vector<double>(f_n.size()+2, 0.); // f_{n+2}
+ double n_const = 1./(n+1.)/(n+2.);
+ for (size_t j = 0; j < f_n.size(); ++j) {
+ f_np2[j] += f_n[j]*n_const/c_l/c_l;
+ }
+ for (size_t j = 0; j < f_n.size(); ++j) {
+ f_np2[j+2] += f_n[j]*n_const;
+ }
+ f_m.push_back(f_np2);
+ }
+ // generate g_{n+2} and h_{n+2} term
+ for(size_t n = 0; n+1 <= maxOrder_m; n += 2) {
+ // n denotes the subscript on f_n
+ // n+1 is the subscript on g_{n+1} and terms proportional to y^{n+1}
+ size_t f_index = n/2;
+ std::vector<double> f_n = f_m[f_index];
+ std::vector<double> g_np1 = std::vector<double>(f_n.size()+1, 0.);
+ std::vector<double> h_np1 = std::vector<double>(f_n.size(), 0.);
+ for (size_t j = 0; j < f_n.size(); ++j) {
+ g_np1[j+1] = -1./(n+1.)*f_n[j];
+ h_np1[j] = -1./c_l/c_l/(n+1.)*f_n[j];
+ }
+ g_m.push_back(g_np1);
+ h_m.push_back(h_np1);
+ }
+}
+
+void VariableRFCavityFringeField::printCoefficients(std::ostream& out) const {
+ out << "f_m" << std::endl;
+ printVector(out, f_m);
+ out << "g_m" << std::endl;
+ printVector(out, g_m);
+ out << "h_m" << std::endl;
+ printVector(out, h_m);
+ out << std::endl;
+}
+
+
+
+void VariableRFCavityFringeField::setEndField(
+ std::shared_ptr<endfieldmodel::EndFieldModel> end) {
+ endField_m = end;
+}
diff --git a/src/Classic/AbsBeamline/VariableRFCavityFringeField.h b/src/Classic/AbsBeamline/VariableRFCavityFringeField.h
new file mode 100644
index 0000000000000000000000000000000000000000..a0480e820251f51c546fdf789a63ce281d6602f8
--- /dev/null
+++ b/src/Classic/AbsBeamline/VariableRFCavityFringeField.h
@@ -0,0 +1,210 @@
+/*
+ * Copyright (c) 2014, Chris Rogers
+ * All rights reserved.
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ * 3. Neither the name of STFC nor the names of its contributors may be used to
+ * endorse or promote products derived from this software without specific
+ * prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef CLASSIC_ABSBEAMLINE_VariableRFCavityFringeField_HH
+#define CLASSIC_ABSBEAMLINE_VariableRFCavityFringeField_HH
+
+#include <iostream>
+#include "AbsBeamline/VariableRFCavity.h"
+
+class Fieldmap;
+namespace endfieldmodel {
+ class EndFieldModel;
+}
+
+/** @class VariableRFCavityFringeField
+ *
+ * Generates a field like
+ * Ey = E0*a(t)*y^{2n+1} g_n(z) sin{f(t)*t-q(t)}
+ * Ez = E0*a(t)*y^{2n} f_n(z) sin{f(t)*t-q(t)}
+ * Bx = B0*a(t)*y^{2n+1} h_n(z) cos{f(t)*t-q(t)}
+ * where E0, is a user-defined field and B0 is the corresponding magnetic field
+ * magnitude; f_n is a user-defined axial field dependence and g_n, h_n are the
+ * corresponding off-axis field dependencies. a(t), f(t), q(t) are time
+ * dependent amplitude, frequency, phase respectively; it is assumed that these
+ * quantities vary sufficiently slowly that Maxwell is satisfied.
+ *
+ * Inherits from VariableRFCavity; inheritance is used here to share code from
+ * VariableRFCavity.
+ *
+ * Set/get methods use metres; but internally we store as mm (for tracking)
+ *
+ * Field units are kG and GV/mm
+ */
+class VariableRFCavityFringeField : public VariableRFCavity {
+ public:
+ /// Constructor with given name.
+ explicit VariableRFCavityFringeField(const std::string &name);
+ /** Copy Constructor; performs deepcopy on time-dependence models */
+ VariableRFCavityFringeField(const VariableRFCavityFringeField &);
+ /** Default constructor */
+ VariableRFCavityFringeField();
+ /** Assignment operator; performs deepcopy on time-dependence models*/
+ VariableRFCavityFringeField& operator=(const VariableRFCavityFringeField &);
+ /** Destructor does nothing
+ *
+ * The shared_ptrs will self-destruct when reference count goes to 0
+ */
+ virtual ~VariableRFCavityFringeField();
+
+ /** Apply visitor to RFCavity.
+ *
+ * The RF cavity finds the "time dependence" models by doing a string
+ * lookup against a list held by AbstractTimeDependence at accept time.
+ */
+ virtual void accept(BeamlineVisitor &) const;
+
+ /** Inheritable deepcopy method */
+ virtual ElementBase* clone() const;
+
+ /** Calculate the field at the position of the i^th particle
+ *
+ * @param i indexes the particle whose field we need
+ * @param t the time at which the field is calculated
+ * @param E return value with electric field strength
+ * @param B return value with magnetic field strength
+ *
+ * @returns True if particle is outside the boundaries; else False
+ */
+ virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
+
+ /** Calculate the field at a given position
+ *
+ * @param R the position at which the field is calculated
+ * @param P the momentum (not used)
+ * @param t the time at which the field is calculated
+ * @param E return value; filled with electric field strength
+ * @param B return value; filled with magnetic field strength
+ *
+ * @returns True if particle is outside the boundaries; else False
+ */
+ virtual bool apply(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B);
+
+
+ /** Calculate the field at a given position. This is identical to "apply".
+ *
+ * @param R the position at which the field is calculated
+ * @param P the momentum (not used)
+ * @param t the time at which the field is calculated
+ * @param E return value; filled with electric field strength
+ * @param B return value; filled with magnetic field strength
+ *
+ * @returns True if particle is outside the boundaries; else False
+ */
+ virtual bool applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B);
+
+ /** Set the endFieldModel
+ *
+ * @param endField the new endFieldModel. VariableRFCavityFringe takes
+ * ownership of the memory allocated to endField.
+ */
+ virtual void setEndField
+ (std::shared_ptr<endfieldmodel::EndFieldModel> endField);
+
+ /** Get the endFieldModel
+ *
+ * @returns the endFieldModel; VariableRFCavityFringe retains ownership of
+ * the memory allocated to the endFieldModel.
+ */
+ virtual inline std::shared_ptr<endfieldmodel::EndFieldModel>
+ getEndField() const;
+
+ /** Initialise ready for tracking
+ *
+ * Does any setup on VirtualRFCavity and sets field expansion coefficients
+ */
+ virtual void initialise(PartBunchBase<double, 3> *bunch,
+ double &startField,
+ double &endField);
+
+ /** Get the offset of centre of the cavity field from the element start in metres */
+ virtual inline void setCavityCentre(double zCentre);
+ /** Set the offset of centre of the cavity field from the element start in metres */
+ virtual inline double getCavityCentre() const;
+ /** Set the maximum y power */
+ virtual inline void setMaxOrder(size_t maxOrder);
+ /** Get the maximum y power that will be used in field calculations */
+ virtual inline size_t getMaxOrder() const;
+ /** Set the coefficients for calculating the field expansion */
+ void initialiseCoefficients();
+ /** Print the coefficients to ostream out */
+ void printCoefficients(std::ostream& out) const;
+ /** Get the coefficients for Ez */
+ inline std::vector<std::vector<double> > getEzCoefficients() const;
+ /** Get the coefficients for Ey */
+ inline std::vector<std::vector<double> > getEyCoefficients() const;
+ /** Get the coefficients for Bx */
+ inline std::vector<std::vector<double> > getBxCoefficients() const;
+protected:
+ double zCentre_m; // offsets this element
+ size_t maxOrder_m;
+ std::shared_ptr<endfieldmodel::EndFieldModel> endField_m;
+ std::vector<std::vector<double> > f_m;
+ std::vector<std::vector<double> > g_m;
+ std::vector<std::vector<double> > h_m;
+private:
+ void initNull();
+};
+
+void VariableRFCavityFringeField::setCavityCentre(double zCentre) {
+ zCentre_m = zCentre*1000.; // stored internally as mm
+}
+
+double VariableRFCavityFringeField::getCavityCentre() const {
+ return zCentre_m/1000.; // stored internally as mm
+}
+
+void VariableRFCavityFringeField::setMaxOrder(size_t maxOrder) {
+ maxOrder_m = maxOrder;
+ initialiseCoefficients();
+}
+
+size_t VariableRFCavityFringeField::getMaxOrder() const {
+ return maxOrder_m;
+}
+
+std::shared_ptr<endfieldmodel::EndFieldModel>
+ VariableRFCavityFringeField::getEndField() const {
+ return endField_m;
+}
+
+std::vector<std::vector<double> >
+ VariableRFCavityFringeField::getEzCoefficients() const {
+ return f_m;
+}
+
+std::vector<std::vector<double> >
+ VariableRFCavityFringeField::getEyCoefficients() const {
+ return g_m;
+}
+
+std::vector<std::vector<double> >
+ VariableRFCavityFringeField::getBxCoefficients() const {
+ return h_m;
+}
+
+#endif // #ifdef CLASSIC_VirtualRFCavityFringeField_HH
\ No newline at end of file
diff --git a/src/Classic/Algorithms/AbstractTimeDependence.cpp b/src/Classic/Algorithms/AbstractTimeDependence.cpp
index 118634d80a53edee2f894c026de43bccb32e4d25..188230a5f43e773cb14be0c1514072502d6121a7 100644
--- a/src/Classic/Algorithms/AbstractTimeDependence.cpp
+++ b/src/Classic/Algorithms/AbstractTimeDependence.cpp
@@ -1,3 +1,30 @@
+/*
+ * Copyright (c) 2014, Chris Rogers
+ * All rights reserved.
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ * 3. Neither the name of STFC nor the names of its contributors may be used to
+ * endorse or promote products derived from this software without specific
+ * prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
#include "Algorithms/AbstractTimeDependence.h"
#include "Utilities/GeneralClassicException.h"
@@ -17,9 +44,6 @@ std::shared_ptr<AbstractTimeDependence> AbstractTimeDependence::getTimeDependenc
void AbstractTimeDependence::setTimeDependence(std::string name,
std::shared_ptr<AbstractTimeDependence> time_dep) {
- // if (td_map.find(name) != td_map.end()) {
- // delete td_map[name];
- // }
td_map[name] = time_dep;
}
diff --git a/src/Classic/Algorithms/AbstractTimeDependence.h b/src/Classic/Algorithms/AbstractTimeDependence.h
index f63f34ff5c3f87fd8fdfe7f7d304a30b9905692b..d36671870c0eb70877d5eaf7251677c650ba3c67 100644
--- a/src/Classic/Algorithms/AbstractTimeDependence.h
+++ b/src/Classic/Algorithms/AbstractTimeDependence.h
@@ -1,3 +1,30 @@
+/*
+ * Copyright (c) 2014, Chris Rogers
+ * All rights reserved.
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ * 3. Neither the name of STFC nor the names of its contributors may be used to
+ * endorse or promote products derived from this software without specific
+ * prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
#ifndef _CLASSIC_SRC_ALGORITHMS_ABSTRACTTIMEDEPENDENCE_H_
#define _CLASSIC_SRC_ALGORITHMS_ABSTRACTTIMEDEPENDENCE_H_
@@ -5,22 +32,64 @@
#include <map>
#include <memory>
-/** Time dependence abstraction for field parameters that vary slowly with time
+/** @class AbstractTimeDependence
+ *
+ * Time dependence abstraction for field parameters that vary slowly with time;
+ * for example, in RF cavities for synchrotrons the RF frequency varies to
+ * match the time-of-flight of the particles
+ *
+ * This base class stores a mapping of string to time dependence. At Visit
+ * time, we do a map look-up to assign the appropriate TimeDependence to the
+ * relevant field Elements.
*/
class AbstractTimeDependence {
// Unit tests for this class are in with PolynomialTimeDependence
public:
+ /** Destructor does nothing */
virtual ~AbstractTimeDependence() {}
- virtual AbstractTimeDependence* clone() = 0;
+ /** Inheritable copy constructor
+ *
+ * @returns new AbstractTimeDependence that is a copy of this. User owns
+ * returned memory.
+ */
+ virtual AbstractTimeDependence* clone() = 0;
+
+ /** getValue(time) returns the value as a function of time.
+ *
+ * This could represent RF voltage or frequency, magnetic field strength
+ * etc
+ */
virtual double getValue(double time) = 0;
- static std::shared_ptr<AbstractTimeDependence> getTimeDependence(std::string name);
+ /** Look up the time dependence that has a given name
+ *
+ * @param name name of the time dependence
+ *
+ * @returns shared_ptr to the appropriate time dependence.
+ * @throws GeneralClassicException if name is not recognised
+ */
+ static std::shared_ptr<AbstractTimeDependence>
+ getTimeDependence(std::string name);
+
+ /** Add a value to the lookup table
+ *
+ * @param name name of the time dependence. If name already exists in the
+ * map, it is overwritten with the new value.
+ * @param time_dep shared_ptr to the time dependence.
+ */
static void setTimeDependence(std::string name,
- std::shared_ptr<AbstractTimeDependence> time_dep);
+ std::shared_ptr<AbstractTimeDependence> time_dep);
- /** Reverse the mapping - slow
+ /** Get the name corresponding to a given time_dep
+ *
+ * @param time_dep time dependence to lookup
+ *
+ * @returns name corresponding to the time dependence. Note that this
+ * just does a dumb loop over the stored map values; so O(N).
+ * @throws GeneralClassicException if time_dep is not recognised
*/
- static std::string getName(std::shared_ptr<AbstractTimeDependence> time_dep);
+ static std::string getName
+ (std::shared_ptr<AbstractTimeDependence> time_dep);
private:
static std::map<std::string, std::shared_ptr<AbstractTimeDependence> > td_map;
diff --git a/src/Classic/Algorithms/DefaultVisitor.cpp b/src/Classic/Algorithms/DefaultVisitor.cpp
index 2125b2f0f51c84d4bdfc24dc3a7c9b96bf073cd9..d47e6d0a0d710b33ec1fb71b0a85a18cd8597227 100644
--- a/src/Classic/Algorithms/DefaultVisitor.cpp
+++ b/src/Classic/Algorithms/DefaultVisitor.cpp
@@ -42,6 +42,7 @@
#include "AbsBeamline/RBend3D.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/VariableRFCavity.h"
+#include "AbsBeamline/VariableRFCavityFringeField.h"
#include "AbsBeamline/TravelingWave.h"
#include "AbsBeamline/RFQuadrupole.h"
#include "AbsBeamline/SBend.h"
@@ -180,6 +181,11 @@ void DefaultVisitor::visitVariableRFCavity(const VariableRFCavity &vcav) {
applyDefault(vcav);
}
+void DefaultVisitor::visitVariableRFCavityFringeField
+ (const VariableRFCavityFringeField &vcav) {
+ applyDefault(vcav);
+}
+
void DefaultVisitor::visitRFCavity(const RFCavity &cav) {
applyDefault(cav);
}
diff --git a/src/Classic/Algorithms/DefaultVisitor.h b/src/Classic/Algorithms/DefaultVisitor.h
index 3006f0e9381430d99e654d05229abf2e85266166..270ff1d6d33835277c3c5a8d33e52a1519da70de 100644
--- a/src/Classic/Algorithms/DefaultVisitor.h
+++ b/src/Classic/Algorithms/DefaultVisitor.h
@@ -115,6 +115,9 @@ public:
/// Apply the algorithm to a RF cavity.
virtual void visitVariableRFCavity(const VariableRFCavity &vcav);
+ /// Apply the algorithm to a RF cavity.
+ virtual void visitVariableRFCavityFringeField(const VariableRFCavityFringeField &vcav);
+
/// Apply the algorithm to a RF cavity.
virtual void visitRFCavity(const RFCavity &);
diff --git a/src/Classic/Algorithms/PolynomialTimeDependence.cpp b/src/Classic/Algorithms/PolynomialTimeDependence.cpp
index d0d270ac977457b4b6400a72ae4794c710a64fab..00a07a7bffbec68b1785b771f2d03c3087be900d 100644
--- a/src/Classic/Algorithms/PolynomialTimeDependence.cpp
+++ b/src/Classic/Algorithms/PolynomialTimeDependence.cpp
@@ -1,3 +1,30 @@
+/*
+ * Copyright (c) 2014, Chris Rogers
+ * All rights reserved.
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ * 3. Neither the name of STFC nor the names of its contributors may be used to
+ * endorse or promote products derived from this software without specific
+ * prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
#include "Algorithms/PolynomialTimeDependence.h"
Inform &PolynomialTimeDependence::print(Inform &os) {
diff --git a/src/Classic/Algorithms/PolynomialTimeDependence.h b/src/Classic/Algorithms/PolynomialTimeDependence.h
index f7d11aa9fb58ccbeae0eb3813f14f9bc824824ad..ac2e7abaa31ee5dc91cf7eb83241169dbd1fd016 100644
--- a/src/Classic/Algorithms/PolynomialTimeDependence.h
+++ b/src/Classic/Algorithms/PolynomialTimeDependence.h
@@ -1,3 +1,30 @@
+/*
+ * Copyright (c) 2014, Chris Rogers
+ * All rights reserved.
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ * 3. Neither the name of STFC nor the names of its contributors may be used to
+ * endorse or promote products derived from this software without specific
+ * prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
#ifndef _CLASSIC_SRC_ALGORITHMS_POLYNOMIALTIMEDEPENDENCE_H_
#define _CLASSIC_SRC_ALGORITHMS_POLYNOMIALTIMEDEPENDENCE_H_
@@ -6,18 +33,42 @@
#include "Ippl.h"
#include "Algorithms/AbstractTimeDependence.h"
+
+/** @class PolynomialTimeDependence
+ *
+ * Time dependence that follows a polynomial, like
+ * p_0 + p_1*t + p_2*t^2 + ... + p_i*t^i + ...
+ */
class PolynomialTimeDependence : public AbstractTimeDependence {
public:
+ /** Constructor
+ *
+ * @param ptd the polynomial coefficients p_i; can be of arbitrary length
+ * (user is responsible for issues like floating point precision).
+ */
PolynomialTimeDependence(std::vector<double> ptd) : coeffs(ptd) {}
+
+ /** Default Constructor makes a 0 length polynomial */
PolynomialTimeDependence() {}
+ /** Destructor does nothing */
~PolynomialTimeDependence() {}
+ /** Return the polynomial Sum_i p_i t^i; returns 0 if p is of 0 length */
inline double getValue(double time);
+ /** Inheritable copy constructor
+ *
+ * @returns new PolynomialTimeDependence that is a copy of this. User owns
+ * returned memory.
+ */
PolynomialTimeDependence* clone() {
std::vector<double> temp(coeffs);
PolynomialTimeDependence* d = new PolynomialTimeDependence(temp);
return d;
}
+ /** Print the polynomials
+ *
+ * @param os "Inform" stream to which the polynomials are printed.
+ */
Inform &print(Inform &os);
private:
diff --git a/src/Elements/CMakeLists.txt b/src/Elements/CMakeLists.txt
index 826582a49f5e4b3ac462f1e1acb873e7a303bdff..5677cc467488f6979ab1c134a30cfd604253f630 100644
--- a/src/Elements/CMakeLists.txt
+++ b/src/Elements/CMakeLists.txt
@@ -45,6 +45,7 @@ set (_SRCS
OpalSource.cpp
OpalSRot.cpp
OpalVariableRFCavity.cpp
+ OpalVariableRFCavityFringeField.cpp
OpalVKicker.cpp
OpalVMonitor.cpp
OpalYRot.cpp
@@ -105,6 +106,7 @@ set (HDRS
OpalStripper.h
OpalTravelingWave.h
OpalVariableRFCavity.h
+ OpalVariableRFCavityFringeField.h
OpalVKicker.h
OpalVMonitor.h
OpalYRot.h
diff --git a/src/Elements/OpalVariableRFCavityFringeField.cpp b/src/Elements/OpalVariableRFCavityFringeField.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..893cc427007dda0b50995ef7b72f0c5608b65d79
--- /dev/null
+++ b/src/Elements/OpalVariableRFCavityFringeField.cpp
@@ -0,0 +1,188 @@
+/*
+ * Copyright (c) 2018, Chris Rogers
+ * All rights reserved.
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ * 3. Neither the name of STFC nor the names of its contributors may be used to
+ * endorse or promote products derived from this software without specific
+ * prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "Physics/Physics.h"
+#include "Utilities/OpalException.h"
+#include "Attributes/Attributes.h"
+#include "Algorithms/AbstractTimeDependence.h"
+#include "Utilities/OpalException.h"
+#include "AbsBeamline/EndFieldModel/Tanh.h"
+#include "AbsBeamline/VariableRFCavityFringeField.h"
+#include "Elements/OpalVariableRFCavityFringeField.h"
+
+extern Inform *gmsg;
+
+const std::string OpalVariableRFCavityFringeField::doc_string =
+ std::string("The \"VARIABLE_RF_CAVITY_FRINGE_FIELD\" element defines an RF cavity ")+
+ std::string("with time dependent frequency, phase and amplitude.");
+
+OpalVariableRFCavityFringeField::OpalVariableRFCavityFringeField():
+ OpalElement(SIZE, "VARIABLE_RF_CAVITY_FRINGE_FIELD", doc_string.c_str()) {
+ itsAttr[PHASE_MODEL] = Attributes::makeString("PHASE_MODEL",
+ "The name of the phase time dependence model, which should give the phase in [rad].");
+ itsAttr[AMPLITUDE_MODEL] = Attributes::makeString("AMPLITUDE_MODEL",
+ "The name of the amplitude time dependence model, which should give the field in [MV/m]");
+ itsAttr[FREQUENCY_MODEL] = Attributes::makeString("FREQUENCY_MODEL",
+ "The name of the frequency time dependence model, which should give the field in [MHz].");
+ itsAttr[WIDTH] = Attributes::makeReal("WIDTH",
+ "Full width of the cavity [m].");
+ itsAttr[HEIGHT] = Attributes::makeReal("HEIGHT",
+ "Full height of the cavity [m].");
+ itsAttr[CENTRE_LENGTH] = Attributes::makeReal("CENTRE_LENGTH",
+ "Length of the cavity field flat top [m].");
+ itsAttr[END_LENGTH] = Attributes::makeReal("END_LENGTH",
+ "Length of the cavity fringe fields [m].");
+ itsAttr[CAVITY_CENTRE] = Attributes::makeReal("CAVITY_CENTRE",
+ "Offset of the cavity centre from the beginning of the cavity [m].");
+ itsAttr[MAX_ORDER] = Attributes::makeReal("MAX_ORDER",
+ "Maximum power of y that will be evaluated in field calculations.");
+
+ registerStringAttribute("PHASE_MODEL");
+ registerStringAttribute("AMPLITUDE_MODEL");
+ registerStringAttribute("FREQUENCY_MODEL");
+ registerRealAttribute("WIDTH");
+ registerRealAttribute("HEIGHT");
+ registerRealAttribute("CENTRE_LENGTH");
+ registerRealAttribute("END_LENGTH");
+ registerRealAttribute("CAVITY_CENTRE");
+ registerRealAttribute("MAX_ORDER");
+
+ registerOwnership();
+
+ setElement((new VariableRFCavityFringeField("VARIABLE_RF_CAVITY_FRINGE_FIELD"))->
+ makeAlignWrapper());
+}
+
+OpalVariableRFCavityFringeField::OpalVariableRFCavityFringeField(
+ const std::string &name,
+ OpalVariableRFCavityFringeField *parent
+ ) : OpalElement(name, parent) {
+ VariableRFCavityFringeField *cavity = dynamic_cast
+ <VariableRFCavityFringeField*>(parent->getElement()->removeWrappers());
+ setElement((new VariableRFCavityFringeField(*cavity))->makeAlignWrapper());
+}
+
+OpalVariableRFCavityFringeField::~OpalVariableRFCavityFringeField() {
+}
+
+OpalVariableRFCavityFringeField *OpalVariableRFCavityFringeField::clone(
+ const std::string &name) {
+ return new OpalVariableRFCavityFringeField(name, this);
+}
+
+OpalVariableRFCavityFringeField *OpalVariableRFCavityFringeField::clone() {
+ return new OpalVariableRFCavityFringeField(this->getOpalName(), this);
+}
+
+void OpalVariableRFCavityFringeField::
+fillRegisteredAttributes(const ElementBase &base, ValueFlag flag) {
+ OpalElement::fillRegisteredAttributes(base, flag);
+ const VariableRFCavityFringeField* cavity =
+ dynamic_cast<const VariableRFCavityFringeField*>(&base);
+ if (cavity == NULL) {
+ throw OpalException("OpalVariableRFCavityFringeField::fillRegisteredAttributes",
+ "Failed to cast ElementBase to a VariableRFCavityFringeField");
+ }
+ std::shared_ptr<endfieldmodel::EndFieldModel> model = cavity->getEndField();
+ endfieldmodel::Tanh* tanh = dynamic_cast<endfieldmodel::Tanh*>(model.get());
+ if (tanh == NULL) {
+ throw OpalException("OpalVariableRFCavityFringeField::fillRegisteredAttributes",
+ "Failed to cast EndField to a Tanh model");
+ }
+
+
+ attributeRegistry["L"]->setReal(cavity->getLength());
+ std::shared_ptr<AbstractTimeDependence> phase_model = cavity->getPhaseModel();
+ std::shared_ptr<AbstractTimeDependence> freq_model = cavity->getFrequencyModel();
+ std::shared_ptr<AbstractTimeDependence> amp_model = cavity->getAmplitudeModel();
+ std::string phase_name = AbstractTimeDependence::getName(phase_model);
+ std::string amp_name = AbstractTimeDependence::getName(amp_model);
+ std::string freq_name = AbstractTimeDependence::getName(freq_model);
+ attributeRegistry["PHASE_MODEL"]->setString(phase_name);
+ attributeRegistry["AMPLITUDE_MODEL"]->setString(amp_name);
+ attributeRegistry["FREQUENCY_MODEL"]->setString(freq_name);
+ attributeRegistry["WIDTH"]->setReal(cavity->getWidth());
+ attributeRegistry["HEIGHT"]->setReal(cavity->getHeight());
+ // flat top length is 2*x0
+ attributeRegistry["CENTRE_LENGTH"]->setReal(tanh->getX0()/2.);
+ attributeRegistry["END_LENGTH"]->setReal(tanh->getLambda());
+ attributeRegistry["CAVITY_CENTRE"]->setReal(cavity->getCavityCentre());
+ attributeRegistry["MAX_ORDER"]->setReal(cavity->getMaxOrder());
+}
+
+
+
+void OpalVariableRFCavityFringeField::update() {
+ OpalElement::update();
+
+ VariableRFCavityFringeField *cavity = dynamic_cast
+ <VariableRFCavityFringeField*>(getElement()->removeWrappers());
+ double length = Attributes::getReal(itsAttr[LENGTH]);
+ cavity->setLength(length);
+ std::string phaseName = Attributes::getString(itsAttr[PHASE_MODEL]);
+ cavity->setPhaseName(phaseName);
+ std::string ampName = Attributes::getString(itsAttr[AMPLITUDE_MODEL]);
+ cavity->setAmplitudeName(ampName);
+ std::string freqName = Attributes::getString(itsAttr[FREQUENCY_MODEL]);
+ cavity->setFrequencyName(freqName);
+ double width = Attributes::getReal(itsAttr[WIDTH]);
+ cavity->setWidth(width);
+ double height = Attributes::getReal(itsAttr[HEIGHT]);
+ cavity->setHeight(height);
+ double maxOrderReal = Attributes::getReal(itsAttr[MAX_ORDER]);
+ size_t maxOrder = convertToUnsigned(maxOrderReal, "MAX_ORDER");
+ cavity->setMaxOrder(maxOrder);
+ double cavity_centre = Attributes::getReal(itsAttr[CAVITY_CENTRE]);
+ cavity->setCavityCentre(cavity_centre); // mm
+ // convert to mm; x0 is double length of flat top so divide 2
+ double centreLength = Attributes::getReal(itsAttr[CENTRE_LENGTH])*1e3;
+ double endLength = Attributes::getReal(itsAttr[END_LENGTH])*1e3;
+ endfieldmodel::Tanh* tanh = new endfieldmodel::Tanh(centreLength/2.,
+ endLength,
+ (maxOrder+1)/2);
+ std::shared_ptr<endfieldmodel::EndFieldModel> end(tanh);
+ cavity->setEndField(end);
+
+ setElement(cavity->makeAlignWrapper());
+}
+
+
+size_t OpalVariableRFCavityFringeField::convertToUnsigned(double value,
+ std::string name) {
+ value += unsignedTolerance; // prevent rounding error
+ if (fabs(floor(value) - value) > 2*unsignedTolerance) {
+ throw OpalException("OpalVariableRFCavityFringeField::convertToUnsigned",
+ "Value for "+name+
+ " should be an unsigned int but a real value was found");
+ }
+ if (floor(value) < -0.5) {
+ throw OpalException("OpalVariableRFCavityFringeField::convertToUnsigned",
+ "Value for "+name+" should be 0 or more");
+ }
+ size_t ret(floor(value));
+ return ret;
+}
\ No newline at end of file
diff --git a/src/Elements/OpalVariableRFCavityFringeField.h b/src/Elements/OpalVariableRFCavityFringeField.h
new file mode 100644
index 0000000000000000000000000000000000000000..192e8c6bc8a149f0492312d0cca327de958d309d
--- /dev/null
+++ b/src/Elements/OpalVariableRFCavityFringeField.h
@@ -0,0 +1,94 @@
+/*
+ * Copyright (c) 2018, Chris Rogers
+ * All rights reserved.
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ * 3. Neither the name of STFC nor the names of its contributors may be used to
+ * endorse or promote products derived from this software without specific
+ * prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef OPAL_OPALVARIABLERFCAVITYFRINGEFIELD_H
+#define OPAL_OPALVARIABLERFCAVITYFRINGEFIELD_H
+
+#include "Elements/OpalElement.h"
+
+/** OpalVariableRFCavityFringeField
+ */
+class OpalVariableRFCavityFringeField : public OpalElement {
+ public:
+ /** enum maps string to integer value for UI definitions */
+ enum {
+ PHASE_MODEL = COMMON,
+ AMPLITUDE_MODEL,
+ FREQUENCY_MODEL,
+ WIDTH,
+ HEIGHT,
+ CENTRE_LENGTH,
+ END_LENGTH,
+ CAVITY_CENTRE,
+ MAX_ORDER,
+ SIZE // size of the enum
+ };
+
+ /** Copy constructor **/
+ OpalVariableRFCavityFringeField(const std::string &name,
+ OpalVariableRFCavityFringeField *parent);
+
+ /** Default constructor **/
+ OpalVariableRFCavityFringeField();
+
+ /** Inherited copy constructor
+ *
+ * Call on a base class to instantiate an object of derived class's type
+ **/
+ OpalVariableRFCavityFringeField* clone();
+
+ /** Inherited copy constructor
+ *
+ * Call on a base class to instantiate an object of derived class's type
+ */
+ virtual OpalVariableRFCavityFringeField *clone(const std::string &name);
+
+ /** Destructor does nothing */
+ virtual ~OpalVariableRFCavityFringeField();
+
+ /** Fill in all registered attributes
+ *
+ * This updates the registered attributed with values from the ElementBase
+ */
+ virtual void fillRegisteredAttributes(const ElementBase &, ValueFlag);
+
+ /** Update the OpalVariableRFCavity with new parameters from UI parser */
+ virtual void update();
+ private:
+ // Not implemented.
+ OpalVariableRFCavityFringeField(const OpalVariableRFCavityFringeField &);
+ void operator=(const OpalVariableRFCavityFringeField &);
+
+ /** Check for conversion to unsigned int */
+ inline static size_t convertToUnsigned(double value, std::string name);
+
+ static const std::string doc_string;
+ static constexpr double unsignedTolerance = 1e-9;
+};
+
+
+#endif // OPAL_OPALVARIABLERFCAVITY_H
diff --git a/src/OpalConfigure/Configure.cpp b/src/OpalConfigure/Configure.cpp
index 93d0c942469c6293be4a998ae26a6d5b070fb3c0..08da087b89442a5a2c479d23bd18f5731dc03c4f 100644
--- a/src/OpalConfigure/Configure.cpp
+++ b/src/OpalConfigure/Configure.cpp
@@ -137,6 +137,7 @@
#include "Elements/OpalStripper.h"
#include "Elements/OpalRingDefinition.h"
#include "Elements/OpalVariableRFCavity.h"
+#include "Elements/OpalVariableRFCavityFringeField.h"
// Structure-related commands.
#include "Lines/Line.h"
@@ -278,6 +279,7 @@ namespace {
opal->create(new OpalSRot());
opal->create(new OpalTravelingWave());
opal->create(new OpalVariableRFCavity());
+ opal->create(new OpalVariableRFCavityFringeField());
opal->create(new OpalVKicker());
opal->create(new OpalVMonitor());
// opal->create(new OpalWire());
diff --git a/tests/classic_src/AbsBeamline/CMakeLists.txt b/tests/classic_src/AbsBeamline/CMakeLists.txt
index 0dfd61844fd510a059dbe1ee7ac09887d378fdae..78a4b0fd05bfa4597cf5951c848bafa38ebef067 100644
--- a/tests/classic_src/AbsBeamline/CMakeLists.txt
+++ b/tests/classic_src/AbsBeamline/CMakeLists.txt
@@ -7,6 +7,7 @@ set (_SRCS
ScalingFFAGMagnetTest.cpp
TrimCoilTest.cpp
VariableRFCavityTest.cpp
+ VariableRFCavityFringeFieldTest.cpp
)
include_directories (
diff --git a/tests/classic_src/AbsBeamline/VariableRFCavityFringeFieldTest.cpp b/tests/classic_src/AbsBeamline/VariableRFCavityFringeFieldTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b23f5c586f731f4565d6b8af2c60190cbbcf53f8
--- /dev/null
+++ b/tests/classic_src/AbsBeamline/VariableRFCavityFringeFieldTest.cpp
@@ -0,0 +1,304 @@
+/*
+ * Copyright (c) 2014, Chris Rogers
+ * All rights reserved.
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ * 3. Neither the name of STFC nor the names of its contributors may be used to
+ * endorse or promote products derived from this software without specific
+ * prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include <memory>
+
+#include "gtest/gtest.h"
+
+#include "opal_test_utilities/SilenceTest.h"
+#include "Physics/Physics.h"
+
+#include "Algorithms/PolynomialTimeDependence.h"
+#include "AbsBeamline/EndFieldModel/Tanh.h"
+#include "AbsBeamline/VariableRFCavityFringeField.h"
+
+class VariableRFCavityFringeFieldTest : public ::testing::Test {
+public:
+ VariableRFCavityFringeFieldTest() {
+ // OpalTestUtilities::SilenceTest silencer;
+
+ cav1 = VariableRFCavityFringeField("bob");
+ // centre, end, max_order
+ shared = std::shared_ptr<endfieldmodel::EndFieldModel>
+ (new endfieldmodel::Tanh(250., 50., 10));
+ cav1.setCavityCentre(0.500);
+ cav1.setEndField(shared);
+ PolynomialTimeDependence* time =
+ new PolynomialTimeDependence(std::vector<double>(1, 1.));
+ std::shared_ptr<AbstractTimeDependence> timePtr(time);
+ cav1.setAmplitudeModel(timePtr);
+ cav1.setPhaseModel(timePtr);
+ cav1.setFrequencyModel(timePtr);
+ cav1.setHeight(0.500);
+ cav1.setWidth(0.200);
+ cav1.setLength(1.000);
+ cav1.initialiseCoefficients();
+
+ cav2 = cav1;
+ PolynomialTimeDependence* phase2 =
+ new PolynomialTimeDependence(std::vector<double>(1, 0.));
+ PolynomialTimeDependence* freq2 =
+ new PolynomialTimeDependence(std::vector<double>(1, 1000.));
+ std::shared_ptr<AbstractTimeDependence> phase2Ptr(phase2);
+ std::shared_ptr<AbstractTimeDependence> freq2Ptr(freq2);
+ cav2.setPhaseModel(phase2Ptr);
+ cav2.setFrequencyModel(freq2Ptr);
+ }
+
+ VariableRFCavityFringeField cav1;
+ VariableRFCavityFringeField cav2;
+ std::shared_ptr<endfieldmodel::EndFieldModel> shared;
+};
+
+TEST_F(VariableRFCavityFringeFieldTest, TestConstructor) {
+ OpalTestUtilities::SilenceTest silencer;
+ VariableRFCavityFringeField cav("bob");
+ EXPECT_FLOAT_EQ(cav.getCavityCentre(), 0.);
+ endfieldmodel::EndFieldModel* null = NULL;
+ EXPECT_EQ(&(*(cav.getEndField())), null);
+}
+
+TEST_F(VariableRFCavityFringeFieldTest, TestSetGet) {
+ OpalTestUtilities::SilenceTest silencer;
+ EXPECT_FLOAT_EQ(cav1.getCavityCentre(), 0.5); // mm
+ EXPECT_FLOAT_EQ(cav1.getWidth(), 0.2); // metres
+ EXPECT_EQ(&(*(cav1.getEndField())), &(*shared));
+ long int count = shared.use_count();
+ EXPECT_EQ(count, 3);
+ PolynomialTimeDependence* time = NULL;
+ EXPECT_NE(&(*(cav1.getFrequencyModel())), time);
+}
+
+TEST_F(VariableRFCavityFringeFieldTest, TestClone) {
+ VariableRFCavityFringeField* cav2 =
+ dynamic_cast<VariableRFCavityFringeField*>(cav1.clone());
+ EXPECT_FLOAT_EQ(cav2->getCavityCentre(), 0.5);
+ EXPECT_EQ(&(*(cav2->getEndField())), &(*shared));
+ long int count = shared.use_count();
+ EXPECT_EQ(count, 4);
+}
+
+TEST_F(VariableRFCavityFringeFieldTest, TestAssignment) {
+ VariableRFCavityFringeField cav3 = cav1;
+ EXPECT_FLOAT_EQ(cav3.getCavityCentre(), 0.5);
+ EXPECT_EQ(&(*(cav3.getEndField())), &(*shared));
+ long int count = shared.use_count();
+ EXPECT_EQ(count, 4);
+}
+
+TEST_F(VariableRFCavityFringeFieldTest, TestCopyConstructor) {
+ VariableRFCavityFringeField cav4(cav1);
+ EXPECT_FLOAT_EQ(cav4.getCavityCentre(), 0.5);
+ EXPECT_EQ(&(*(cav4.getEndField())), &(*shared));
+ long int count = shared.use_count();
+ EXPECT_EQ(count, 4);
+}
+
+TEST_F(VariableRFCavityFringeFieldTest, TestApplyBoundingBox) {
+ Vector_t centroid(0., 0., 0.);
+ Vector_t B(0., 0., 0.);
+ Vector_t E(0., 0., 0.);
+ std::vector<Vector_t> rVectorFalse = {
+ Vector_t(-99.9, -249.9, 0.1),
+ Vector_t(-99.9, 249.9, 999.9),
+ };
+ for (size_t i = 0; i < rVectorFalse.size(); ++i) {
+ bool outOfBounds = cav1.apply(rVectorFalse[i], centroid, 0., B, E);
+ EXPECT_FALSE(outOfBounds);
+ }
+
+ std::vector<Vector_t> rVectorTrue = {
+ Vector_t(-100.1, 0., 500.),
+ Vector_t(+100.1, 0., 500.),
+ Vector_t(0., -250.1, 500.),
+ Vector_t(0., +250.1, 500.),
+ Vector_t(0., 0., -0.1),
+ Vector_t(0., 0., 1000.1),
+ };
+ for (size_t i = 0; i < rVectorTrue.size(); ++i) {
+ bool outOfBounds = cav1.apply(rVectorTrue[i], centroid, 0., B, E);
+ EXPECT_TRUE(outOfBounds) << rVectorTrue[i];
+ }
+}
+
+void testFieldLookup(VariableRFCavityFringeField& cav, Vector_t R, double t, Vector_t E, Vector_t B) {
+ Vector_t centroid(0., 0., 0.);
+ Vector_t Btest(0., 0., 0.);
+ Vector_t Etest(0., 0., 0.);
+ cav.apply(R, centroid, t, Etest, Btest);
+ for (size_t i = 0; i < 3; ++i) {
+ EXPECT_FLOAT_EQ(E[i], Etest[i]) << "\nR:" << R << " t: " << t
+ << "\nE expected: " << E << " " << " E meas: " << Etest
+ << "\nB expected: " << B << " " << " B meas: " << Btest << std::endl;
+ EXPECT_FLOAT_EQ(B[i], Btest[i]) << "\nR:" << R
+ << "\nE expected: " << E << " " << " E meas: " << Etest
+ << "\nB expected: " << B << " " << " B meas: " << Btest << std::endl;
+ }
+}
+
+void partial(VariableRFCavityFringeField& cav, Vector_t pos, double t, double delta, int var, Vector_t& dE, Vector_t& dB) {
+ bool verbose = false;
+ Vector_t centroid(0., 0., 0.);
+ Vector_t Bplus(0., 0., 0.);
+ Vector_t Bminus(0., 0., 0.);
+ Vector_t Eplus(0., 0., 0.);
+ Vector_t Eminus(0., 0., 0.);
+ Vector_t posPlus(pos), posMinus(pos);
+ double tPlus(t), tMinus(t);
+ if (var == 3) {
+ tMinus -= delta;
+ tPlus += delta;
+ } else if (var >= 0 or var < 3) {
+ posMinus[var] -= delta;
+ posPlus[var] += delta;
+ }
+ cav.apply(posPlus, centroid, tPlus, Eplus, Bplus);
+ cav.apply(posMinus, centroid, tMinus, Eminus, Bminus);
+ dE = (Eplus - Eminus)/2/delta;
+ dB = (Bplus - Bminus)/2/delta;
+ if (verbose) {
+ std::cerr << "Partial " << var << " " << delta << std::endl;
+ std::cerr << "R: " << posPlus << " " << posMinus << std::endl;
+ std::cerr << "t: " << tPlus << " " << tMinus << std::endl;
+ std::cerr << "B: " << Bplus << " " << Bminus << std::endl;
+ std::cerr << "E: " << Eplus << " " << Eminus << std::endl;
+ }
+}
+
+// curl B = 1/c^2 dE/dt
+Vector_t testMaxwell4(VariableRFCavityFringeField& cav, Vector_t pos, double t, double deltaPos, double deltaT) {
+ bool verbose = false;
+ Vector_t dummy;
+ Vector_t dBdx(0, 0, 0);
+ Vector_t dBdy(0, 0, 0);
+ Vector_t dBdz(0, 0, 0);
+ Vector_t dEdt(0, 0, 0);
+ partial(cav, pos, t, deltaPos, 0, dummy, dBdx);
+ partial(cav, pos, t, deltaPos, 1, dummy, dBdy);
+ partial(cav, pos, t, deltaPos, 2, dummy, dBdz);
+ partial(cav, pos, t, deltaT, 3, dEdt, dummy);
+ Vector_t curlB(
+ dBdy[2] - dBdz[1],
+ dBdz[0] - dBdx[2],
+ dBdx[1] - dBdy[0]
+ );
+ double c_l = Physics::c*1e-6; // 3e8 m/s = 300 mm/ns
+ Vector_t result = dEdt - curlB*1e-4*c_l*c_l;
+
+ if (verbose) {
+ std::cerr << "dBdx " << dBdx*1e-3 << std::endl;
+ std::cerr << "dBdy " << dBdy*1e-3 << std::endl;
+ std::cerr << "dBdz " << dBdz*1e-3 << std::endl;
+ std::cerr << "dEdt " << dEdt << std::endl;
+ std::cerr << "c^2 curlB " << curlB*1e-3*c_l*c_l << std::endl;
+ std::cerr << "dEdt-c^2 curlB " << result << std::endl << std::endl;
+ }
+ return result;
+}
+
+Vector_t testMaxwell3(VariableRFCavityFringeField& cav, Vector_t pos, double t, double deltaPos, double deltaT) {
+ bool verbose = false;
+ Vector_t dummy;
+ Vector_t dEdx(0, 0, 0);
+ Vector_t dEdy(0, 0, 0);
+ Vector_t dEdz(0, 0, 0);
+ Vector_t dBdt(0, 0, 0);
+ partial(cav, pos, t, deltaPos, 0, dEdx, dummy);
+ partial(cav, pos, t, deltaPos, 1, dEdy, dummy);
+ partial(cav, pos, t, deltaPos, 2, dEdz, dummy);
+ partial(cav, pos, t, deltaT, 3, dummy, dBdt);
+ Vector_t curlE(
+ dEdy[2] - dEdz[1],
+ dEdz[0] - dEdx[2],
+ dEdx[1] - dEdy[0]
+ );
+ Vector_t result = dBdt*1e-4 - curlE;
+
+ if (verbose) {
+ std::cerr << "maxwell3 at R: " << pos << " t: " << t << " with dx: "
+ << deltaPos << " dt: " << deltaT << std::endl;
+ std::cerr << " dEdx " << dEdx << std::endl;
+ std::cerr << " dEdy " << dEdy << std::endl;
+ std::cerr << " dEdz " << dEdz << std::endl;
+ std::cerr << " dBdt " << dBdt << std::endl;
+ std::cerr << " curlE " << curlE << std::endl;
+ std::cerr << " dBdt-curlE " << result << std::endl << std::endl;
+ }
+ return result;
+}
+
+
+std::vector<double> testMaxwell1and2(VariableRFCavityFringeField& cav, Vector_t pos, double t, double deltaPos) {
+ bool verbose = false;
+ Vector_t dummy;
+ Vector_t dEdx(0, 0, 0);
+ Vector_t dEdy(0, 0, 0);
+ Vector_t dEdz(0, 0, 0);
+ Vector_t dBdx(0, 0, 0);
+ Vector_t dBdy(0, 0, 0);
+ Vector_t dBdz(0, 0, 0);
+ partial(cav, pos, t, deltaPos, 0, dEdx, dBdx);
+ partial(cav, pos, t, deltaPos, 1, dEdy, dBdy);
+ partial(cav, pos, t, deltaPos, 2, dEdz, dBdz);
+
+ double divE = dEdx[0] + dEdy[1] + dEdz[2];
+ double divB = dBdx[0] + dBdy[1] + dBdz[2];
+
+ if (verbose) {
+ std::cerr << "maxwell1+2 at R: " << pos << " t: " << t << " with dx: "
+ << deltaPos << std::endl;
+ std::cerr << " dEidi " << dEdx[0] << " " << dEdy[1] << " " << dEdz[2] << std::endl;
+ std::cerr << " dBidi " << dBdx[0] << " " << dBdy[1] << " " << dBdz[2] << std::endl;
+ std::cerr << " divE " << divE << std::endl;
+ std::cerr << " divB " << divB << std::endl;
+ }
+ std::vector<double> result = {divE, divB};
+ return result;
+}
+
+TEST_F(VariableRFCavityFringeFieldTest, TestMaxwell) {
+ //double pi = Physics::pi;
+ Vector_t centroid(0., 0., 0.);
+ Vector_t B(0., 0., 0.);
+ Vector_t E(0., 0., 0.);
+ Vector_t R(0., 0., 500.);
+ double t = 0.;
+ std::cerr << "\nOff midplane, 45 degree phase, in fringe field" << std::endl;
+ std::cerr << "order B E max1 max2 maxwell3 maxwell4" << std::endl;
+ R = Vector_t(0., 1., 750.);
+ t = 0.125;
+ for (size_t i = 0; i < 5; ++i) {
+ cav2.setMaxOrder(i);
+ cav2.apply(R, centroid, t, E, B);
+ Vector_t result1 = testMaxwell3(cav2, R, t, 0.01, 0.0001);
+ Vector_t result2 = testMaxwell4(cav2, R, t, 0.01, 0.0001);
+ std::vector<double> div = testMaxwell1and2(cav2, R, t, 0.01);
+ std::cerr << i << " ** " << B << " " << E << " " << div[0] << " " << div[1] << " " << result1 << " " << result2 << std::endl;
+
+
+ }
+}