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

+#include "Hamiltonian.h"
+
+Hamiltonian::Hamiltonian(int truncOrder) : truncOrder_m(truncOrder)
+{
+        
+    series_t::setGlobalTruncOrder(truncOrder+1);
+        
+    x     = series_t::makeVariable(0);
+    px    = series_t::makeVariable(1);
+    y     = series_t::makeVariable(2);
+    py    = series_t::makeVariable(3);
+    z     = series_t::makeVariable(4);
+    delta = series_t::makeVariable(5);
+}
+
+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 )
+                    - ( py*py )
+                    - 1./( beta0 * beta0 * gamma0 * gamma0 ),truncOrder_m+1);
+}
+
+    
+Hamiltonian::series_t Hamiltonian::rbend(double& beta0,
+                            double& gamma0,
+                            double& q,
+                            double& h,
+                            double& k0)
+{
+
+
+
+        return ( delta / beta0 )
+        - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
+                    - ( px*px )
+                    - ( py*py )
+                    - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m+1
+            ))
+    - (h * x)
+    * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
+                    - ( px*px )
+                    - ( py*py )
+                    - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m
+            ))
+    + k0 * x * (1. + 0.5 * h* x);
+
+
+}
+    
+    
+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)
+                                    - ( px*px )
+                                    - ( py*py )
+                                    - 1./( beta0*beta0 * gamma0*gamma0 ),(truncOrder_m+1)
+                            ))
+                    - (h * x)
+                    * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
+                                    - ( px*px )
+                                    - ( py*py )
+                                    - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m
+                            ))
+                    + K0 * x * (1. + 0.5 * h* x);
+}
+
+
+Hamiltonian::series_t Hamiltonian::bendFringe(
+                double& beta0,
+                double& gamma0,
+                double& h,
+                double& k0,
+                series_t& ax,
+                series_t& az)
+{
+    if (truncOrder_m == 2){
+        return ( delta / beta0 )
+                        - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
+                                - ( px*px)
+                                - ( py*py )
+                                - 1./( beta0 * beta0 * gamma0 * gamma0 ),truncOrder_m+1
+                        ) - az;
+    }else{
+        return ( delta / beta0 )
+            - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
+                    - ( px*px - 2.0*px*ax - ax*ax)
+                    - ( py*py )
+                    - 1./( beta0 * beta0 * gamma0 * gamma0 ),truncOrder_m+1
+            ) - az;
+    }
+    //std::cout << H << std::endl;
+    //std::cout << H.getMaxOrder() << std::endl;
+    /*H=( delta / beta0 )
+                        - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
+                                        - ( px - ax )*( px - ax )
+                                        - ( py*py )
+                                        - 1./( beta0*beta0 * gamma0*gamma0 ),(truncOrder_m+1)))
+                        - (h * x)
+                        * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
+                                        - ( px - ax )*( px - ax )
+                                        - ( py*py )
+                                        - 1./( beta0*beta0 * gamma0*gamma0 ), truncOrder_m))
+                         + h * x * az;*/
+}
+
+
+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);
+}
+
+
+Hamiltonian::series_t Hamiltonian::fringeField(const double& phi,
+                                               const double& k0)
+{
+    if ( truncOrder_m > 1 ) {
+        //TODO higher order thin lens approximation
+    }
+    // else
+    
+    double angle = phi;
+    if ( k0 < 0.0 )
+        angle = -phi;
+    return -0.5 * (x * x - y * y) * k0 * std::tan(angle);
+}
+
+
+double Hamiltonian::getBeta_m(const double& gamma) {
+    return std::sqrt(1.0 - 1.0 / (gamma * gamma) );
+}