Add IDs Passmethod

INSTALL.md

@@ -19,7 +19,7 @@ INSTRUCTIONS
 
 REQUIREMENTS
 
-- a C++ compiler with C++11 support (tested with GCC 4.8)
+- a C++ compiler with C++17 support (tested with GCC 4.8)
 - pthread
 - blas
 - gsl (GNU Scientific Library)
@@ -49,7 +49,7 @@ INSTRUCTIONS
 
 REQUIREMENTS
 
-- a C++ compiler with C++11 support (tested with GCC 4.8)
+- a C++ compiler with C++17 support (tested with GCC 4.8)
 - blas
 - gsl (GNU Scientific Library)
 - SWIG 2 (tested with SWIG 2.0.11)

Makefile

@@ -36,8 +36,8 @@ CC          = gcc
 CXX         = g++
 AR          = ar
 MACHINE     = -m64
-OPT_FLAG    = -O3 -std=c++11 -fPIC
-DBG_FLAG    = -O0 -g3 -std=c++11 -fPIC
+OPT_FLAG    = -O3 -std=c++17 -fPIC
+DBG_FLAG    = -O0 -g3 -std=c++17 -fPIC
 ARFLAGS     = rcs
 DFLAGS      = -DVERSION=$(VERSION)
 LIBSOURCES_CPP  =	lattice.cpp \
@@ -58,7 +58,7 @@ LIBSOURCES_CPP  =	lattice.cpp \
 					auxiliary.cpp
 BINSOURCES_CPP =	exec.cpp \
 					tests.cpp \
-					commands.cpp \
+					commands.cpp
 
 
 AUXFILES  = VERSION

include/trackcpp/auxiliary.h

@@ -37,7 +37,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_kickpoly_pass                  = 11;
+    static const int pm_nr_pms                         = 12;  // counter for number of passmethods
     PassMethodsClass() {
         passmethods.push_back("identity_pass");
         passmethods.push_back("drift_pass");
@@ -50,6 +51,7 @@ class PassMethodsClass {
         passmethods.push_back("kicktable_pass");
         passmethods.push_back("matrix_pass");
         passmethods.push_back("drift_g2l_pass");
+        passmethods.push_back("kickpoly_pass");
     }
     int size() const { return passmethods.size(); }
     std::string operator[](const int i) const { return passmethods[i]; }
@@ -71,7 +73,8 @@ struct PassMethod {
         pm_kickmap_pass                   = 8,
         pm_matrix_pass                    = 9,
         pm_drift_g2l_pass                 = 10,
-        pm_nr_pms                         = 11,
+        pm_kickpoly_pass                  = 11,
+        pm_nr_pms                         = 12,
     };
 };
 
@@ -89,6 +92,7 @@ const std::vector<std::string> pm_dict = {
     "kicktable_pass",
     "matrix_pass",
     "drift_g2l_pass",
+    "kickpoly_pass",
 };
 
 struct RadiationState {

include/trackcpp/elements.h

@@ -61,11 +61,13 @@ 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
+  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;
+  std::vector<double> polynom_kickx = default_polynom;  // for kickpoly
+  std::vector<double> polynom_kicky = default_polynom;  // for kickpoly
   Matrix              matrix66 = Matrix(6);
 
   double t_in[6];
@@ -111,7 +113,7 @@ 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 kickpoly   (const std::string& fam_name_, const double& length_, const int nr_steps_ = 20, const double& rescale_kicks_ = 1.0);
   bool operator==(const Element& o) const;
   bool operator!=(const Element& o) const { return !(*this == o); };
 
@@ -132,5 +134,6 @@ 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_kickpoly(Element& element, const int& nr_steps, const double &rescale_kicks);
 
 #endif

include/trackcpp/flat_file.h

@@ -38,8 +38,6 @@ struct FlatFileType {
     };
   };
 
-//Status::type read_flat_file(const std::string& filename, Accelerator& accelerator);
-
 Status::type read_flat_file(std::string& filename, Accelerator& accelerator, bool file_flag = true);
 Status::type write_flat_file(std::string& filename, const Accelerator& accelerator, bool file_flag = true);
 

include/trackcpp/kicktable.h

@@ -18,108 +18,159 @@
 #define _KICKTABLE_H
 
 #include "auxiliary.h"
+// #include "../alglib/interpolation.h"
 #include <string>
 #include <vector>
+#include <list>
+#include <algorithm>
 
 
 class Kicktable {
 
-  // Kicktable: assumes x_kick and y_kick tables sampled on the same regular (x,y)-point grid.
-
 public:
 
   std::string         filename;
   double              length;
-  unsigned int        x_nrpts, y_nrpts;
-  double              x_min, x_max;
-  double              y_min, y_max;
+  std::vector<double> x_pos, y_pos;
   std::vector<double> x_kick,  y_kick;
-
+  static std::list<Kicktable> kicktable_list;
+
+  Kicktable(
+    const std::vector<double>& x_pos,
+    const std::vector<double>& y_pos,
+    const std::vector<double>& x_kick,
+    const std::vector<double>& y_kick,
+    const double length = 1
+  );
   Kicktable(const std::string& filename_ = "");
   Kicktable(const Kicktable &) = default;
 
-  Status::type load_from_file(const std::string& filename_);
-
-  unsigned int get_idx(unsigned int ix, unsigned int iy) const { return iy*x_nrpts+ix; }
-  double       get_x(unsigned int ix) const { return x_min + ix * (x_max - x_min) / (x_nrpts - 1.0); }
-  double       get_y(unsigned int iy) const { return y_min + iy * (y_max - y_min) / (y_nrpts - 1.0); }
-  template <typename T> unsigned int get_ix(const T& x) const { return (int) ((x - x_min) / ((x_max - x_min) / (x_nrpts - 1))); }
-  template <typename T> unsigned int get_iy(const T& y) const { return (int) ((y - y_min) / ((y_max - y_min) / (y_nrpts - 1))); }
-
-  bool operator==(const Kicktable& o) const;
-  bool operator!=(const Kicktable& o) const { return !(*this == o); }
-
-};
-
-
-extern std::vector<Kicktable> kicktable_list;
-
-
-int add_kicktable(const std::string& filename);
-void clear_kicktables(std::vector<Kicktable>& kicktable_list);
-
-template <typename T>
-Status::type kicktable_getkicks_bilinear(const int& kicktable_idx, const T& rx, const T& ry, T& hkick, T& vkick) {
-
-  const Kicktable* kicktable = &kicktable_list[kicktable_idx];
+  Status::type load_from_file(
+    const std::string& filename_, bool file_flag=true
+  );
+  Status::type _dump_to_stream(
+    std::ostream& fp, const std::string author_name = " "
+  );
+  Status::type save_to_file(
+    const std::string filename_, const std::string author_name = " "
+  );
+  std::string save_to_string(const std::string author_name = " ");
+
+  bool is_valid_kicktable() const;
+  unsigned int get_idx(unsigned int ix, unsigned int iy) const
+  {
+    return iy * x_pos.size() + ix;
+  }
+  double get_x(unsigned int ix) const {return x_pos[ix];}
+  double get_y(unsigned int iy) const {return y_pos[iy];}
+  template <typename T> unsigned int get_ix(const T& x) const
+  {
+    int idx;
+    idx = std::lower_bound(x_pos.begin(), x_pos.end(), x) - x_pos.begin() - 1;
+    return idx >= 0 ? (unsigned int) idx : 0u;
+  }
+  template <typename T> unsigned int get_iy(const T& y) const
+  {
+    int idx;
+    idx = std::lower_bound(y_pos.begin(), y_pos.end(), y) - y_pos.begin() - 1;
+    return idx >= 0 ? (unsigned int) idx : 0u;
+  }
+  unsigned int get_ix(const double& x) const {return get_ix<double>(x);}
+  unsigned int get_iy(const double& y) const {return get_iy<double>(y);}
+
+  template <typename T> Status::type getkicks_bilinear(
+    const T& rx, const T& ry, T& hkick, T& vkick
+  ) const
+  {
+
+    // gets indices
+    unsigned int ix = get_ix(rx);
+    unsigned int iy = get_iy(ry);
+    unsigned int ixp1, iyp1;
+    if (ix >= x_pos.size()-1) ixp1 = ix--;
+    else ixp1 = ix + 1;
+    if (iy >= y_pos.size()-1) iyp1 = iy--;
+    else iyp1 = iy + 1;
+    const unsigned int i00 = get_idx(ix, iy);
+    const unsigned int i01 = get_idx(ix, iyp1);
+    const unsigned int i10 = get_idx(ixp1, iy);
+    const unsigned int i11 = get_idx(ixp1, iyp1);
+
+    /* coordinates */
+    const double x1 = get_x(ix);
+    const double x2 = get_x(ixp1);
+    const double y1 = get_y(iy);
+    const double y2 = get_y(iyp1);
+    const double denom = (x2 - x1) * (y2 - y1);
+    const T dx1 = rx - x1;
+    const T dx2 = x2 - rx;
+    const T dy1 = ry - y1;
+    const T dy2 = y2 - ry;
+
+    // Bilinear interpolation
+    // https://en.wikipedia.org/wiki/Bilinear_interpolation
+    {  /* hkick */
+      const double& f11 = x_kick[i00];
+      const double& f12 = x_kick[i01];
+      const double& f21 = x_kick[i10];
+      const double& f22 = x_kick[i11];
+      hkick = (f11 * dx2 + f21 * dx1) * dy2 + (f12 * dx2 + f22 * dx1) * dy1;
+      hkick /= denom;
+    }
+    {  /* vkick */
+      const double& f11 = y_kick[i00];
+      const double& f12 = y_kick[i01];
+      const double& f21 = y_kick[i10];
+      const double& f22 = y_kick[i11];
+      vkick = (f11 * dx2 + f21 * dx1) * dy2 + (f12 * dx2 + f22 * dx1) * dy1;
+      vkick /= denom;
+    }
+
+    return Status::success;
+  }
 
-  // checks x limits
-  const double& xmin = kicktable->x_min;
-  const double& xmax = kicktable->x_max;
-  if ((rx < xmin) or (rx > xmax)) {
-    //std::cout << "rx: " << double(rx) << std::endl;
-    hkick = nan("");
-    return Status::kicktable_out_of_range;
+  template <typename T> Status::type getkicks(
+    const T& rx, const T& ry, T& hkick, T& vkick
+  ) const
+  {
+    return getkicks_bilinear(rx, ry, hkick, vkick);
   }
 
-  // checks y limits
-  const double& ymin = kicktable->y_min;
-  const double& ymax = kicktable->y_max;
-  if ((ry < ymin) or (ry > ymax)) {
-    //std::cout << "ry: " << double(ry) << std::endl;
-    vkick = nan("");
-    return Status::kicktable_out_of_range;
+  Status::type getkicks(
+    const double& rx, const double& ry, double& hkick__, double& vkick__
+  ) const
+  {
+    return getkicks<double>(rx, ry, hkick__, vkick__);
   }
 
-  // gets indices
-  const unsigned int ix = kicktable->get_ix(rx);
-  const unsigned int iy = kicktable->get_iy(ry);
-
-  /* coordinates */
-  const double x1  = kicktable->get_x(ix);
-  const double x2  = kicktable->get_x(ix+1);
-  const double y1  = kicktable->get_y(iy);
-  const double y2  = kicktable->get_y(iy+1);
-
-  // Bilinear interpolation - https://en.wikipedia.org/wiki/Bilinear_interpolation
-
-  {  /* hkick */
-    const double fq11 = kicktable->x_kick[kicktable->get_idx(ix+0, iy+0)];
-    const double fq12 = kicktable->x_kick[kicktable->get_idx(ix+0, iy+1)];
-    const double fq21 = kicktable->x_kick[kicktable->get_idx(ix+1, iy+0)];
-    const double fq22 = kicktable->x_kick[kicktable->get_idx(ix+1, iy+1)];
-    const T fR1 = ((x2 - rx) * fq11 + (rx - x1) * fq21) / (x2 - x1);
-    const T fR2 = ((x2 - rx) * fq12 + (rx - x1) * fq22) / (x2 - x1);
-    const T fP  = ((y2 - ry) * fR1  + (ry - y1) * fR2 ) / (y2 - y1);
-    hkick = fP;
+  static int add_kicktable(
+    const std::vector<double>& x_pos,
+    const std::vector<double>& y_pos,
+    const std::vector<double>& x_kick,
+    const std::vector<double>& y_kick,
+    const double length=1
+  );
+  static int add_kicktable(const std::string filename);
+  static int add_kicktable(const Kicktable &new_ktab);
+  static Kicktable& get_kicktable(const int& kicktable_idx)
+  {
+    if (not is_valid_kicktable_index(kicktable_idx))
+      throw std::out_of_range("kicktable_idx is out of range.");
+
+    auto it = kicktable_list.begin();
+    std::advance(it, kicktable_idx);
+    return *it;
   }
-  {  /* vkick */
-    const double fq11 = kicktable->y_kick[kicktable->get_idx(ix+0, iy+0)];
-    const double fq12 = kicktable->y_kick[kicktable->get_idx(ix+0, iy+1)];
-    const double fq21 = kicktable->y_kick[kicktable->get_idx(ix+1, iy+0)];
-    const double fq22 = kicktable->y_kick[kicktable->get_idx(ix+1, iy+1)];
-    const T fR1 = ((x2 - rx) * fq11 + (rx - x1) * fq21) / (x2 - x1);
-    const T fR2 = ((x2 - rx) * fq12 + (rx - x1) * fq22) / (x2 - x1);
-    const T fP  = ((y2 - ry) * fR1  + (ry - y1) * fR2 ) / (y2 - y1);
-    vkick = fP;
+  static bool is_valid_kicktable_index(const int idx)
+  {
+    return ((idx>=0) & (idx < kicktable_list.size()));
   }
+  static void clear_kicktables();
+  static size_t get_kicktable_list_size(){return kicktable_list.size();}
 
-  return Status::success;
-}
+  bool operator==(const Kicktable& o) const;
+  bool operator!=(const Kicktable& o) const { return !(*this == o); }
 
-template <typename T>
-Status::type kicktable_getkicks(const int& kicktable_idx, const T& rx, const T& ry, T& hkick, T& vkick) {
-  return kicktable_getkicks_bilinear(kicktable_idx, rx, ry, hkick, vkick);
-}
+};
 
 #endif

include/trackcpp/lattice.h

@@ -37,6 +37,8 @@ std::vector<int>     latt_findcells_angle       (const std::vector<Element>& lat
 std::vector<int>     latt_findcells_frequency   (const std::vector<Element>& lattice, const double& value,      bool reverse = false);
 std::vector<int>     latt_findcells_polynom_b   (const std::vector<Element>& lattice, unsigned int n, const double& value, bool reverse = false);
 std::vector<int>     latt_findcells_polynom_a   (const std::vector<Element>& lattice, unsigned int n, const double& value, bool reverse = false);
+std::vector<int>     latt_findcells_polynom_kickx(const std::vector<Element>& lattice, unsigned int n, const double& value, bool reverse = false);
+std::vector<int>     latt_findcells_polynom_kicky(const std::vector<Element>& lattice, unsigned int n, const double& value, bool reverse = false);
 std::vector<int>     latt_findcells_pass_method (const std::vector<Element>& lattice, const std::string& value, bool reverse = false);
 
 template <typename T>

include/trackcpp/passmethods.h

@@ -66,6 +66,7 @@ template <typename T> Status::type pm_thinsext_pass              (Pos<T> &pos, c
 template <typename T> Status::type pm_kickmap_pass               (Pos<T> &pos, const Element &elem, const Accelerator& accelerator);
 template <typename T> Status::type pm_matrix_pass                (Pos<T> &pos, const Element &elem, const Accelerator& accelerator);
 template <typename T> Status::type pm_drift_g2l_pass             (Pos<T> &pos, const Element &elem, const Accelerator& accelerator);
+template <typename T> Status::type pm_kickpoly_pass              (Pos<T> &pos, const Element &elem, const Accelerator& accelerator);
 
 #include "passmethods.hpp"