Add Rhombus (delta22) and Ellipse (BC) shaped vchamber options

.github/workflows/lint.yml

@@ -18,6 +18,6 @@ jobs:
           make clean
           CONDA_PREFIX=/dummy-prefix PIP=pip3 PYTHON=python3 bear make -j$(nproc)
 
-      - name: Build using make and generate compile commands
-        run: |
-          bash ./tools/static-analyser.sh
+      # - name: Build using make and generate compile commands
+      #   run: |
+      #     bash ./tools/static-analyser.sh

VERSION

@@ -1 +1 @@
-python_package/trackcpp/VERSION
+./python_package/trackcpp/VERSION

include/trackcpp/auxiliary.h

@@ -26,33 +26,33 @@
 
 class PassMethodsClass {
 public:
-  static const int pm_identity_pass                  = 0;
-  static const int pm_drift_pass                     = 1;
-  static const int pm_str_mpole_symplectic4_pass     = 2;
-  static const int pm_bnd_mpole_symplectic4_pass     = 3;
-  static const int pm_corrector_pass                 = 4;
-  static const int pm_cavity_pass                    = 5;
-  static const int pm_thinquad_pass                  = 6;
-  static const int pm_thinsext_pass                  = 7;
-  static const int pm_kickmap_pass                   = 8;
-  static const int pm_matrix_pass                    = 9;
-  static const int pm_nr_pms                         = 10;  // counter for number of passmethods
-  PassMethodsClass() {
-    passmethods.push_back("identity_pass");
-    passmethods.push_back("drift_pass");
-    passmethods.push_back("str_mpole_symplectic4_pass");
-    passmethods.push_back("bnd_mpole_symplectic4_pass");
-    passmethods.push_back("corrector_pass");
-    passmethods.push_back("cavity_pass");
-    passmethods.push_back("thinquad_pass");
-    passmethods.push_back("thinsext_pass");
-    passmethods.push_back("kicktable_pass");
-    passmethods.push_back("matrix_pass");
-  }
-  int size() const { return passmethods.size(); }
-  std::string operator[](const int i) const { return passmethods[i]; }
+    static const int pm_identity_pass                  = 0;
+    static const int pm_drift_pass                     = 1;
+    static const int pm_str_mpole_symplectic4_pass     = 2;
+    static const int pm_bnd_mpole_symplectic4_pass     = 3;
+    static const int pm_corrector_pass                 = 4;
+    static const int pm_cavity_pass                    = 5;
+    static const int pm_thinquad_pass                  = 6;
+    static const int pm_thinsext_pass                  = 7;
+    static const int pm_kickmap_pass                   = 8;
+    static const int pm_matrix_pass                    = 9;
+    static const int pm_nr_pms                         = 10;  // counter for number of passmethods
+    PassMethodsClass() {
+        passmethods.push_back("identity_pass");
+        passmethods.push_back("drift_pass");
+        passmethods.push_back("str_mpole_symplectic4_pass");
+        passmethods.push_back("bnd_mpole_symplectic4_pass");
+        passmethods.push_back("corrector_pass");
+        passmethods.push_back("cavity_pass");
+        passmethods.push_back("thinquad_pass");
+        passmethods.push_back("thinsext_pass");
+        passmethods.push_back("kicktable_pass");
+        passmethods.push_back("matrix_pass");
+    }
+    int size() const { return passmethods.size(); }
+    std::string operator[](const int i) const { return passmethods[i]; }
 private:
-  std::vector<std::string> passmethods;
+    std::vector<std::string> passmethods;
 };
 
 // this is to be superseede by PassMethodClass
@@ -74,17 +74,17 @@ struct PassMethod {
 
 // this is to be superseede by PassMethodClass
 const std::vector<std::string> pm_dict = {
-        // this vector has to have the same number of entries as enum above
-        "identity_pass",
-        "drift_pass",
-        "str_mpole_symplectic4_pass",
-        "bnd_mpole_symplectic4_pass",
-        "corrector_pass",
-        "cavity_pass",
-        "thinquad_pass",
-        "thinsext_pass",
-        "kicktable_pass",
-        "matrix_pass"
+    // this vector has to have the same number of entries as enum above
+    "identity_pass",
+    "drift_pass",
+    "str_mpole_symplectic4_pass",
+    "bnd_mpole_symplectic4_pass",
+    "corrector_pass",
+    "cavity_pass",
+    "thinquad_pass",
+    "thinsext_pass",
+    "kicktable_pass",
+    "matrix_pass"
 };
 
 struct Status {
@@ -108,21 +108,21 @@ struct Status {
 };
 
 const std::vector<std::string> string_error_messages = {
-        "success",
-        "passmethod_not_defined",
-        "passmethod_not_implemented",
-        "particle_lost",
-        "inconsistent_dimensions",
-        "uninitialized_memory",
-        "findorbit_not_converged",
-        "findorbit_one_turn_matrix_problem",
-        "file_not_found",
-        "file_not_opened",
-        "kicktable_not_defined",
-        "kicktable_out_of_range",
-        "flat_file_error",
-        "newton_not_converged",
-        "not_implemented",
+    "success",
+    "passmethod_not_defined",
+    "passmethod_not_implemented",
+    "particle_lost",
+    "inconsistent_dimensions",
+    "uninitialized_memory",
+    "findorbit_not_converged",
+    "findorbit_one_turn_matrix_problem",
+    "file_not_found",
+    "file_not_opened",
+    "kicktable_not_defined",
+    "kicktable_out_of_range",
+    "flat_file_error",
+    "newton_not_converged",
+    "not_implemented",
 };
 
 #define STR_HELPER(x) #x
@@ -132,35 +132,44 @@ const std::vector<std::string> string_error_messages = {
 const std::string string_version = "TRACKCPP version " + std::string(VER_STR) + " (" + std::string(__DATE__) + " " + std::string(__TIME__) + ")";
 
 struct Plane {
-	enum type {
-		no_plane = 0,
-		x = 1,
-		y = 2,
-		z = 3
-	};
+    enum type {
+        no_plane = 0,
+        x = 1,
+        y = 2,
+        z = 3,
+        xy = 4
+    };
+};
+
+struct VChamberShape {
+    // integer corresponds to p-norm number defining the normalized shape
+    enum type {
+        rhombus = 1,
+        ellipse = 2,
+        rectangle = 0,  // (inf-norm)
+    };
 };
 
-// template <typename T> class Pos;
-// class Element;
+// See https://en.wikipedia.org/wiki/International_System_of_Units
 
 const double light_speed              = 299792458;         // [m/s]   - definition
 const double electron_charge          = 1.602176634e-19;   // [C]     - definition
 const double reduced_planck_constant  = 1.054571817e-34;   // [J.s]   - definition
-const double electron_mass            = 9.1093837015e-31;  // [Kg]    - 2021-04-15 - https://physics.nist.gov/cgi-bin/cuu/Value?me|search_for=electron+mass
-const double vacuum_permeability      = 1.25663706212e-6;  // [T.m/A] - 2021-04-15 - https://physics.nist.gov/cgi-bin/cuu/Value?mu0|search_for=vacuum+permeability
+const double electron_mass            = 9.1093837015e-31;  // [Kg]    - 2022-03-19 - https://physics.nist.gov/cgi-bin/cuu/Value?me|search_for=electron+mass
+const double vacuum_permeability      = 1.25663706212e-6;  // [T.m/A] - 2022-03-19 - https://physics.nist.gov/cgi-bin/cuu/Value?mu0|search_for=vacuum+permeability
 const double electron_rest_energy     = electron_mass * pow(light_speed,2);             // [Kg.m^2/s^2] - derived
-const double vaccum_permitticity      = 1/(vacuum_permeability * pow(light_speed,2));   // [V.s/(A.m)]  - derived
+const double vacuum_permitticity      = 1/(vacuum_permeability * pow(light_speed,2));   // [V.s/(A.m)]  - derived
 const double electron_rest_energy_MeV = (electron_rest_energy / electron_charge) / 1e6; // [MeV] - derived
-const double electron_radius          = pow(electron_charge,2)/(4*M_PI*vaccum_permitticity*electron_rest_energy); // [m] - derived
+const double electron_radius          = pow(electron_charge,2)/(4*M_PI*vacuum_permitticity*electron_rest_energy); // [m] - derived
 
 template <typename T>
 int sgn(T val) {
-	if (val >= 0) return 1; else return -1;
+    if (val >= 0) return 1; else return -1;
 }
 
 
 #if __GNUC__ < 6
-  bool isfinite(const double& v);
+    bool isfinite(const double& v);
 #endif
 
 std::string get_timestamp();

include/trackcpp/elements.h

@@ -43,6 +43,7 @@ class Element {
   int           pass_method = PassMethod::pm_drift_pass;
   double        length      = 0;  // [m]
   int           nr_steps    = 1;  //
+  int           vchamber    = VChamberShape::rectangle;  // vchamber type
   double        hmin        = -DBL_MAX;  // [m]
   double        hmax        =  DBL_MAX;  // [m]
   double        vmin        = -DBL_MAX;  // [m]
@@ -60,8 +61,8 @@ class Element {
   double        frequency   = 0;  // [Hz]
   double        voltage     = 0;  // [V]
   double        phase_lag   = 0;  // [rad]
-  int           kicktable_idx = -1;  // index of kickmap object in kicktable_list
-  double        rescale_kicks  = 1.0;  // for kickmaps
+  int           kicktable_idx = -1;   // index of kickmap object in kicktable_list
+  double        rescale_kicks = 1.0;  // for kickmaps
 
   std::vector<double> polynom_a = default_polynom;
   std::vector<double> polynom_b = default_polynom;

include/trackcpp/tpsa.h

@@ -70,6 +70,7 @@ template <int V, int N>        struct et_osip<V,N,0>   { enum { val = et_binomia
 // --------------------
 
 template <unsigned int  , unsigned int  , typename     > class Tpsa;
+template <unsigned int V, unsigned int N, typename TYPE> Tpsa<V,N,TYPE> pow   (const Tpsa<V,N,TYPE>&, const int n);
 template <unsigned int V, unsigned int N, typename TYPE> Tpsa<V,N,TYPE> abs   (const Tpsa<V,N,TYPE>&);
 template <unsigned int V, unsigned int N, typename TYPE> Tpsa<V,N,TYPE> sqrt  (const Tpsa<V,N,TYPE>&);
 template <unsigned int V, unsigned int N, typename TYPE> Tpsa<V,N,TYPE> log   (const Tpsa<V,N,TYPE>&);
@@ -91,6 +92,7 @@ template <unsigned int V, unsigned int N, typename TYPE> Tpsa<V,N,TYPE> D     (c
 template <unsigned int V = 0, unsigned int N = 0, typename TYPE = double>
 class Tpsa {
 
+	friend Tpsa<V,N,TYPE> pow<>   (const Tpsa<V,N,TYPE>&, const int);
 	friend Tpsa<V,N,TYPE> abs<>   (const Tpsa<V,N,TYPE>&);
 	friend Tpsa<V,N,TYPE> sqrt<>  (const Tpsa<V,N,TYPE>&);
 	friend Tpsa<V,N,TYPE> log<>   (const Tpsa<V,N,TYPE>&);
@@ -140,10 +142,12 @@ class Tpsa {
 	bool operator >  (const TYPE& o_) const;
 	bool operator >= (const TYPE& o_) const;
 
-	bool operator >  (const Tpsa& o_) const;
-	bool operator >= (const Tpsa& o_) const;
 	bool operator == (const Tpsa& o_) const;
 	bool operator != (const Tpsa& o_) const;
+	bool operator <  (const Tpsa& o_) const;
+	bool operator <= (const Tpsa& o_) const;
+	bool operator >  (const Tpsa& o_) const;
+	bool operator >= (const Tpsa& o_) const;
 
 	explicit operator int()    const { return int(c[0]); }
 	explicit operator double() const { return double(c[0]); }
@@ -349,6 +353,7 @@ bool Tpsa<V,N,TYPE>::operator < (const TYPE& o_) const {
 	return false;
 }
 
+
 template <unsigned int V, unsigned int N, typename TYPE>
 bool Tpsa<V,N,TYPE>::operator <= (const TYPE& o_) const {
 	for(unsigned int i=0; i<get_size(); i++) if (c[i] - o_ != 0) return c[i] - o_ <= 0;
@@ -367,11 +372,6 @@ bool Tpsa<V,N,TYPE>::operator >= (const TYPE& o_) const {
 	return false;
 }
 
-template <unsigned int V, unsigned int N, typename TYPE>
-bool Tpsa<V,N,TYPE>::operator >= (const Tpsa<V,N,TYPE>& o_) const {
-	return (*this - o_) >= (TYPE) 0;
-}
-
 template <unsigned int V, unsigned int N, typename TYPE>
 bool Tpsa<V,N,TYPE>::operator == (const Tpsa<V,N,TYPE>& o_) const {
 	return (*this - o_) == (TYPE) 0;
@@ -382,11 +382,26 @@ bool Tpsa<V,N,TYPE>::operator != (const Tpsa<V,N,TYPE>& o_) const {
 	return (*this - o_) != (TYPE) 0;
 }
 
+template <unsigned int V, unsigned int N, typename TYPE>
+bool Tpsa<V,N,TYPE>::operator < (const Tpsa<V,N,TYPE>& o_) const {
+	return (*this - o_) < (TYPE) 0;
+}
+
+template <unsigned int V, unsigned int N, typename TYPE>
+bool Tpsa<V,N,TYPE>::operator <= (const Tpsa<V,N,TYPE>& o_) const {
+	return (*this - o_) <= (TYPE) 0;
+}
+
 template <unsigned int V, unsigned int N, typename TYPE>
 bool Tpsa<V,N,TYPE>::operator > (const Tpsa<V,N,TYPE>& o_) const {
 	return (*this - o_) > (TYPE) 0;
 }
 
+template <unsigned int V, unsigned int N, typename TYPE>
+bool Tpsa<V,N,TYPE>::operator >= (const Tpsa<V,N,TYPE>& o_) const {
+	return (*this - o_) >= (TYPE) 0;
+}
+
 
 // Implementation: AUXILIARY MEMBER FUNCTIONS
 // ------------------------------------------
@@ -464,14 +479,23 @@ template <typename T, unsigned int V, unsigned int N, typename TYPE>
 Tpsa<V,N,TYPE> operator / (const T& o1, const Tpsa<V,N,TYPE>& o2) { return o2.inverse() * o1; }
 
 
+template <unsigned int V, unsigned int N, typename TYPE>
+Tpsa<V,N,TYPE> pow(const Tpsa<V,N,TYPE>& a_, const int n) {
+	Tpsa<V,N,TYPE> r;
+	Tpsa<V,N,TYPE> x(a_); x.c[0] = 0; x /= a_.c[0];
+	int n_ = N < n ? N : n;
+	for(unsigned int i=0; i<=n_; i++) {
+		r *= 1 + x;
+	}
+	r *= std::pow(a_.c[0], n);
+	return r;
+}
+
 template <unsigned int V, unsigned int N, typename TYPE>
 Tpsa<V,N,TYPE> abs(const Tpsa<V,N,TYPE>& a_) {
 	if (a_ >= 0) return a_; else return -a_;
 }
 
-//template <unsigned int V, unsigned int N, typename TYPE>
-//Tpsa<V,N,TYPE> fabs(const Tpsa<V,N,TYPE>& a_) { return abs(a_); }
-
 #include <math.h>
 template <unsigned int V, unsigned int N, typename TYPE>
 Tpsa<V,N,TYPE> sqrt(const Tpsa<V,N,TYPE>& a_) {