VariableRFCavityFringeField - passes tests, needs a bit of clean up

.gitignore

@@ -3,6 +3,7 @@
 *.o
 *.a
 *.aux
+optimizer/Tests/*.exe
 build
 CMakeCache.txt
 CMakeFiles

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");
+}
+
 /**
  *
  *

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();

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;
 }
 

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;
 

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

src/Classic/AbsBeamline/EndFieldModel/EndFieldModel.h

@@ -32,7 +32,7 @@
 #include <vector>
 
 namespace endfieldmodel {
-
+  
 class EndFieldModel {
  public:
   virtual ~EndFieldModel() {;}

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;
+}
 }
 

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;
 

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);

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",

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 {

src/Classic/AbsBeamline/VariableRFCavityFringeField.cpp

+/*
+ *  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);