Add IDs Passmethod

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

@@ -66,6 +66,8 @@ class Element {
 
   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/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"
 

include/trackcpp/passmethods.hpp

@@ -95,13 +95,32 @@ T b2_perp(const T& bx, const T& by, const T& px, const T& py, const T& curv=1) {
 }
 
 template <typename T>
-Status::type kicktablethinkick(Pos<T>& pos, const int& kicktable_idx,
-                               const double& brho, const int nr_steps, const double& rescale_kicks) {
-
+Status::type kicktablethinkick(
+  Pos<T>& pos,
+  const int& kicktable_idx,
+  const double& brho2,
+  const int nr_steps,
+  const double& rescale_kicks
+)
+{
   T hkick, vkick;
   Status::type status = kicktable_getkicks(kicktable_idx, pos.rx, pos.ry, hkick, vkick);
-  pos.px += rescale_kicks * hkick / (brho * brho) / nr_steps;
-  pos.py += rescale_kicks * vkick / (brho * brho) / nr_steps;
+  // According to Ellaune's theory of kick maps:
+  // https://accelconf.web.cern.ch/e92/PDF/EPAC1992_0661.PDF
+  //
+  // there should be a dependency of the kicks with 1/(1+delta)² for kick maps
+  // generated with the potential function calculated from the field
+  // integrals. However, most kicks maps we use nowadays comes from
+  // Runge-Kutta integration of the equations of motion. Tests with different
+  // energies with this integration setup have shown a very complicated
+  // behavior of the kicks as function of energy, corroborating with no of
+  // very weak dependence. For this reason the kicks here are influenced by
+  // the energy of each particle. Looking at the code of AT and Elegant, I
+  // noticed that AT normalizes the kicks by 1/(1+delta) while Elegant
+  // normalizes by 1/(1+delta)². None of them, however, adds the necessary
+  // terms to dl to make the map symplectic.
+  pos.px += rescale_kicks * hkick / brho2 / nr_steps;
+  pos.py += rescale_kicks * vkick / brho2 / nr_steps;
   if (status == Status::kicktable_out_of_range) {
     if (not isfinite(pos.px)) {
       pos.rx = nan("");
@@ -113,6 +132,52 @@ Status::type kicktablethinkick(Pos<T>& pos, const int& kicktable_idx,
   return status;
 }
 
+template <typename T>
+void kickpolythinkick(
+  Pos<T>& pos,
+  const std::vector<double> polyx,
+  const std::vector<double> polyy,
+  const double rescale,
+  const double& brho2,
+  const unsigned int nr_steps,
+  const unsigned int order
+)
+{
+  T hkick = 0;
+  T vkick = 0;
+  T rx_p = 1;
+  unsigned int k = 0;
+  for (auto i = 0; i <= order; ++i)
+  {
+    T ry_p = 1;
+    for (auto j = 0; j <= order; ++j)
+    {
+      if (i + j > order) break;
+      hkick += polyx[k] * rx_p * ry_p;
+      vkick += polyy[k] * rx_p * ry_p;
+      ry_p *= pos.ry;
+      ++k;
+    }
+    rx_p *= pos.rx;
+  }
+  // According to Ellaune's theory of kick maps:
+  // https://accelconf.web.cern.ch/e92/PDF/EPAC1992_0661.PDF
+  //
+  // there should be a dependency of the kicks with 1/(1+delta)² for kick maps
+  // generated with the potential function calculated from the field
+  // integrals. However, most kicks maps we use nowadays comes from
+  // Runge-Kutta integration of the equations of motion. Tests with different
+  // energies with this integration setup have shown a very complicated
+  // behavior of the kicks as function of energy, corroborating with no of
+  // very weak dependence. For this reason the kicks here are influenced by
+  // the energy of each particle. Looking at the code of AT and Elegant, I
+  // noticed that AT normalizes the kicks by 1/(1+delta) while Elegant
+  // normalizes by 1/(1+delta)². None of them, however, adds the necessary
+  // terms to dl to make the map symplectic.
+  pos.px += rescale * hkick / brho2 / nr_steps;
+  pos.py += rescale * vkick / brho2 / nr_steps;
+}
+
 template <typename T>
 void matthinkick(Pos<T> &pos, const Matrix &m) {
 
@@ -450,15 +515,22 @@ Status::type pm_kickmap_pass(Pos<T> &pos, const Element &elem,
   Status::type status = Status::success;
 
   double sl   = elem.length / float(elem.nr_steps);
-  const double brho = get_magnetic_rigidity(accelerator.energy);
+  double brho2 = get_magnetic_rigidity(accelerator.energy);
+  brho2 *= brho2;
 
   global_2_local(pos, elem);
   if (elem.kicktable_idx < 0) {
     drift<T>(pos, sl);
   } else {
     for(unsigned int i=0; i<elem.nr_steps; ++i) {
       drift<T>(pos, sl / 2);
-      Status::type status = kicktablethinkick(pos, elem.kicktable_idx, brho, elem.nr_steps, elem.rescale_kicks);
+      Status::type status = kicktablethinkick(
+        pos,
+        elem.kicktable_idx,
+        brho2,
+        elem.nr_steps,
+        elem.rescale_kicks
+      );
       if (status != Status::success) return status;
       drift<T>(pos, sl / 2);
     }
@@ -468,6 +540,44 @@ Status::type pm_kickmap_pass(Pos<T> &pos, const Element &elem,
   return status;
 }
 
+template <typename T>
+Status::type pm_kickpoly_pass(
+  Pos<T> &pos,
+  const Element &elem,
+  const Accelerator& accelerator
+){
+  Status::type status = Status::success;
+
+  double sl   = elem.length / float(elem.nr_steps);
+  double brho2 = get_magnetic_rigidity(accelerator.energy);
+  brho2 *= brho2;
+
+  // the polynom_kickx and polynom_kicky variables are the coefficients
+  // of 2D polynoms ordered in the following way:
+  // f(x,y) = a0 + a1*x + a2*y + a3*x² + a4*x*y + a5*y² + a6*x³ + a7*x²*y +...
+  //
+  // Find the order of the polynom, given its length:
+  const unsigned int order = (sqrt(1 + 8*elem.polynom_kickx.size()) - 1) / 2;
+
+  global_2_local(pos, elem);
+  for(unsigned int i=0; i<elem.nr_steps; ++i) {
+    drift<T>(pos, sl / 2);
+    kickpolythinkick(
+      pos,
+      elem.polynom_kickx,
+      elem.polynom_kicky,
+      elem.rescale_kicks,
+      brho2,
+      elem.nr_steps,
+      order
+    );
+    drift<T>(pos, sl / 2);
+  }
+  local_2_global(pos, elem);
+
+  return status;
+}
+
 template <typename T>
 Status::type pm_matrix_pass(Pos<T> &pos, const Element &elem,
                             const Accelerator& accelerator) {

include/trackcpp/tracking.h

@@ -94,6 +94,9 @@ Status::type track_elementpass (
     case PassMethod::pm_drift_g2l_pass:
         if ((status = pm_drift_g2l_pass<T>(orig_pos, el, accelerator)) != Status::success) return status;
         break;
+    case PassMethod::pm_kickpoly_pass:
+        if ((status = pm_kickpoly_pass<T>(orig_pos, el, accelerator)) != Status::success) return status;
+        break;
     default:
         return Status::passmethod_not_defined;
     }

python_package/interface.cpp

@@ -237,6 +237,10 @@ Element kickmap_wrapper(const std::string& fam_name_,  const std::string& kickta
     return Element::kickmap(fam_name_, kicktable_fname_, nr_steps_, rescale_length_, rescale_kicks_);
 }
 
+Element kickpoly_wrapper(const std::string& fam_name_, const double& length_, const int nr_steps_, const double& rescale_kicks_) {
+    return Element::kickpoly(fam_name_, length_, nr_steps_, rescale_kicks_);
+}
+
 Status::type read_flat_file_wrapper(String& fname, Accelerator& accelerator, bool file_flag) {
   return read_flat_file(fname.data, accelerator, file_flag);
 }

python_package/interface.h

@@ -91,6 +91,7 @@ Element quadrupole_wrapper(const std::string& fam_name_, const double& length_,
 Element sextupole_wrapper(const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5);
 Element rfcavity_wrapper(const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag);
 Element kickmap_wrapper(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);
+Element kickpoly_wrapper(const std::string& fam_name_, const double& length, const int nr_steps_ = 20, const double& rescale_kicks = 1.0);
 Element rbend_wrapper(const std::string& fam_name_, const double& length_,
                       const double& angle_, const double& angle_in_, const double& angle_out_,
                       const double& gap_, const double& fint_in_, const double& fint_out_,

src/elements.cpp

@@ -152,10 +152,15 @@ Element Element::kickmap (const std::string& fam_name_, const std::string& kickt
 
   const Kicktable& kicktable = kicktable_list[idx];
   Element e = Element(fam_name_, rescale_length_ * kicktable.length);
-    initialize_kickmap(e, idx, nr_steps_, rescale_kicks_);
+  initialize_kickmap(e, idx, nr_steps_, rescale_kicks_);
   return e;
 }
 
+Element Element::kickpoly (const std::string& fam_name_, const double& length_, const int nr_steps_, const double &rescale_kicks_) {
+  Element e = Element(fam_name_, length_);
+  initialize_kickpoly(e, nr_steps_, rescale_kicks_);
+  return e;
+}
 
 void print_polynom(std::ostream& out, const std::string& label, const std::vector<double>& polynom) {
   int order = 0;
@@ -201,6 +206,8 @@ bool Element::operator==(const Element& o) const {
     if (this->phase_lag != o.phase_lag) return false;
     if (this->polynom_a != o.polynom_a) return false;
     if (this->polynom_b != o.polynom_b) return false;
+    if (this->polynom_kickx != o.polynom_kickx) return false;
+    if (this->polynom_kicky != o.polynom_kicky) return false;
     const Matrix& m = this->matrix66;
     const Matrix& mo = o.matrix66;
     for(unsigned int i=0; i<m.size(); ++i){
@@ -240,6 +247,8 @@ std::ostream& operator<< (std::ostream &out, const Element& el) {
   }
   print_polynom(        out,                "polynom_a     : ", el.polynom_a);
   print_polynom(        out,                "polynom_b     : ", el.polynom_b);
+  print_polynom(        out,                "polynom_kickx : ", el.polynom_kickx);
+  print_polynom(        out,                "polynom_kicky : ", el.polynom_kicky);
   if (el.frequency != 0)out << std::endl << "frequency     : " << el.frequency;
   if (el.voltage != 0)  out << std::endl << "voltage       : " << el.voltage;
   if (el.phase_lag != 0)out << std::endl << "phase_lag     : " << el.phase_lag;
@@ -308,3 +317,12 @@ void initialize_kickmap(Element& element, const int& kicktable_idx, const int& n
     element.kicktable_idx = kicktable_idx;
     element.rescale_kicks = rescale_kicks;
 }
+
+void initialize_kickpoly(
+  Element& element, const int& nr_steps, const double &rescale_kicks
+)
+{
+    element.pass_method = PassMethod::pm_kickpoly_pass;
+    element.nr_steps = nr_steps;
+    element.rescale_kicks = rescale_kicks;
+}

src/flat_file.cpp

@@ -36,6 +36,7 @@ static void write_6d_vector(std::ostream& fp, const std::string& label, const do
 static void write_6d_vector(std::ostream& fp, const std::string& label, const std::vector<double>& t);
 static void write_polynom(std::ostream& fp, const std::string& label, const std::vector<double>& p);
 static void synchronize_polynomials(Element& e);
+static void synchronize_polynomials_kick(Element& e);
 static void read_polynomials(std::ifstream& fp, Element& e);
 static void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator);
 static Status::type read_flat_file_trackcpp(std::istream&, Accelerator& accelerator);
@@ -112,6 +113,8 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator)
     }
     if (has_polynom(e.polynom_a)) write_polynom(fp, "polynom_a", e.polynom_a);
     if (has_polynom(e.polynom_b)) write_polynom(fp, "polynom_b", e.polynom_b);
+    if (has_polynom(e.polynom_kickx)) write_polynom(fp, "polynom_kickx", e.polynom_kickx);
+    if (has_polynom(e.polynom_kicky)) write_polynom(fp, "polynom_kicky", e.polynom_kicky);
     if (e.vchamber != VChamberShape::rectangle) { fp << std::setw(pw) << "vchamber" << e.vchamber << '\n'; }
     if (e.hmin != -DBL_MAX) { fp << std::setw(pw) << "hmin" << e.hmin << '\n'; }
     if (e.hmax != DBL_MAX) { fp << std::setw(pw) << "hmax" << e.hmax << '\n'; }
@@ -305,6 +308,40 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator)
       synchronize_polynomials(e);
       continue;
     }
+    if (cmd.compare("polynom_kickx") == 0) {
+      std::vector<unsigned int> order;
+      std::vector<double> multipole;
+      unsigned int size = 0;
+      while (not ss.eof()) {
+        unsigned int o; double m; ss >> o >> m;
+        if (ss.eof()) break;
+        order.push_back(o); multipole.push_back(m);
+        if (o+1 > size) size = o+1;
+      }
+      if (size > 0) {
+        e.polynom_kickx.resize(size, 0);
+        for(unsigned int i=0; i<order.size(); ++i) e.polynom_kickx[order[i]] = multipole[i];
+      }
+      synchronize_polynomials_kick(e);
+      continue;
+    }
+    if (cmd.compare("polynom_kicky") == 0) {
+      std::vector<unsigned int> order;
+      std::vector<double> multipole;
+      unsigned int size = 0;
+      while (not ss.eof()) {
+        unsigned int o; double m; ss >> o >> m;
+        if (ss.eof()) break;
+        order.push_back(o); multipole.push_back(m);
+        if (o+1 > size) size = o+1;
+      }
+      if (size > 0) {
+        e.polynom_kicky.resize(size, 0);
+        for(unsigned int i=0; i<order.size(); ++i) e.polynom_kicky[order[i]] = multipole[i];
+      }
+      synchronize_polynomials_kick(e);
+      continue;
+    }
     if (line.size()<2) continue;
     return Status::flat_file_error;
   }
@@ -433,6 +470,15 @@ static void synchronize_polynomials(Element& e) {
 
 }
 
+static void synchronize_polynomials_kick(Element& e) {
+
+  unsigned int size = (e.polynom_kickx.size() > e.polynom_kicky.size()) ? e.polynom_kickx.size() : e.polynom_kicky.size();
+  e.polynom_kickx.resize(size, 0);
+  e.polynom_kicky.resize(size, 0);
+
+}
+
+//This method is only used to read tracy3 files.
 static void read_polynomials(std::ifstream& fp, Element& e) {
   unsigned int nr_monomials, n_design, order;
   e.polynom_a = std::vector<double>(Element::default_polynom);