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/Hamiltonian.cpp

@@ -13,9 +13,10 @@ Hamiltonian::Hamiltonian(int truncOrder) : truncOrder_m(truncOrder)
     delta = series_t::makeVariable(5);
 }
 
-Hamiltonian::series_t Hamiltonian::drift(double& beta0,
-                            double& gamma0)
+Hamiltonian::series_t Hamiltonian::drift(const double& gamma0)
 {
+    double beta0 = this->getBeta_m(gamma0);
+    
     return ( delta / beta0 )
             - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
                     - ( px*px )
@@ -28,7 +29,7 @@ Hamiltonian::series_t Hamiltonian::rbend(double& beta0,
                             double& gamma0,
                             double& q,
                             double& h,
-                            double& K0)
+                            double& k0)
 {
 
 
@@ -45,18 +46,17 @@ Hamiltonian::series_t Hamiltonian::rbend(double& beta0,
                     - ( py*py )
                     - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m
             ))
-    + K0 * x * (1. + 0.5 * h* x);
+    + k0 * x * (1. + 0.5 * h* x);
 
 
 }
     
     
-Hamiltonian::series_t Hamiltonian::sbend(double& beta0,
-                                         double& gamma0,
-                                         double& h,
-                                         double& K0)
+Hamiltonian::series_t Hamiltonian::sbend(const double& gamma0,
+                                         const double& h,
+                                         const double& K0)
 {
-
+    double beta0 = this->getBeta_m(gamma0);
 
     return ( delta / beta0 )
                     - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
@@ -113,17 +113,22 @@ Hamiltonian::series_t Hamiltonian::bendFringe(
 }
 
 
-Hamiltonian::series_t Hamiltonian::quadrupole(double& beta0,
-                double& gamma0,
-                double& q,
-                double& K1)
+Hamiltonian::series_t Hamiltonian::quadrupole(const double& gamma0,
+                                              const double& q,
+                                              const double& k1)
 {
-
+    double beta0 = this->getBeta_m(gamma0);
+    
     return ( delta / beta0 )
     - sqrt ((1./ beta0 + delta ) *(1./ beta0 + delta)
             - ( px*px )
             - ( py*py )
             - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m+1
     )
-    + 0.5 * K1 * (x*x - y*y);
+    + 0.5 * k1 * (x*x - y*y);
+}
+
+
+double Hamiltonian::getBeta_m(const double& gamma) {
+    return std::sqrt(1.0 - 1.0 / (gamma * gamma) );
 }