ParallelCyclotronTracker.cpp

//
// Class ParallelCyclotronTracker
//   Tracker for OPAL-Cycl
//
// Copyright (c) 2007 - 2014, Jianjun Yang, Andreas Adelmann and Matthias Toggweiler,
//                            Paul Scherrer Institut, Villigen PSI, Switzerland
// Copyright (c) 2014,        Daniel Winklehner, MIT, Cambridge, MA, USA
// Copyright (c) 2012 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Beam dynamics in high intensity cyclotrons including neighboring bunch effects"
// and the paper
// "Beam dynamics in high intensity cyclotrons including neighboring bunch effects:
// Model, implementation, and application"
// (https://journals.aps.org/prab/pdf/10.1103/PhysRevSTAB.13.064201)
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Algorithms/ParallelCyclotronTracker.h"

#include <fstream>
#include <iostream>
#include <limits>
#include <vector>
#include <numeric>
#include <cmath>

#include "AbstractObjects/Element.h"
#include "AbstractObjects/OpalData.h"

#include "AbsBeamline/CCollimator.h"
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Cyclotron.h"
#include "AbsBeamline/Degrader.h"
#include "AbsBeamline/Drift.h"
#include "AbsBeamline/Marker.h"
#include "AbsBeamline/Monitor.h"
#include "AbsBeamline/Multipole.h"
#include "AbsBeamline/MultipoleT.h"
#include "AbsBeamline/MultipoleTBase.h"
#include "AbsBeamline/MultipoleTStraight.h"
#include "AbsBeamline/MultipoleTCurvedConstRadius.h"
#include "AbsBeamline/MultipoleTCurvedVarRadius.h"
#include "AbsBeamline/Offset.h"
#include "AbsBeamline/PluginElement.h"
#include "AbsBeamline/Probe.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/RFCavity.h"
#include "AbsBeamline/Ring.h"
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/SBend3D.h"
#include "AbsBeamline/ScalingFFAMagnet.h"
#include "AbsBeamline/Septum.h"
#include "AbsBeamline/Solenoid.h"
#include "AbsBeamline/Stripper.h"
#include "AbsBeamline/Vacuum.h"
#include "AbsBeamline/VariableRFCavity.h"
#include "AbsBeamline/VariableRFCavityFringeField.h"
#include "AbsBeamline/VerticalFFAMagnet.h"

#include "Algorithms/Ctunes.h"
#include "Algorithms/PolynomialTimeDependence.h"

#include "Beamlines/Beamline.h"
#include "Beamlines/FlaggedBeamline.h"

#include "Elements/OpalBeamline.h"

#include "Physics/Physics.h"

#include "Utilities/OpalException.h"
#include "Utilities/Options.h"

#include "BasicActions/DumpFields.h"
#include "BasicActions/DumpEMFields.h"

#include "Structure/BoundaryGeometry.h"
#include "Structure/DataSink.h"
#include "Structure/LossDataSink.h"


constexpr double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns

Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0);
Vector_t const ParallelCyclotronTracker::yaxis = Vector_t(0.0, 1.0, 0.0);
Vector_t const ParallelCyclotronTracker::zaxis = Vector_t(0.0, 0.0, 1.0);

extern Inform *gmsg;

/**
 * Constructor ParallelCyclotronTracker
 *
 * @param beamline
 * @param bunch
 * @param ds
 * @param reference
 * @param revBeam
 * @param revTrack
 * @param maxSTEPS
 * @param timeIntegrator
 */
ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
                                                   PartBunchBase<double, 3> *bunch,
                                                   DataSink &ds,
                                                   const PartData &reference,
                                                   bool revBeam, bool revTrack,
                                                   int maxSTEPS, int timeIntegrator,
                                                   const int& numBunch,
                                                   const double& mbEta,
                                                   const double& mbPara,
                                                   const std::string& mbMode,
                                                   const std::string& mbBinning)
    : Tracker(beamline, bunch, reference, revBeam, revTrack)
    , bgf_m(nullptr)
    , maxSteps_m(maxSTEPS)
    , lastDumpedStep_m(0)
    , myNode_m(Ippl::myNode())
    , initialLocalNum_m(bunch->getLocalNum())
    , initialTotalNum_m(bunch->getTotalNum())
    , opalRing_m(nullptr)
    , itsStepper_mp(nullptr)
{
    itsBeamline = dynamic_cast<Beamline *>(beamline.clone());
    itsDataSink = &ds;

    if ( numBunch > 1 ) {
        mbHandler_m = std::unique_ptr<MultiBunchHandler>(
            new MultiBunchHandler(bunch, numBunch, mbEta,
                                  mbPara, mbMode, mbBinning)
        );
    }

    IntegrationTimer_m = IpplTimings::getTimer("Integration");
    TransformTimer_m   = IpplTimings::getTimer("Frametransform");
    DumpTimer_m        = IpplTimings::getTimer("Dump");
    BinRepartTimer_m   = IpplTimings::getTimer("Binaryrepart");
    PluginElemTimer_m  = IpplTimings::getTimer("PluginElements");
    DelParticleTimer_m = IpplTimings::getTimer("DeleteParticles");

    // FIXME Change track command
    if ( initialTotalNum_m == 1 ) {
        mode_m = MODE::SINGLE;
    } else if ( initialTotalNum_m == 2 ) {
        mode_m = MODE::SEO;
    } else if ( initialTotalNum_m > 2 ) {
        mode_m = MODE::BUNCH;
    } else
        mode_m = MODE::UNDEFINED;

    if ( timeIntegrator == 0 ) {
        stepper_m = stepper::INTEGRATOR::RK4;
    } else if ( timeIntegrator == 1) {
        stepper_m = stepper::INTEGRATOR::LF2;
    } else if ( timeIntegrator == 2) {
        stepper_m = stepper::INTEGRATOR::MTS;
    } else
        stepper_m = stepper::INTEGRATOR::UNDEFINED;
}

/**
 * Destructor ParallelCyclotronTracker
 *
 */
ParallelCyclotronTracker::~ParallelCyclotronTracker() {
    if(bgf_m)
        lossDs_m->save();
    for(Component* component : myElements) {
        delete(component);
    }
    for(auto fd : FieldDimensions) {
        delete(fd);
    }
    delete itsBeamline;
    // delete opalRing_m;
}


void ParallelCyclotronTracker::bgf_main_collision_test() {
    if(!bgf_m) return;

    Inform msg("bgf_main_collision_test ");

    /**
     *Here we check if a particle is outside the domain, flag it for deletion
     */

    Vector_t intecoords = 0.0;

    // This has to match the dT in the rk4 pusher
    double dtime = itsBunch_m->getdT() * getHarmonicNumber();

    int triId = 0;
    for(size_t i = 0; i < itsBunch_m->getLocalNum(); i++) {
        int res = bgf_m->partInside(itsBunch_m->R[i], itsBunch_m->P[i],
                                    dtime, intecoords, triId);
        if(res >= 0) {
            lossDs_m->addParticle(OpalParticle(itsBunch_m->ID[i],
                                               itsBunch_m->R[i], itsBunch_m->P[i],
                                               itsBunch_m->getT(),
                                               itsBunch_m->Q[i], itsBunch_m->M[i]),
                                  std::make_pair(turnnumber_m, itsBunch_m->bunchNum[i]));
            itsBunch_m->Bin[i] = -1;
            Inform gmsgALL("OPAL", INFORM_ALL_NODES);
            gmsgALL << level4 << "* Particle " << itsBunch_m->ID[i]
                    << " lost on boundary geometry" << endl;
        }
    }
}

// only used for dumping into stat file
void ParallelCyclotronTracker::dumpAngle(const double& theta,
                                         double& prevAzimuth,
                                         double& azimuth,
                                         const short& bunchNr)
{
    if ( prevAzimuth < 0.0 ) { // only at first occurrence
        double plus = 0.0;
        if ( OpalData::getInstance()->inRestartRun() ) {
            plus = 360.0 * (turnnumber_m - bunchNr);
        }
        azimuth = theta + plus;
    } else {
        double dtheta = theta - prevAzimuth;
        if ( dtheta < 0 ) {
            dtheta += 360.0;
        }
        if ( dtheta > 180 ) { // rotating clockwise, reduce angle
            dtheta -= 360;
        }
        azimuth += dtheta;
    }
    prevAzimuth = theta;
}


double ParallelCyclotronTracker::computeRadius(const Vector_t& meanR) const {
    // New OPAL 2.0: m --> mm
    return 1000.0 * std::sqrt(meanR(0) * meanR(0) + meanR(1) * meanR(1));
}


void ParallelCyclotronTracker::computePathLengthUpdate(std::vector<double>& dl,
                                                       const double& dt)
{
    // the last element in dotP is the dot-product over all particles
    std::vector<double> dotP(dl.size());
    if ( Options::psDumpFrame == Options::BUNCH_MEAN || isMultiBunch()) {

        for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
            dotP[itsBunch_m->bunchNum[i]] += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
        }

        allreduce(dotP.data(), dotP.size(), std::plus<double>());

        // dot-product over all particles
        double sum = std::accumulate(dotP.begin(), dotP.end() - 1, 0);
        dotP.back() = sum / double(itsBunch_m->getTotalNum());

        // bunch specific --> multi-bunches only
        for (short b = 0; b < (short)dotP.size() - 1; ++b) {
            dotP[b] /= double(itsBunch_m->getTotalNumPerBunch(b));
        }

    } else if ( itsBunch_m->getLocalNum() == 0 ) {
        // here we are in Options::GLOBAL mode
        dotP[0] = 0.0;
    } else {
        // here we are in Options::GLOBAL mode
        dotP[0] = dot(itsBunch_m->P[0], itsBunch_m->P[0]);
    }

    for (size_t i = 0; i < dotP.size(); ++i) {
        double const gamma = std::sqrt(1.0 + dotP[i]);
        double const beta  = std::sqrt(dotP[i]) / gamma;
        dl[i] = c_mmtns * dt * 1.0e-3 * beta; // unit: m
    }
}



/**
 *
 *
 * @param fn Base file name
 */
void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {

    for (unsigned int i=0; i<numFiles; i++) {
        std::string SfileName2 = SfileName;
        if (i<=2) {
            SfileName2 += std::string("-Angle" + std::to_string(int(i)) + ".dat");
        }
        else {
            // for single Particle Mode, output after each turn, to define matched initial phase ellipse.
            SfileName2 += std::string("-afterEachTurn.dat");
        }

        outfTheta_m.emplace_back(new std::ofstream(SfileName2.c_str()));
        outfTheta_m.back()->precision(8);
        outfTheta_m.back()->setf(std::ios::scientific, std::ios::floatfield);
        *outfTheta_m.back() << "# r [mm]        beta_r*gamma       "
                            << "theta [deg]     beta_theta*gamma        "
                            << "z [mm]          beta_z*gamma" << std::endl;
    }
}

/**
 * Close all files related to
 * special output in the Cyclotron
 * mode.
 */
void ParallelCyclotronTracker::closeFiles() {
    for (auto & file : outfTheta_m) {
        file->close();
    }
}


/**
 *
 *
 * @param cycl
 */
void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {

    *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;

    cycl_m = dynamic_cast<Cyclotron *>(cycl.clone());
    myElements.push_back(cycl_m);

    // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
    // useful information
    spiral_flag = cycl_m->getSpiralFlag();

    if (spiral_flag) {
        *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
        *gmsg         << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
        *gmsg         << "*     (Use BANDRF type cyclotron and use RFMAPFN to load both magnetic" << endl;
        *gmsg         << "*     and electric fields, setting SUPERPOSE to an array of TRUE values.)" << endl;
        *gmsg         << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
        *gmsg         << "*     FFT does not give the correct results (boundary conditions are missing)." << endl;
        *gmsg         << "* 3.) The whole geometry will be meshed and used for the fieldsolver." << endl;
        *gmsg         << "*     There will be no transformations of the bunch into a local frame und consequently," << endl;
        *gmsg         << "*     the problem will be treated non-relativistically!" << endl;
        *gmsg         << "*     (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" << endl;
        *gmsg << endl << "* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" << endl;
        *gmsg         << "* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." << endl;
        if (isMultiBunch()) {
            mbHandler_m = nullptr;
        }
    }

    // Fresh run (no restart):
    if (!OpalData::getInstance()->inRestartRun()) {

        // Get reference values from cyclotron element
        // For now, these are still stored in mm. should be the only ones. -DW
        referenceR     = cycl_m->getRinit();
        referenceTheta = cycl_m->getPHIinit();
        referenceZ     = cycl_m->getZinit();
        referencePr    = cycl_m->getPRinit();
        referencePz    = cycl_m->getPZinit();

        if (referenceTheta <= -180.0 || referenceTheta > 180.0) {
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
        }

        referencePtot =  itsReference.getGamma() * itsReference.getBeta();

        // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
        float insqrt = referencePtot * referencePtot - \
            referencePr * referencePr - referencePz * referencePz;

        if (insqrt < 0) {

            if (insqrt > -1.0e-10) {
                referencePt = 0.0;
            } else {
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "Pt imaginary!");
            }

        } else {
            referencePt = std::sqrt(insqrt);
        }

        if (referencePtot < 0.0) {
            referencePt *= -1.0;
        }
        // End calculate referencePt

        // Restart a run:
    } else {

        // If the user wants to save the restarted run in local frame,
        // make sure the previous h5 file was local too
        if (Options::psDumpFrame != Options::GLOBAL) {
            if (!previousH5Local) {
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a local restart from a global h5 file!");
            }
            // Else, if the user wants to save the restarted run in global frame,
            // make sure the previous h5 file was global too
        } else {
            if (previousH5Local) {
                throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                    "You are trying a global restart from a local h5 file!");
            }
        }

        // Adjust some of the reference variables from the h5 file
        referencePhi *= Physics::deg2rad;
        referencePsi *= Physics::deg2rad;
        referencePtot = bega;
        if (referenceTheta <= -180.0 || referenceTheta > 180.0) {
            throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
                                "PHIINIT is out of [-180, 180)!");
        }
    }

    sinRefTheta_m = std::sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = std::cos(referenceTheta * Physics::deg2rad);

    *gmsg << endl;
    *gmsg << "* Bunch global starting position:" << endl;
    *gmsg << "* RINIT = " << referenceR  << " [mm]" << endl;
    *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl;
    *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl;
    *gmsg << endl;
    *gmsg << "* Bunch global starting momenta:" << endl;
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
    *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl;
    *gmsg << "* Reference azimuthal momentum (Pt) = " << referencePt * 1000.0 << " [MCU]" << endl;
    *gmsg << "* Reference radial momentum (Pr) = " << referencePr * 1000.0 << " [MCU]" << endl;
    *gmsg << "* Reference axial momentum (Pz) = " << referencePz * 1000.0 << " [MCU]" << endl;
    *gmsg << endl;

    double sym = cycl_m->getSymmetry();
    *gmsg << "* " << sym << "-fold field symmetry " << endl;

    // ckr: this just returned the default value as defined in Component.h
    // double rff = cycl_m->getRfFrequ(0);
    // *gmsg << "* Rf frequency= " << rff << " [MHz]" << endl;

    std::string fmfn = cycl_m->getFieldMapFN();
    *gmsg << "* Field map file name = " << fmfn << " " << endl;

    std::string type = cycl_m->getCyclotronType();
    *gmsg << "* Type of cyclotron = " << type << " " << endl;

    double rmin = cycl_m->getMinR();
    double rmax = cycl_m->getMaxR();
    *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;

    double zmin = cycl_m->getMinZ();
    double zmax = cycl_m->getMaxZ();
    *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [m]"<< endl;

    double h = cycl_m->getCyclHarm();
    *gmsg << "* Number of trimcoils = " << cycl_m->getNumberOfTrimcoils() << endl;
    *gmsg << "* Harmonic number h = " << h << " " << endl;

    if (type == std::string("BANDRF")) {
        double escale = cycl_m->getEScale(0);
        *gmsg << "* RF field scale factor = " << escale << endl;
        double rfphi= cycl_m->getRfPhi(0);
        *gmsg << "* RF inital phase = " << rfphi * Physics::rad2deg << " [deg]" << endl;
        bool superpose = cycl_m->getSuperpose(0);
        *gmsg << std::boolalpha << "* Superpose electric field maps -> " << superpose << endl;
    }

    // Read in cyclotron field maps
    cycl_m->initialise(itsBunch_m, cycl_m->getBScale());

    double BcParameter[8] = {};

    BcParameter[0] = 0.001 * cycl_m->getRmin();
    BcParameter[1] = 0.001 * cycl_m->getRmax();

    // Store inner radius and outer radius of cyclotron field map in the list
    buildupFieldList(BcParameter, ElementBase::CYCLOTRON, cycl_m);
}

/**
 *
 *
 * @param coll
 */
void ParallelCyclotronTracker::visitCCollimator(const CCollimator &coll) {

    *gmsg << "* ------------------------------ Collimator ------------------------------" << endl;

    CCollimator* elptr = dynamic_cast<CCollimator *>(coll.clone());
    myElements.push_back(elptr);

    double xstart = elptr->getXStart();
    *gmsg << "* Xstart  = " << xstart << " [m]" << endl;

    double xend = elptr->getXEnd();
    *gmsg << "* Xend    = " << xend << " [m]" << endl;

    double ystart = elptr->getYStart();
    *gmsg << "* Ystart  = " << ystart << " [m]" << endl;

    double yend = elptr->getYEnd();
    *gmsg << "* Yend    = " << yend << " [m]" << endl;

    double zstart = elptr->getZStart();
    *gmsg << "* Zstart  = " << zstart << " [m]" << endl;

    double zend = elptr->getZEnd();
    *gmsg << "* Zend    = " << zend << " [m]" << endl;

    double width = elptr->getWidth();
    *gmsg << "* Width   = " << width << " [m]" << endl;

    elptr->initialise(itsBunch_m);

    double BcParameter[8] = {};

    BcParameter[0] = xstart;
    BcParameter[1] = xend;
    BcParameter[2] = ystart;
    BcParameter[3] = yend;
    BcParameter[4] = width;

    buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr);
}

/**
 *
 *
 * @param corr
 */
void ParallelCyclotronTracker::visitCorrector(const Corrector &corr) {
    *gmsg << "In Corrector; L= " << corr.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Corrector *>(corr.clone()));
}

/**
 *
 *
 * @param degrader
 */
void ParallelCyclotronTracker::visitDegrader(const Degrader &deg) {
    *gmsg << "In Degrader; L= " << deg.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Degrader *>(deg.clone()));

}

/**
 *
 *
 * @param drift
 */
void ParallelCyclotronTracker::visitDrift(const Drift &drift) {
    *gmsg << "In drift L= " << drift.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Drift *>(drift.clone()));
}

/**
 *
 *
 *  @param
 */
void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {

}

/**
 *
 *
 * @param off
 */
void ParallelCyclotronTracker::visitOffset(const Offset& off) {
    if (opalRing_m == NULL)
        throw OpalException(
                            "ParallelCylcotronTracker::visitOffset",
                            "Attempt to place an offset when Ring not defined");
    Offset* offNonConst = const_cast<Offset*>(&off);
    offNonConst->updateGeometry(opalRing_m->getNextPosition(),
                                opalRing_m->getNextNormal());
    opalRing_m->appendElement(off);
}

/**
 *
 *
 * @param marker
 */
void ParallelCyclotronTracker::visitMarker(const Marker &marker) {
    //   *gmsg << "In Marker; L= " << marker.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Marker *>(marker.clone()));
    // Do nothing.
}

/**
 *
 *
 * @param corr
 */
void ParallelCyclotronTracker::visitMonitor(const Monitor &corr) {
    //   *gmsg << "In Monitor; L= " << corr.getElementLength() << endl;
    myElements.push_back(dynamic_cast<Monitor *>(corr.clone()));
    //   applyDrift(flip_s * corr.getElementLength());
}

/**
 *
 *
 * @param mult
 */
void ParallelCyclotronTracker::visitMultipole(const Multipole &mult) {
    *gmsg << "In Multipole; L= " << mult.getElementLength() << " however the element is missing " << endl;
    myElements.push_back(dynamic_cast<Multipole *>(mult.clone()));
}

/**
 *
 *
 * @param multT
 */
void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) {
    *gmsg << "Adding MultipoleT" << endl;
    if (opalRing_m != NULL) {
        opalRing_m->appendElement(multT);
    } else {
        throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
                            "Need to define a RINGDEFINITION to use MultipoleT element");
    }
    myElements.push_back(dynamic_cast<MultipoleT *>(multT.clone()));
}

/**
 *
 *
 * @param multTstraight
 */
void ParallelCyclotronTracker::visitMultipoleTStraight(const MultipoleTStraight &multTstraight) {
    *gmsg << "Adding MultipoleTStraight" << endl;
    if (opalRing_m != NULL) {
        opalRing_m->appendElement(multTstraight);
    } else {
        throw OpalException("ParallelCyclotronTracker::visitMultipoleTStraight",
                            "Need to define a RINGDEFINITION to use MultipoleTStraight element");
    }
    myElements.push_back(dynamic_cast<MultipoleTStraight *>(multTstraight.clone()));
}

/**
 *
 *
 * @param multTccurv
 */
void ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius(const MultipoleTCurvedConstRadius &multTccurv) {
    *gmsg << "Adding MultipoleTCurvedConstRadius" << endl;
    if (opalRing_m != NULL) {
        opalRing_m->appendElement(multTccurv);
    } else {
        throw OpalException("ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius",
                            "Need to define a RINGDEFINITION to use MultipoleTCurvedConstRadius element");
    }
    myElements.push_back(dynamic_cast<MultipoleTCurvedConstRadius *>(multTccurv.clone()));
}

/**
 *
 *
 * @param multTvcurv
 */
void ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius(const MultipoleTCurvedVarRadius &multTvcurv) {
    *gmsg << "Adding MultipoleTCurvedVarRadius" << endl;
    if (opalRing_m != NULL) {
        opalRing_m->appendElement(multTvcurv);
    } else {
        throw OpalException("ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius",
                            "Need to define a RINGDEFINITION to use MultipoleTCurvedVarRadius element");
    }
    myElements.push_back(dynamic_cast<MultipoleTCurvedVarRadius *>(multTvcurv.clone()));
}

/**
 *
 *
 * @param prob
 */
void ParallelCyclotronTracker::visitProbe(const Probe &prob) {
    *gmsg << "* ------------------------------  Probe ------------------------------" << endl;
    Probe *elptr = dynamic_cast<Probe *>(prob.clone());
    myElements.push_back(elptr);

    *gmsg << "* Name    = " << elptr->getName() << endl;

    double xstart = elptr->getXStart();
    *gmsg << "* XStart  = " << xstart << " [m]" << endl;

    double xend = elptr->getXEnd();
    *gmsg << "* XEnd    = " << xend << " [m]" << endl;

    double ystart = elptr->getYStart();
    *gmsg << "* YStart  = " << ystart << " [m]" << endl;

    double yend = elptr->getYEnd();
    *gmsg << "* YEnd    = " << yend << " [m]" << endl;

    // initialise, do nothing
    elptr->initialise(itsBunch_m);

    double BcParameter[8] = {};

    BcParameter[0] = xstart;
    BcParameter[1] = xend;
    BcParameter[2] = ystart;
    BcParameter[3] = yend;
    BcParameter[4] = 1 ; // width

    // store probe parameters in the list
    buildupFieldList(BcParameter, ElementBase::PROBE, elptr);
}

/**
 *
 *
 * @param bend
 */
void ParallelCyclotronTracker::visitRBend(const RBend &bend) {
    *gmsg << "In RBend; L= " << bend.getElementLength() << " however the element is missing " << endl;
    myElements.push_back(dynamic_cast<RBend *>(bend.clone()));
}

/**
 *
 *
 * @param as
 */
void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) {

    *gmsg << "* ------------------------------ RFCavity ------------------------------" << endl;

    RFCavity *elptr = dynamic_cast<RFCavity *>(as.clone());
    myElements.push_back(elptr);

    if ( elptr->getComponentType() != "SINGLEGAP" ) {
        *gmsg << (elptr->getComponentType()) << endl;
        throw OpalException("ParallelCyclotronTracker::visitRFCavity",
                            "The ParallelCyclotronTracker can only play with cyclotron type RF system currently ...");
    }

    double rmin = elptr->getRmin();
    *gmsg << "* Minimal radius of cavity= " << rmin << " [mm]" << endl;

    double rmax = elptr->getRmax();
    *gmsg << "* Maximal radius of cavity= " << rmax << " [mm]" << endl;

    double rff = elptr->getCycFrequency();
    *gmsg << "* RF frequency (2*pi*f)= " << rff << " [rad/s]" << endl;

    std::string fmfn = elptr->getFieldMapFN();
    *gmsg << "* RF Field map file name= " << fmfn << endl;

    double angle = elptr->getAzimuth();
    *gmsg << "* Cavity azimuth position= " << angle << " [deg] " << endl;

    double gap = elptr->getGapWidth();
    *gmsg << "* Cavity gap width= " << gap << " [mm] " << endl;

    double pdis = elptr->getPerpenDistance();
    *gmsg << "* Cavity Shift distance= " << pdis << " [mm] " << endl;

    double phi0 = elptr->getPhi0();
    *gmsg << "* Initial RF phase (t=0)= " << phi0 << " [deg] " << endl;

    /*
      Setup time dependence and in case of no
      timedependence use a polynom with  a_0 = 1 and a_k = 0, k = 1,2,3.
    */

    std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
    std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;

    dvector_t  unityVec;
    unityVec.push_back(1.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);
    unityVec.push_back(0.);

    if (elptr->getFrequencyModelName() != "") {
        freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName());
        *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl;
    }
    else
        freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));

    if (elptr->getAmplitudeModelName() != "") {
        ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName());
        *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl;
    }
    else
        ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));

    if (elptr->getPhaseModelName() != "") {
        phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName());
        *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl;
    }
    else
        phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));

    // read cavity voltage profile data from file.
    elptr->initialise(itsBunch_m, freq_atd, ampl_atd, phase_atd);

    double BcParameter[8] = {};

    BcParameter[0] = 0.001 * rmin;
    BcParameter[1] = 0.001 * rmax;
    BcParameter[2] = 0.001 * pdis;
    BcParameter[3] = angle;

    buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr);
}

/**
 *
 * @param ring
 */
void ParallelCyclotronTracker::visitRing(const Ring &ring) {

    *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl;

    delete opalRing_m;

    opalRing_m = dynamic_cast<Ring*>(ring.clone());

    myElements.push_back(opalRing_m);

    opalRing_m->initialise(itsBunch_m);

    referenceR = opalRing_m->getBeamRInit();
    referencePr = opalRing_m->getBeamPRInit();
    referenceTheta = opalRing_m->getBeamPhiInit();

    if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
        throw OpalException("Error in ParallelCyclotronTracker::visitRing",
                            "PHIINIT is out of [-180, 180)!");
    }

    referenceZ = 0.0;
    referencePz = 0.0;

    referencePtot = itsReference.getGamma() * itsReference.getBeta();
    referencePt = std::sqrt(referencePtot * referencePtot - referencePr * referencePr);

    if(referencePtot < 0.0)
        referencePt *= -1.0;

    sinRefTheta_m = std::sin(referenceTheta * Physics::deg2rad);
    cosRefTheta_m = std::cos(referenceTheta * Physics::deg2rad);

    double BcParameter[8] = {}; // zero initialise array

    buildupFieldList(BcParameter, ElementBase::RING, opalRing_m);

    // Finally print some diagnostic
    *gmsg << "* Initial beam radius = " << referenceR << " [mm] " << endl;
    *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
    *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
    *gmsg << "* Total reference momentum   = " << referencePtot * 1000.0
          << " [MCU]" << endl;
    *gmsg << "* Reference azimuthal momentum  = " << referencePt * 1000.0
          << " [MCU]" << endl;
    *gmsg << "* Reference radial momentum     = " << referencePr * 1000.0
          << " [MCU]" << endl;
    *gmsg << "* " << opalRing_m->getSymmetry() << " fold field symmetry "
          << endl;
    *gmsg << "* Harmonic number h= " << opalRing_m->getHarmonicNumber() << " "
          << endl;
}