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

@@ -112,7 +112,7 @@ public:
     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.
@@ -141,7 +141,7 @@ public:
                           const std::vector<double> &zstop,
                           const std::vector<double> &dt,
                           const int& truncOrder);
-    
+
     virtual ~ThickTracker();
 
 
@@ -229,48 +229,48 @@ public:
              double& gamma0, double& P0, double& q, std::array<double,2>& entrFringe, std::string e);
      */
 private:
-    
+
 
     // Not implemented.
     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
@@ -278,12 +278,12 @@ private:
      */
     void updateParticle_m(particle_t& particle,
                           const map_t& map);
-    
+
     /*!
      * Dumps bunch in .stat or .h5 files
      */
     void dump_m();
-    
+
 
     /*!
      * Updates itsBunch_m
@@ -292,25 +292,25 @@ private:
      */
     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 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.
@@ -332,7 +332,7 @@ private:
     void advanceDispersion_m(fMatrix_t tempMatrix,
                              FMatrix<double, 1, 4> initialVal,
                              double pos);
-    
+
     /*! :TODO:
      * Fringe fields for entrance of SBEND.
      *
@@ -358,15 +358,15 @@ private:
 
     Hamiltonian hamiltonian_m;
     MapAnalyser mapAnalyser_m;
-    
-    
+
+
     Vector_t RefPartR_m;
     Vector_t RefPartP_m;
-    
+
     DataSink *itsDataSink_m;
 
     OpalBeamline itsOpalBeamline_m;
-    
+
     double zstart_m;        ///< Start of beam line
     double zstop_m;         ///< End of beam line
     double threshold_m;     ///< Threshold for element overlaps and gaps
@@ -419,12 +419,12 @@ inline void ThickTracker::visitDiagnostic(const Diagnostic &diag) {
 
 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));
@@ -457,7 +457,7 @@ inline void ThickTracker::visitMonitor(const Monitor &mon) {
 
 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();
@@ -466,9 +466,9 @@ inline void ThickTracker::visitMultipole(const Multipole &mult) {
     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));
@@ -504,47 +504,46 @@ inline void ThickTracker::visitRFQuadrupole(const RFQuadrupole &rfq) {
 
 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 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
-    double e1 = bend.getEntranceAngle();
-    
-    elements_m.push_back(std::make_tuple(hamiltonian_m.fringeField(e1, h),
-                                         1, 0.0));
-    
+    //double e1 = bend.getEntranceAngle();
+    //elements_m.push_back(std::make_tuple(hamiltonian_m.fringeField(e1, h),
+    //                                     1, 0.0));
+
     elements_m.push_back(std::make_tuple(hamiltonian_m.sbend(gamma, h, k0),
                                          nSlices,
                                          arclength));
-    
+
     //FIXME e2 not initialised
-    double e2 = bend.getExitAngle();
-    
-    elements_m.push_back(std::make_tuple(hamiltonian_m.fringeField(e2, h),
-                                         1, 0.0));
+    //double e2 = bend.getExitAngle();
+    //elements_m.push_back(std::make_tuple(hamiltonian_m.fringeField(e2, h),
+    //                                     1, 0.0));
 }