Add ID symplectic passmethod

include/trackcpp/auxiliary.h

@@ -38,7 +38,8 @@ class PassMethodsClass {
     static const int pm_kickmap_pass                   = 8;
     static const int pm_matrix_pass                    = 9;
     static const int pm_drift_g2l_pass                 = 10;
-    static const int pm_nr_pms                         = 11;  // counter for number of passmethods
+    static const int pm_field3d_pass                   = 11;
+    static const int pm_nr_pms                         = 12;  // counter for number of passmethods
     static const std::vector<int> time_aware_passmethods;
     PassMethodsClass() {
         passmethods.push_back("identity_pass");
@@ -52,6 +53,7 @@ class PassMethodsClass {
         passmethods.push_back("kicktable_pass");
         passmethods.push_back("matrix_pass");
         passmethods.push_back("drift_g2l_pass");
+        passmethods.push_back("field3d_pass");
     }
     int size() const { return passmethods.size(); }
     std::string operator[](const int i) const { return passmethods[i]; }
@@ -76,7 +78,8 @@ struct PassMethod {
         pm_kickmap_pass                   = 8,
         pm_matrix_pass                    = 9,
         pm_drift_g2l_pass                 = 10,
-        pm_nr_pms                         = 11,
+        pm_field3d_pass                   = 11,
+        pm_nr_pms                         = 12,
     };
 };
 
@@ -94,6 +97,7 @@ const std::vector<std::string> pm_dict = {
     "kicktable_pass",
     "matrix_pass",
     "drift_g2l_pass",
+    "field3d_pass",
 };
 
 struct RadiationState {

include/trackcpp/coefmatrix.h

@@ -0,0 +1,124 @@
+// TRACKCPP - Particle tracking code
+// Copyright (C) 2015  LNLS Accelerator Physics Group
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef COEFMATRIX_H
+#define COEFMATRIX_H
+
+#include <cstddef>
+#include <vector>
+
+#ifndef SWIG
+template<typename T>
+class MatrixRow {
+public:
+    explicit MatrixRow(T* ptr)
+        : ptr_(ptr)
+    {}
+
+    T& operator[](size_t j)
+    {
+        return ptr_[j];
+    }
+
+private:
+    T* ptr_;
+};
+#endif
+
+
+class CoefMatrix {
+public:
+
+#ifndef SWIG
+    using Row = MatrixRow<double>;
+    using ConstRow = MatrixRow<const double>;
+#endif
+
+    CoefMatrix() = default;
+
+    CoefMatrix(size_t rows, size_t cols)
+        : rows_(rows),
+          cols_(cols),
+          data_(rows * cols, 0.0)
+    {}
+
+    void resize(size_t rows, size_t cols)
+    {
+        rows_ = rows;
+        cols_ = cols;
+        data_.resize(rows * cols);
+    }
+
+    double& operator()(size_t i, size_t j)
+    {
+        return data_[i * cols_ + j];
+    }
+
+    const double& operator()(size_t i, size_t j) const
+    {
+        return data_[i * cols_ + j];
+    }
+
+#ifndef SWIG
+    Row operator[](size_t i)
+    {
+        return Row(data_.data() + i * cols_);
+    }
+
+    ConstRow operator[](size_t i) const
+    {
+        return ConstRow(data_.data() + i * cols_);
+    }
+#endif
+
+    size_t rows() const
+    {
+        return rows_;
+    }
+
+    size_t cols() const
+    {
+        return cols_;
+    }
+
+    size_t size() const
+    {
+        return data_.size();
+    }
+
+    bool empty() const
+    {
+        return data_.empty();
+    }
+
+    double* data()
+    {
+        return data_.data();
+    }
+
+    const double* data() const
+    {
+        return data_.data();
+    }
+
+private:
+
+    size_t rows_ = 0;
+    size_t cols_ = 0;
+    std::vector<double> data_;
+};
+
+#endif

include/trackcpp/elements.h

@@ -24,6 +24,7 @@
 #include <string>
 #include <fstream>
 #include <cfloat>
+#include "coefmatrix.h"
 
 
 class Element {
@@ -63,7 +64,12 @@ class Element {
   double        phase_lag   = 0;  // [rad]
   int           kicktable_idx = -1;   // index of kickmap object in kicktable_list
   double        rescale_kicks = 1.0;  // for kickmaps
-
+  double        field3d_hori_ks = 0;       // [1/m]
+  double        field3d_hori_kx = 0;       // [1/m]
+  double        field3d_hori_s0 = 0;       // [m]
+  CoefMatrix field3d_hori_coefs_cos;
+  CoefMatrix field3d_hori_coefs_sin;
+
   std::vector<double> polynom_a = default_polynom;
   std::vector<double> polynom_b = default_polynom;
   Matrix              matrix66 = Matrix(6);
@@ -111,6 +117,10 @@ class Element {
   static Element sextupole  (const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5);
   static Element rfcavity   (const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag_);
   static Element kickmap    (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length_ = 1.0, const double& rescale_kicks_ = 1.0);
+  static Element field3d    (const std::string& fam_name_, const double& length_, const double& field3d_hori_s0_, const double& field3d_hori_kx_, const double& field3d_hori_ks_,
+                             const CoefMatrix& field3d_hori_coefs_cos_,
+                             const CoefMatrix& field3d_hori_coefs_sin_,
+                             const int nr_steps_ = 40);
 
   bool operator==(const Element& o) const;
   bool operator!=(const Element& o) const { return !(*this == o); };
@@ -132,5 +142,9 @@ void initialize_quadrupole(Element& element, const double& K, const int& nr_step
 void initialize_sextupole(Element& element, const double& S, const int& nr_steps);
 void initialize_rfcavity(Element& element, const double& frequency, const double& voltage, const double& phase_lag);
 void initialize_kickmap(Element& element, const int& kicktable_idx, const int& nr_steps, const double &rescale_kicks);
+void initialize_field3d(Element& element, const double& field3d_hori_s0_, const double& field3d_hori_kx_, const double& field3d_hori_ks_,
+                        const CoefMatrix& field3d_hori_coefs_cos_,
+                        const CoefMatrix& field3d_hori_coefs_sin_,
+                        const int nr_steps_);
 
 #endif