Opal maps

This MR will implement Map tracking for:

• Quadrupole
• Dipole
• Drifts
• Bends

TODO list:

• Code implementation (@ganz_p): ThickTracker, Hamiltonian, MapAnalyser
• Doxygen documentation (@ganz_p)
• Documentation (@ganz_p)
• Code cleaning (@frey_m)
• Regression tests (@snuverink_j)
• Examples (@ganz_p)
• Jupyter notebook (@ganz_p)
• Use beam line order

Remaining issues:

• Output dumps in local coordinates only
• Dipole until 1st order
• Tune calculation not stable
• Space charge calculations needs to be verified
• Fringe fields for Bends

src/Algorithms/ThickTracker.h

@@ -2,10 +2,6 @@
 #define OPAL_ThickTracker_HH
 
 // ------------------------------------------------------------------------
-// $RCSfile: ThickTracker.h,v $
-// ------------------------------------------------------------------------
-// $Revision: 1.1.2.1 $
-// ------------------------------------------------------------------------
 // Copyright: see Copyright.readme
 // ------------------------------------------------------------------------
 //
@@ -13,25 +9,23 @@
 //
 // ------------------------------------------------------------------------
 //
-// $Date: 2004/11/12 20:10:11 $
-// $Author: adelmann $
+// $Author: ganz_p $
 //
 // ------------------------------------------------------------------------
 
 #include "Algorithms/Tracker.h"
-#include "Steppers/BorisPusher.h"
 #include "Structure/DataSink.h"
 
-#include "Classic/FixedAlgebra/FTps.h"
-#include "Classic/FixedAlgebra/FVps.h"
+#include "Hamiltonian.h"
 
+#include "MapAnalyser.h"
 
-#include "Algorithms/OrbitThreader.h"
 #include "Algorithms/IndexMap.h"
-#include "AbsBeamline/AlignWrapper.h"
 #include "AbsBeamline/BeamBeam.h"
+
 #include "AbsBeamline/CCollimator.h"
 #include "AbsBeamline/Corrector.h"
+#include "AbsBeamline/CyclotronValley.h"
 #include "AbsBeamline/Diagnostic.h"
 #include "AbsBeamline/Degrader.h"
 #include "AbsBeamline/Drift.h"
@@ -41,33 +35,32 @@
 #include "AbsBeamline/Marker.h"
 #include "AbsBeamline/Monitor.h"
 #include "AbsBeamline/Multipole.h"
+#include "AbsBeamline/ParallelPlate.h"
 #include "AbsBeamline/Probe.h"
-#include "AbsBeamline/RBend.h"
-#include "AbsBeamline/RBend3D.h"
 #include "AbsBeamline/RFCavity.h"
-#include "AbsBeamline/TravelingWave.h"
 #include "AbsBeamline/RFQuadrupole.h"
+#include "AbsBeamline/RBend.h"
+#include "AbsBeamline/RBend3D.h"
 #include "AbsBeamline/SBend.h"
 #include "AbsBeamline/Separator.h"
 #include "AbsBeamline/Septum.h"
 #include "AbsBeamline/Solenoid.h"
-#include "AbsBeamline/ParallelPlate.h"
-#include "AbsBeamline/CyclotronValley.h"
+#include "AbsBeamline/TravelingWave.h"
 
 #include "Elements/OpalBeamline.h"
 
-#include <cmath>
+#include "Structure/Beam.h"
 
-#define PSdim 6
+//#include <array>
+#include <cmath>
 
-#define PHIL_WRITE 1
+#include <tuple>
 
 class BMultipoleField;
 
 template <class T, unsigned Dim>
 class PartBunchBase;
 
-class PlanarArcGeometry;
 
 
 // Class ThickTracker
@@ -100,9 +93,9 @@ class PlanarArcGeometry;
 // [/tab][p]
 // Approximations used:
 // [ul]
-// [li] blah
-// [li] blah
-// [li] blah
+// [li]
+// [li]
+// [li]
 // [/ul]
 //
 // On going through an element, we use the following steps:
@@ -111,6 +104,15 @@ class PlanarArcGeometry;
 class ThickTracker: public Tracker {
 
 public:
+    typedef Hamiltonian::series_t                       series_t;
+    typedef MapAnalyser::fMatrix_t                      fMatrix_t;
+    typedef MapAnalyser::cfMatrix_t                     cfMatrix_t;
+    typedef FVps<double, 6>                             map_t;
+    typedef FVector<double, 6>                          particle_t;
+    typedef std::tuple<series_t, std::size_t, double>   tuple_t;
+    typedef std::list<tuple_t>                          beamline_t;
+
+public:
 
     /// Constructor.
     //  The beam line to be tracked is "bl".
@@ -127,16 +129,18 @@ public:
     //  The particle bunch tracked is taken from [b]bunch[/b].
     //  If [b]revBeam[/b] is true, the beam runs from s = C to s = 0.
     //  If [b]revTrack[/b] is true, we track against the beam.
-    explicit ThickTracker(const Beamline &bl, PartBunchBase<double, 3> *bunch,
-			  DataSink &ds,
+    explicit ThickTracker(const Beamline &bl,
+                          PartBunchBase<double, 3> *bunch,
+                          Beam &beam,
+                          DataSink &ds,
                           const PartData &data,
-			  bool revBeam, bool revTrack,
-			  const std::vector<unsigned long long> &maxSTEPS,
-			  double zstart,
-			  const std::vector<double> &zstop,
-			  const std::vector<double> &dt,
+                          bool revBeam, bool revTrack,
+                          const std::vector<unsigned long long> &maxSTEPS,
+                          double zstart,
+                          const std::vector<double> &zstop,
+                          const std::vector<double> &dt,
                           const int& truncOrder);
-    
+
     virtual ~ThickTracker();
 
 
@@ -217,73 +221,144 @@ public:
     //  overwrite the execute-methode from DefaultVisitor
     virtual void execute();
 
-    void updateReferenceParticle(const BorisPusher &pusher);
-    void updateRFElement(std::string elName, double maxPhase);
     void prepareSections();
-    void saveCavityPhases();
-    void restoreCavityPhases();
-    void changeDT();
-    void selectDT();
-    void autophaseCavities(const BorisPusher &pusher);
-    void findStartPosition(const BorisPusher &pusher);
-
-
-
-
-    typedef FTps<double, PSdim> series_t;
-    typedef FVps<double, PSdim> map_t;
-    typedef FMatrix<double, PSdim, PSdim> fMatrix_t;
-    typedef FMatrix<std::complex<double>, PSdim, PSdim> cfMatrix_t;
-
-    series_t x, px, y, py, delta; //z
-
-
-    struct structMapTracking {
-                	std::string elementName;
-                	map_t elementMap;
-                	std::size_t nSlices;
-                	double elementPos;
-                	double stepSize;
-                };
-
-
-    series_t createHamiltonian(std::shared_ptr<Component> element, double& stepSize, std::size_t& nSlices);
-    void defMapTrackingElement(std::shared_ptr<Component> element, structMapTracking& elSrct, std::list<structMapTracking>& mBL);
-    void fillDrift(std::list<structMapTracking>& mapBeamLine,double& elementPos, double& undefSpace);
-
-    void setHamiltonianDrift(series_t& H, double& beta0, double& gamma0);
-    void setHamiltonianSBend(series_t& H, double& beta0, double& gamma0, double& q, double& h, double& K0 );
-    void setHamiltonianRBend(series_t& H, double& beta0, double& gamma0, double& q, double& h, double& K0 );
-    void setHamiltonianQuadrupole(series_t& H, double& beta0, double& gamma0, double& q, double& K1 );
-
-    void eigenDecomp(fMatrix_t& M, cfMatrix_t& eigenVal, cfMatrix_t& eigenVec, cfMatrix_t& invEigenVec);
-    void linTAnalyze(fMatrix_t& tMatrix);
-    void linSigAnalyze();
-    cfMatrix_t getBlockDiagonal(fMatrix_t& M, cfMatrix_t& eigenVecM, cfMatrix_t& invEigenVecM);
-
-    void dumpStats(long long step, bool psDump, bool statDump);
-    void writePhaseSpace(const long long step, bool psDump, bool statDump);
-    void setTime();
 
+    /*
+    void insertFringeField(SBend* pSBend, std::list<structMapTracking>& mBL, double& beta0,
+             double& gamma0, double& P0, double& q, std::array<double,2>& entrFringe, std::string e);
+     */
 private:
 
+
     // Not implemented.
-    ThickTracker();
-    ThickTracker(const ThickTracker &);
-    void operator=(const ThickTracker &);
-    
-    void trackParticles_m(
-#ifdef PHIL_WRITE
-        std::ofstream& outfile,
-#endif
-                          const std::list<structMapTracking>& mapBeamLine);
-
-    // Fringe fields for entrance and exit of magnetic elements.
+    ThickTracker() = delete;
+    ThickTracker(const ThickTracker &)   = delete;
+    void operator=(const ThickTracker &) = delete;
+
+    void throwElementError_m(std::string element) {
+        throw LogicalError("ThickTracker::execute()",
+                           "Element '" + element + "' not supported.");
+    }
+
+
+
+    /*!
+     * Tests the order of the elements in the beam line according to their position
+     */
+    void checkElementOrder_m();
+
+    /*!
+     * Inserts Drift maps in undefined beam line sections
+     */
+    void fillGaps_m();
+
+    /*!
+     * Tracks itsBunch_m trough beam line
+     */
+    void track_m();
+
+    particle_t particleToVector_m(const Vector_t& R,
+                                  const Vector_t& P) const;
+
+    void vectorToParticle_m(const particle_t& particle,
+                            Vector_t& R,
+                            Vector_t& P) const;
+
+    /*!
+     * Advances itsBunch_m trough map
+     * @param map Map of slice
+     */
+    void advanceParticles_m(const map_t& map);
+
+    /*!
+     * Applies map on particle
+     * @param particle tracked particle
+     * @param map Map of slice
+     */
+    void updateParticle_m(particle_t& particle,
+                          const map_t& map);
+
+    /*!
+     * Dumps bunch in .stat or .h5 files
+     */
+    void dump_m();
+
+
+    /*!
+     * Updates itsBunch_m
+     * @param spos position of tracking
+     * @param step stepsize of applied map
+     */
+    void update_m(const double& spos,
+                  const std::size_t& step);
+
+
+    /*!
+     * Writes map (Transfermap) in .map file
+     * @map text for .map file
+     */
+    void write_m(const map_t& map);
+
+    /*!
+     * Concatenate map x and y
+     * \f[
+     *  y = x \circ y
+     * \f]
+     * @param x
+     * @param y is result
+     */
+    void concatenateMaps_m(const map_t& x, map_t& y);
+
+
+    /*!
+     * Tracks Dispersion along beam line
+     * writes it in .dispersion file.
+     *
+     *  \f[
+     *  \begin{pmatrix} \eta_{x} \\ \eta_{p_x} \end{pmatrix}_{s_1}
+     *  =
+     *  \begin{pmatrix} R_{11} & R_{12} \\ R_{21} & R_{22} \end{pmatrix}
+     *  \cdot
+     *  \begin{pmatrix} \eta_{x} \\ \eta_{p_x} \end{pmatrix}_{s_0}
+     *  +
+     *  \begin{pmatrix} R_{16} \\ R_{26} \end{pmatrix}
+     *  \f]
+     *
+     * @param tempMatrix accumulated Transfer map \f$R\f$ at pos
+     * @param initialVal initial Dispersion { \f$\eta_{x0},\, \eta_{p_x0},\, \eta_{y0},\, \eta_{p_y0} \f$}
+     * @param pos position of tracking
+     */
+    void advanceDispersion_m(fMatrix_t tempMatrix,
+                             FMatrix<double, 1, 4> initialVal,
+                             double pos);
+
+    /*! :TODO:
+     * Fringe fields for entrance of SBEND.
+     *
+     * @param edge
+     * @param curve
+     * @param field
+     * @param scale
+     */
     void applyEntranceFringe(double edge, double curve,
                              const BMultipoleField &field, double scale);
+
+    /*! :TODO:
+     * Fringe fields for exit of SBEND.
+     *
+     * @param edge
+     * @param curve
+     * @param field
+     * @param scale
+     */
     void applyExitFringe(double edge, double curve,
                          const BMultipoleField &field, double scale);
 
+
+    Hamiltonian hamiltonian_m;
+    MapAnalyser mapAnalyser_m;
+
+
     Vector_t RefPartR_m;
     Vector_t RefPartP_m;
 
@@ -291,136 +366,220 @@ private:
 
     OpalBeamline itsOpalBeamline_m;
 
-    double pathLength_m;
-
-    /// where to start
-    double zstart_m;
-
-
-    /// where to stop
-    std::queue<double> zStop_m;
+    double zstart_m;        ///< Start of beam line
+    double zstop_m;         ///< End of beam line
+    double threshold_m;     ///< Threshold for element overlaps and gaps
+    beamline_t elements_m;  ///< elements in beam line
 
-    double dtCurrentTrack_m;
-    std::queue<double> dtAllTracks_m;
+    CoordinateSystemTrafo referenceToLabCSTrafo_m;
 
-    /// The maximal number of steps the system is integrated per TRACK
-    std::queue<unsigned long long> localTrackSteps_m;
 
-    CoordinateSystemTrafo referenceToLabCSTrafo_m;
+    int truncOrder_m; ///< truncation order of map tracking
 
-    bool globalEOL_m;
-    
-    
-    int truncOrder_m;
+    IpplTimings::TimerRef mapCreation_m;    ///< creation of elements_m
+    IpplTimings::TimerRef mapCombination_m; ///< map accumulation along elements_m -> Transfermap
+    IpplTimings::TimerRef mapTracking_m;    ///< track particles trough maps of elements_m
 };
 
 inline void ThickTracker::visitAlignWrapper(const AlignWrapper &wrap) {
-    itsOpalBeamline_m.visit(wrap, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(wrap, *this, itsBunch_m);
+    this->throwElementError_m("AlignWrapper");
 }
 
 inline void ThickTracker::visitBeamBeam(const BeamBeam &bb) {
-    itsOpalBeamline_m.visit(bb, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(bb, *this, itsBunch_m);
+    this->throwElementError_m("BeamBeam");
 }
 
 
 inline void ThickTracker::visitCCollimator(const CCollimator &coll) {
-    itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
+    this->throwElementError_m("CCollimator");
 }
 
 
 inline void ThickTracker::visitCorrector(const Corrector &corr) {
-    itsOpalBeamline_m.visit(corr, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(corr, *this, itsBunch_m);
+    this->throwElementError_m("Corrector");
 }
 
 
 inline void ThickTracker::visitDegrader(const Degrader &deg) {
-    itsOpalBeamline_m.visit(deg, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(deg, *this, itsBunch_m);
+    this->throwElementError_m("Degrader");
 }
 
 
 inline void ThickTracker::visitDiagnostic(const Diagnostic &diag) {
-    itsOpalBeamline_m.visit(diag, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(diag, *this, itsBunch_m);
+    this->throwElementError_m("Diagnostic");
 }
 
 
 inline void ThickTracker::visitDrift(const Drift &drift) {
     itsOpalBeamline_m.visit(drift, *this, itsBunch_m);
+
+
+    double gamma = itsReference.getGamma();
+    std::size_t nSlices = drift.getNSlices();
+    double length       = drift.getElementLength();
+
+    elements_m.push_back(std::make_tuple(hamiltonian_m.drift(gamma),
+                                         nSlices,
+                                         length));
 }
 
 
 inline void ThickTracker::visitFlexibleCollimator(const FlexibleCollimator &coll) {
-    itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
+    this->throwElementError_m("FlexibleCollimator");
 }
 
 
 inline void ThickTracker::visitLambertson(const Lambertson &lamb) {
-    itsOpalBeamline_m.visit(lamb, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(lamb, *this, itsBunch_m);
+    this->throwElementError_m("Lambertson");
 }
 
 
 inline void ThickTracker::visitMarker(const Marker &marker) {
-    itsOpalBeamline_m.visit(marker, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(marker, *this, itsBunch_m);
+//     this->throwElementError_m("Marker");
 }
 
 
 inline void ThickTracker::visitMonitor(const Monitor &mon) {
-    itsOpalBeamline_m.visit(mon, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(mon, *this, itsBunch_m);
+    this->throwElementError_m("Monitor");
 }
 
 
 inline void ThickTracker::visitMultipole(const Multipole &mult) {
     itsOpalBeamline_m.visit(mult, *this, itsBunch_m);
+
+    std::size_t nSlices = mult.getNSlices();
+    double length       = mult.getElementLength();
+    double gamma        = itsReference.getGamma();
+    double p0           = itsReference.getP();
+    double q            = itsReference.getQ(); // particle change [e]
+    double gradB        = mult.getField().getNormalComponent(2) * ( Physics::c/ p0 ); // [T / m]
+    //FIXME remove the next line
+    gradB = std::round(gradB*1e6)*1e-6;
+
+    double k1           = gradB * q *Physics::c / p0; // [1 / m^2]
+
+    elements_m.push_back(std::make_tuple(hamiltonian_m.quadrupole(gamma, q, k1),
+                                         nSlices,
+                                         length));
 }
 
 inline void ThickTracker::visitProbe(const Probe &prob) {
-    itsOpalBeamline_m.visit(prob, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(prob, *this, itsBunch_m);
+    this->throwElementError_m("Probe");
 }
 
 
 inline void ThickTracker::visitRBend(const RBend &bend) {
-    itsOpalBeamline_m.visit(bend, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(bend, *this, itsBunch_m);
+    this->throwElementError_m("RBend");
 }
 
 
 inline void ThickTracker::visitRFCavity(const RFCavity &as) {
-    itsOpalBeamline_m.visit(as, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(as, *this, itsBunch_m);
+    this->throwElementError_m("RFCavity");
 }
 
 inline void ThickTracker::visitTravelingWave(const TravelingWave &as) {
-    itsOpalBeamline_m.visit(as, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(as, *this, itsBunch_m);
+    this->throwElementError_m("TravelingWave");
 }
 
 
 inline void ThickTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
-    itsOpalBeamline_m.visit(rfq, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(rfq, *this, itsBunch_m);
+    this->throwElementError_m("RFQuadrupole");
 }
 
 inline void ThickTracker::visitSBend(const SBend &bend) {
     itsOpalBeamline_m.visit(bend, *this, itsBunch_m);
+
+    double q     = itsReference.getQ(); // particle change [e]
+    double ekin = bend.getDesignEnergy();
+
+    double m     = itsReference.getM(); // eV / c^2
+    double gamma = ekin / m + 1.0;
+    double beta  = std::sqrt(1.0 - 1.0 / ( gamma * gamma ) );
+    double p0    = gamma * beta * m; // eV / c
+
+    double B     = bend.getB() * Physics::c / p0; // T = V * s / m^2
+    double r     = std::abs(p0   / ( q * B * Physics::c )); // m
+
+    double k0    = B * q * Physics::c / p0; // V * s * e * m / (m^2 * s * c )
+
+    // [1/m]
+    double h = 1.0 / r;
+
+    double L = bend.getElementLength();
+
+    if ( k0 < 0.0 )
+        h *= -1.0;
+
+    std::size_t nSlices = bend.getNSlices();
+
+    double arclength    = 2.0 * r * std::asin( L / ( 2.0 * r ) ); //bend.getEffectiveLength();
+
+    // Fringe fields currently not working
+    //FIXME e1 not initialised
+    //insert Entrance Fringefield
+    double e1 = bend.getEntranceAngle();
+    if (std::abs(e1) > 1e-6){
+        elements_m.push_back(std::make_tuple(hamiltonian_m.fringeField(e1, h),
+                                             1, 0.0));
+    }
+
+    //insert Dipole "body"
+    elements_m.push_back(std::make_tuple(hamiltonian_m.sbend(gamma, h, k0),
+                                         nSlices,
+                                         arclength));
+
+    //FIXME e2 not initialised
+    //insert Exit Fringe field
+    double e2 = bend.getExitAngle();
+    if (std::abs(e2) > 1e-6){
+        elements_m.push_back(std::make_tuple(hamiltonian_m.fringeField(e2, h),
+                                             1, 0.0));
+    }
+
 }
 
 
 inline void ThickTracker::visitSeparator(const Separator &sep) {
-    itsOpalBeamline_m.visit(sep, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(sep, *this, itsBunch_m);
+    this->throwElementError_m("Separator");
 }
 
 
 inline void ThickTracker::visitSeptum(const Septum &sept) {
-    itsOpalBeamline_m.visit(sept, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(sept, *this, itsBunch_m);
+    this->throwElementError_m("Septum");
 }
 
 
 inline void ThickTracker::visitSolenoid(const Solenoid &solenoid) {
-    itsOpalBeamline_m.visit(solenoid, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(solenoid, *this, itsBunch_m);
+    this->throwElementError_m("Solenoid");
 }
 
 
 inline void ThickTracker::visitParallelPlate(const ParallelPlate &pplate) {
-    itsOpalBeamline_m.visit(pplate, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(pplate, *this, itsBunch_m);
+    this->throwElementError_m("ParallelPlate");
 }
 
 inline void ThickTracker::visitCyclotronValley(const CyclotronValley &cv) {
-    itsOpalBeamline_m.visit(cv, *this, itsBunch_m);
+//     itsOpalBeamline_m.visit(cv, *this, itsBunch_m);
+    this->throwElementError_m("CyclotronValley");
 }
 
 #endif // OPAL_ThickTracker_HH