From 9b16579eaef0593176ae5ba7d61d0e052df78110 Mon Sep 17 00:00:00 2001 From: Fernando Date: Fri, 16 Aug 2024 16:30:50 -0300 Subject: [PATCH 01/28] ENH: Add pm_kickpoly_pass to elements.h and elements.cpp. --- include/trackcpp/elements.h | 5 ++++- src/elements.cpp | 20 +++++++++++++++++++- 2 files changed, 23 insertions(+), 2 deletions(-) diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 94942d9..56dee27 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -66,6 +66,8 @@ class Element { std::vector polynom_a = default_polynom; std::vector polynom_b = default_polynom; + std::vector polynom_kickx = default_polynom; // for kickpoly + std::vector 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 diff --git a/src/elements.cpp b/src/elements.cpp index 5d96fb1..7a75f35 100644 --- a/src/elements.cpp +++ b/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& 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 Date: Fri, 16 Aug 2024 16:33:00 -0300 Subject: [PATCH 02/28] ENH: Add pm_kickpoly_pass to passmethods.h and passmethods.cpp. --- include/trackcpp/passmethods.h | 1 + include/trackcpp/passmethods.hpp | 124 +++++++++++++++++++++++++++++-- 2 files changed, 118 insertions(+), 7 deletions(-) diff --git a/include/trackcpp/passmethods.h b/include/trackcpp/passmethods.h index 51a38aa..b5877b6 100644 --- a/include/trackcpp/passmethods.h +++ b/include/trackcpp/passmethods.h @@ -66,6 +66,7 @@ template Status::type pm_thinsext_pass (Pos &pos, c template Status::type pm_kickmap_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_matrix_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_drift_g2l_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_kickpoly_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); #include "passmethods.hpp" diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 9b315ea..2b38665 100644 --- a/include/trackcpp/passmethods.hpp +++ b/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 -Status::type kicktablethinkick(Pos& pos, const int& kicktable_idx, - const double& brho, const int nr_steps, const double& rescale_kicks) { - +Status::type kicktablethinkick( + Pos& 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& pos, const int& kicktable_idx, return status; } +template +void kickpolythinkick( + Pos& pos, + const std::vector polyx, + const std::vector 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 void matthinkick(Pos &pos, const Matrix &m) { @@ -450,7 +515,8 @@ Status::type pm_kickmap_pass(Pos &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) { @@ -458,7 +524,13 @@ Status::type pm_kickmap_pass(Pos &pos, const Element &elem, } else { for(unsigned int i=0; i(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(pos, sl / 2); } @@ -468,6 +540,44 @@ Status::type pm_kickmap_pass(Pos &pos, const Element &elem, return status; } +template +Status::type pm_kickpoly_pass( + Pos &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(pos, sl / 2); + kickpolythinkick( + pos, + elem.polynom_kickx, + elem.polynom_kicky, + elem.rescale_kicks, + brho2, + elem.nr_steps, + order + ); + drift(pos, sl / 2); + } + local_2_global(pos, elem); + + return status; +} + template Status::type pm_matrix_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { From 1de5d9f1e2741a0cab5ca406cf8cb348ff69984b Mon Sep 17 00:00:00 2001 From: Fernando Date: Fri, 16 Aug 2024 16:34:46 -0300 Subject: [PATCH 03/28] ENH: Add pm_kickpoly_pass to auxiliary.h, tracking.h and flat_file.cpp. --- include/trackcpp/auxiliary.h | 8 +++++-- include/trackcpp/tracking.h | 3 +++ src/flat_file.cpp | 46 ++++++++++++++++++++++++++++++++++++ 3 files changed, 55 insertions(+), 2 deletions(-) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 3cdfdc5..96804a7 100644 --- a/include/trackcpp/auxiliary.h +++ b/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 pm_dict = { "kicktable_pass", "matrix_pass", "drift_g2l_pass", + "kickpoly_pass", }; struct RadiationState { diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index a05c7a4..9735055 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -94,6 +94,9 @@ Status::type track_elementpass ( case PassMethod::pm_drift_g2l_pass: if ((status = pm_drift_g2l_pass(orig_pos, el, accelerator)) != Status::success) return status; break; + case PassMethod::pm_kickpoly_pass: + if ((status = pm_kickpoly_pass(orig_pos, el, accelerator)) != Status::success) return status; + break; default: return Status::passmethod_not_defined; } diff --git a/src/flat_file.cpp b/src/flat_file.cpp index 52ff4f3..d241409 100644 --- a/src/flat_file.cpp +++ b/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& t); static void write_polynom(std::ostream& fp, const std::string& label, const std::vector& 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 order; + std::vector 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; + std::vector 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 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(Element::default_polynom); From 87c436d6c85ef5de0747048554b78720f8c5c65b Mon Sep 17 00:00:00 2001 From: Fernando Date: Fri, 16 Aug 2024 16:35:34 -0300 Subject: [PATCH 04/28] ENH: Add pm_kickpoly_pass to lattice.h and lattice.cpp. --- include/trackcpp/lattice.h | 2 ++ src/lattice.cpp | 20 ++++++++++++++++++++ 2 files changed, 22 insertions(+) diff --git a/include/trackcpp/lattice.h b/include/trackcpp/lattice.h index 97ef182..44e0d75 100644 --- a/include/trackcpp/lattice.h +++ b/include/trackcpp/lattice.h @@ -37,6 +37,8 @@ std::vector latt_findcells_angle (const std::vector& lat std::vector latt_findcells_frequency (const std::vector& lattice, const double& value, bool reverse = false); std::vector latt_findcells_polynom_b (const std::vector& lattice, unsigned int n, const double& value, bool reverse = false); std::vector latt_findcells_polynom_a (const std::vector& lattice, unsigned int n, const double& value, bool reverse = false); +std::vector latt_findcells_polynom_kickx(const std::vector& lattice, unsigned int n, const double& value, bool reverse = false); +std::vector latt_findcells_polynom_kicky(const std::vector& lattice, unsigned int n, const double& value, bool reverse = false); std::vector latt_findcells_pass_method (const std::vector& lattice, const std::string& value, bool reverse = false); template diff --git a/src/lattice.cpp b/src/lattice.cpp index 910a841..b393e07 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -140,6 +140,26 @@ std::vector latt_findcells_polynom_a(const std::vector& lattice, u return r; } +std::vector latt_findcells_polynom_kickx(const std::vector& lattice, unsigned int n, const double& value, bool reverse) { + std::vector r; + for(unsigned int i=0; i latt_findcells_polynom_kicky(const std::vector& lattice, unsigned int n, const double& value, bool reverse) { + std::vector r; + for(unsigned int i=0; i latt_findcells_pass_method(const std::vector& lattice, const std::string& value, bool reverse) { std::vector r; for(unsigned int i=0; i Date: Fri, 16 Aug 2024 16:36:03 -0300 Subject: [PATCH 05/28] ENH: Add pm_kickpoly_pass to interface.h and interface.cpp. --- python_package/interface.cpp | 4 ++++ python_package/interface.h | 1 + 2 files changed, 5 insertions(+) diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 383c7fe..57d435a 100644 --- a/python_package/interface.cpp +++ b/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); } diff --git a/python_package/interface.h b/python_package/interface.h index e35e092..517f26a 100644 --- a/python_package/interface.h +++ b/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_, From 54335b1c61fa1b0b38ff1a4927a3740513a5ac46 Mon Sep 17 00:00:00 2001 From: Fernando Date: Wed, 21 Aug 2024 07:45:23 -0300 Subject: [PATCH 06/28] KICKTAB.ENH: Add nrpts in equality comparison. --- src/kicktable.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/kicktable.cpp b/src/kicktable.cpp index aace606..be9a04b 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -135,6 +135,8 @@ bool Kicktable::operator==(const Kicktable& o) const { if (this->x_max != o.x_max) return false; if (this->y_min != o.y_min) return false; if (this->y_max != o.y_max) return false; + if (this->x_nrpts != o.x_nrpts) return false; + if (this->y_nrpts != o.y_nrpts) return false; if (this->x_kick != o.x_kick) return false; if (this->y_kick != o.y_kick) return false; return true; From 62edf716e792575458777027ead1772e2d6c751b Mon Sep 17 00:00:00 2001 From: Fernando Date: Wed, 21 Aug 2024 07:48:04 -0300 Subject: [PATCH 07/28] KICKTAB.ENH: change comparison in add_kicktable with the == operator instead of the filename. --- src/kicktable.cpp | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/kicktable.cpp b/src/kicktable.cpp index be9a04b..a508e3b 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -109,18 +109,21 @@ int add_kicktable(const std::string& filename) { for(unsigned int i=0; i Date: Wed, 21 Aug 2024 07:49:35 -0300 Subject: [PATCH 08/28] KICKTAB.ENH: Add new method add_kicktable, to add kicktables based on its attributes rather than a filename. --- include/trackcpp/kicktable.h | 11 +++++++++ src/kicktable.cpp | 45 ++++++++++++++++++++++++++++++++---- 2 files changed, 52 insertions(+), 4 deletions(-) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index 9457287..cbdd3fe 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -55,6 +55,17 @@ class Kicktable { extern std::vector kicktable_list; +int add_kicktable( + const double x_min, + const double x_max, + const unsigned int x_nrpts, + const std::vector& x_kick, + const double y_min, + const double y_max, + const unsigned int y_nrpts, + const std::vector& y_kick, + const double length=1 +); int add_kicktable(const std::string& filename); void clear_kicktables(std::vector& kicktable_list); diff --git a/src/kicktable.cpp b/src/kicktable.cpp index a508e3b..8184321 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -103,12 +103,49 @@ Status::type Kicktable::load_from_file(const std::string& filename_) { } -int add_kicktable(const std::string& filename) { - + +int add_kicktable( + const double x_min, + const double x_max, + const unsigned int x_nrpts, + const std::vector& x_kick, + const double y_min, + const double y_max, + const unsigned int y_nrpts, + const std::vector& y_kick, + const double length +) +{ + if (x_kick.size() != x_nrpts * y_nrpts) + return -1; + if (y_kick.size() != x_nrpts * y_nrpts) + return -1; + if ((x_min >= x_max) | (y_min >= y_max)) + return -1; + + Kicktable new_kicktable(""); + new_kicktable.length = length; + + new_kicktable.x_min = x_min; + new_kicktable.x_max = x_max; + new_kicktable.x_nrpts = x_nrpts; + new_kicktable.x_kick = x_kick; + + new_kicktable.y_min = y_min; + new_kicktable.y_max = y_max; + new_kicktable.y_nrpts = y_nrpts; + new_kicktable.y_kick = y_kick; + // looks through vector of kicktables... - for(unsigned int i=0; i Date: Thu, 22 Aug 2024 07:46:51 -0300 Subject: [PATCH 09/28] DOC: few documentation changes. --- include/trackcpp/elements.h | 2 +- include/trackcpp/flat_file.h | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h index 56dee27..97f50ce 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -61,7 +61,7 @@ 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 polynom_a = default_polynom; diff --git a/include/trackcpp/flat_file.h b/include/trackcpp/flat_file.h index adb4e8a..3393ee8 100644 --- a/include/trackcpp/flat_file.h +++ b/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); From 7ab9ea9e95619bb5fe6f6e07f6793834594bf197 Mon Sep 17 00:00:00 2001 From: Fernando Date: Thu, 22 Aug 2024 10:08:12 -0300 Subject: [PATCH 10/28] KICKTAB.ENH: Improve Kicktable class - Now the variable kicktable_list is a static variable of the Kicktable class - the methods to add_kicktable and new Kicktable constructors were added to not necessarily depend on a filename. They are also static methods of the Kicktable class; - methods were added to check if a given Kicktable is valid and if a given index for the kicktable_list is valid --- include/trackcpp/kicktable.h | 200 ++++++++++++++++++------------- include/trackcpp/passmethods.hpp | 21 ++-- python_package/interface.cpp | 14 ++- python_package/interface.h | 8 +- python_package/trackcpp.i | 2 +- src/elements.cpp | 4 +- src/flat_file.cpp | 42 ++++++- src/kicktable.cpp | 97 ++++++++------- src/tests.cpp | 8 +- 9 files changed, 252 insertions(+), 144 deletions(-) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index cbdd3fe..e00d93e 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -34,103 +34,137 @@ class Kicktable { double x_min, x_max; double y_min, y_max; std::vector x_kick, y_kick; - + static std::vector kicktable_list; + + Kicktable( + const double x_min, + const double x_max, + const unsigned int x_nrpts, + const std::vector& x_kick, + const double y_min, + const double y_max, + const unsigned int y_nrpts, + const std::vector& 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 unsigned int get_ix(const T& x) const { return (int) ((x - x_min) / ((x_max - x_min) / (x_nrpts - 1))); } - template 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_list; - - -int add_kicktable( - const double x_min, - const double x_max, - const unsigned int x_nrpts, - const std::vector& x_kick, - const double y_min, - const double y_max, - const unsigned int y_nrpts, - const std::vector& y_kick, - const double length=1 -); -int add_kicktable(const std::string& filename); -void clear_kicktables(std::vector& kicktable_list); - -template -Status::type kicktable_getkicks_bilinear(const int& kicktable_idx, const T& rx, const T& ry, T& hkick, T& vkick) { + bool is_valid_kicktable() const; + 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); + } - const Kicktable* kicktable = &kicktable_list[kicktable_idx]; + template + unsigned int get_ix(const T& x) const + { + return (int) ((x - x_min) / ((x_max - x_min) / (x_nrpts - 1))); + } + template + unsigned int get_iy(const T& y) const + { + return (int) ((y - y_min) / ((y_max - y_min) / (y_nrpts - 1))); + } - // 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 + Status::type getkicks_bilinear( + const T& rx, const T& ry, T& hkick, T& vkick + ) const + { + // checks x limits + if ((rx < x_min) or (rx > x_max)) { + //std::cout << "rx: " << double(rx) << std::endl; + hkick = nan(""); + return Status::kicktable_out_of_range; + } + // checks y limits + if ((ry < y_min) or (ry > y_max)) { + //std::cout << "ry: " << double(ry) << std::endl; + vkick = nan(""); + return Status::kicktable_out_of_range; + } + + // gets indices + const unsigned int ix = get_ix(rx); + const unsigned int iy = get_iy(ry); + + /* coordinates */ + const double x1 = get_x(ix); + const double x2 = get_x(ix+1); + const double y1 = get_y(iy); + const double y2 = get_y(iy+1); + + // Bilinear interpolation + // https://en.wikipedia.org/wiki/Bilinear_interpolation + { /* hkick */ + const double fq11 = x_kick[get_idx(ix+0, iy+0)]; + const double fq12 = x_kick[get_idx(ix+0, iy+1)]; + const double fq21 = x_kick[get_idx(ix+1, iy+0)]; + const double fq22 = x_kick[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; + } + { /* vkick */ + const double fq11 = y_kick[get_idx(ix+0, iy+0)]; + const double fq12 = y_kick[get_idx(ix+0, iy+1)]; + const double fq21 = y_kick[get_idx(ix+1, iy+0)]; + const double fq22 = y_kick[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; + } + + return Status::success; } - // 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; + template + Status::type getkicks(const T& rx, const T& ry, T& hkick, T& vkick) const + { + return getkicks_bilinear(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 double x_min, + const double x_max, + const unsigned int x_nrpts, + const std::vector& x_kick, + const double y_min, + const double y_max, + const unsigned int y_nrpts, + const std::vector& y_kick, + const double length=1 + ); + static int add_kicktable(const std::string& filename); + static int add_kicktable(const Kicktable &new_kicktable); + 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."); + return kicktable_list[kicktable_idx]; } - { /* 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(); - return Status::success; -} + bool operator==(const Kicktable& o) const; + bool operator!=(const Kicktable& o) const { return !(*this == o); } -template -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 diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 2b38665..c4a1c81 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -97,14 +97,14 @@ T b2_perp(const T& bx, const T& by, const T& px, const T& py, const T& curv=1) { template Status::type kicktablethinkick( Pos& pos, - const int& kicktable_idx, + const Kicktable& kicktable, 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); + Status::type status = kicktable.getkicks(pos.rx, pos.ry, hkick, vkick); // According to Ellaune's theory of kick maps: // https://accelconf.web.cern.ch/e92/PDF/EPAC1992_0661.PDF // @@ -507,14 +507,19 @@ Status::type pm_thinsext_pass(Pos &pos, const Element &elem, template -Status::type pm_kickmap_pass(Pos &pos, const Element &elem, - const Accelerator& accelerator) { - - if (elem.kicktable_idx < 0 or elem.kicktable_idx >= kicktable_list.size()) return Status::kicktable_not_defined; +Status::type pm_kickmap_pass( + Pos &pos, + const Element &elem, + const Accelerator& accelerator +) +{ + if (not Kicktable::is_valid_kicktable_index(elem.kicktable_idx)) + return Status::kicktable_not_defined; Status::type status = Status::success; + const Kicktable& kicktable = Kicktable::get_kicktable(elem.kicktable_idx); - double sl = elem.length / float(elem.nr_steps); + double sl = elem.length / float(elem.nr_steps); double brho2 = get_magnetic_rigidity(accelerator.energy); brho2 *= brho2; @@ -526,7 +531,7 @@ Status::type pm_kickmap_pass(Pos &pos, const Element &elem, drift(pos, sl / 2); Status::type status = kicktablethinkick( pos, - elem.kicktable_idx, + kicktable, brho2, elem.nr_steps, elem.rescale_kicks diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 57d435a..caa18db 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -249,8 +249,18 @@ Status::type write_flat_file_wrapper(String& fname, const Accelerator& accelerat return write_flat_file(fname.data, accelerator, file_flag); } -Status::type kicktable_getkicks_wrapper(const int& kicktable_idx, const double& rx, const double& ry, double& hkick__, double& vkick__) { - return kicktable_getkicks(kicktable_idx, rx, ry, hkick__, vkick__); +Status::type kicktable_getkicks_wrapper( + const int& kicktable_idx, + const double& rx, + const double& ry, + double& hkick__, + double& vkick__ +) +{ + if (not Kicktable::is_valid_kicktable_index(kicktable_idx)) + return Status::kicktable_not_defined; + const Kicktable &kicktable = Kicktable::get_kicktable(kicktable_idx); + return kicktable.getkicks(rx, ry, hkick__, vkick__); } Status::type track_findm66_wrapper( diff --git a/python_package/interface.h b/python_package/interface.h index 517f26a..32f3a2b 100644 --- a/python_package/interface.h +++ b/python_package/interface.h @@ -99,7 +99,13 @@ Element rbend_wrapper(const std::string& fam_name_, const double& length_, const double& K_, const double& S_); Status::type write_flat_file_wrapper(String& fname, const Accelerator& accelerator, bool file_flag = true); Status::type read_flat_file_wrapper(String& fname, Accelerator& accelerator, bool file_flag = true); -Status::type kicktable_getkicks_wrapper(const int& kicktable_idx, const double& rx, const double& ry, double& hkick__, double& vkick__); +Status::type kicktable_getkicks_wrapper( + const int& kicktable_idx, + const double& rx, + const double& ry, + double& hkick__, + double& vkick__ +); Status::type track_findm66_wrapper( Accelerator& accelerator, diff --git a/python_package/trackcpp.i b/python_package/trackcpp.i index 766a189..3c2889f 100644 --- a/python_package/trackcpp.i +++ b/python_package/trackcpp.i @@ -110,7 +110,7 @@ double get_double_max() { %apply (double* INPLACE_ARRAY2, int DIM1, int DIM2 ) { (double *im_out, int n1_im_out, int n2_im_out)} -// For kicktable_getkicks +// For kicktable_getkicks_wrapper %apply double *OUTPUT { double &hkick__, double &vkick__ }; diff --git a/src/elements.cpp b/src/elements.cpp index 7a75f35..0f93d9b 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -148,9 +148,9 @@ Element Element::rfcavity (const std::string& fam_name_, const double& length_, Element Element::kickmap (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_, const double &rescale_length_, const double &rescale_kicks_) { // add new kicktable to global list, if necessary. - int idx = add_kicktable(kicktable_fname_); + int idx = Kicktable::add_kicktable(kicktable_fname_); - const Kicktable& kicktable = kicktable_list[idx]; + const Kicktable& kicktable = Kicktable::get_kicktable(idx); Element e = Element(fam_name_, rescale_length_ * kicktable.length); initialize_kickmap(e, idx, nr_steps_, rescale_kicks_); return e; diff --git a/src/flat_file.cpp b/src/flat_file.cpp index d241409..c0e3853 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -103,7 +103,7 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) fp << std::setw(pw) << "pass_method" << pm_dict[e.pass_method] << '\n'; if (pm_dict[e.pass_method].compare("kicktable_pass") == 0) { unsigned int idx = e.kicktable_idx; - const Kicktable& ktable = kicktable_list[idx]; + const Kicktable& ktable = Kicktable::get_kicktable(idx); fp << std::setw(pw) << "kicktable_fname" << ktable.filename << '\n'; } if (e.nr_steps != 1) { @@ -175,6 +175,15 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) bool found_hmin = false; bool found_vmin = false; unsigned int line_count = 0; + // once the counter below reaches 9, all info is available; + // x_max, x_min, x_nrpts, x_kick, y_max, y_min, y_nrpts, y_kick, length + // the filename property is not important: + unsigned int kicktable_ready = 0; + double kicktable_x_max, kicktable_x_min; + double kicktable_y_max, kicktable_y_min, kicktable_length; + unsigned int x_nrpts, y_nrpts; + std::vector kicktable_x_kick, kicktable_y_kick; + std::string kicktable_filename; while (std::getline(fp, line)) { //std::cout << line << std::endl; line_count++; @@ -266,7 +275,7 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) } if (cmd.compare("kicktable_fname") == 0) { std::string fname; ss >> fname; - int idx = add_kicktable(fname); + int idx = Kicktable::add_kicktable(fname); if (idx < 0) { return Status::file_not_found; } else { @@ -274,6 +283,31 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) continue; } } + // if (cmd.compare("kicktable_x_kick") == 0) { + // unsigned int size = 0; + // while (not ss.eof()) { + // double m; + // ss >> m; + // if (ss.eof()) break; + // kicktable_x_kick.push_back(m); + // } + // ++kicktable_ready; + // if (kicktable_ready >= 9) { + // int idx = Kicktable::add_kicktable( + // kicktable_x_min, + // kicktable_x_max, + // kicktable_x_nrpts, + // kicktable_x_kick, + // kicktable_y_min, + // kicktable_y_max, + // kicktable_y_nrpts, + // kicktable_y_kick, + // kicktable_length + // ); + // if (idx >= 0) + // e.kicktable_idx = idx; + // } + // } if (cmd.compare("polynom_a") == 0) { std::vector order; std::vector multipole; @@ -438,10 +472,10 @@ static Status::type read_flat_file_tracy(const std::string& filename, Accelerato e.pass_method = PassMethod::pm_kickmap_pass; double tmpdbl; std::string filename; fp >> tmpdbl >> tmpdbl >> filename; - int idx = add_kicktable(filename); + int idx = Kicktable::add_kicktable(filename); if (idx >= 0) { e.kicktable_idx = idx; - e.length = kicktable_list[e.kicktable_idx].length; + e.length = Kicktable::get_kicktable(e.kicktable_idx).length; //std::cout << accelerator.lattice.size() << " " << e.fam_name << ": " << e.kicktable << " " << e.kicktable->x_nrpts << std::endl; } else return Status::file_not_found; diff --git a/src/kicktable.cpp b/src/kicktable.cpp index 8184321..c510014 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -20,18 +20,52 @@ #include #include -std::vector kicktable_list; + +std::vector Kicktable::kicktable_list; Kicktable::Kicktable(const std::string& filename_) : filename(""), x_nrpts(0), y_nrpts(0), x_min(nan("")), x_max(nan("")), - y_min(nan("")), y_max(nan("")) { - - if (filename_ != "") { + y_min(nan("")), y_max(nan("")) +{ + if (filename_ != "") this->load_from_file(filename_); - } +} + +Kicktable::Kicktable( + const double x_min, + const double x_max, + const unsigned int x_nrpts, + const std::vector& x_kick, + const double y_min, + const double y_max, + const unsigned int y_nrpts, + const std::vector& y_kick, + const double length +) : + filename(""), + x_nrpts(x_nrpts), y_nrpts(y_nrpts), + x_min(x_min), x_max(x_max), + y_min(y_min), y_max(y_max), + x_kick(x_kick), y_kick(y_kick), + length(length) +{} + +bool Kicktable::is_valid_kicktable() const +{ + if (x_nrpts <= 1) + return false; + if (y_nrpts <= 1) + return false; + if (x_kick.size() != x_nrpts * y_nrpts) + return false; + if (y_kick.size() != x_nrpts * y_nrpts) + return false; + if ((x_min >= x_max) | (y_min >= y_max)) + return false; + return true; } Status::type Kicktable::load_from_file(const std::string& filename_) { @@ -104,7 +138,7 @@ Status::type Kicktable::load_from_file(const std::string& filename_) { } -int add_kicktable( +int Kicktable::add_kicktable( const double x_min, const double x_max, const unsigned int x_nrpts, @@ -116,42 +150,26 @@ int add_kicktable( const double length ) { - if (x_kick.size() != x_nrpts * y_nrpts) - return -1; - if (y_kick.size() != x_nrpts * y_nrpts) - return -1; - if ((x_min >= x_max) | (y_min >= y_max)) - return -1; - - Kicktable new_kicktable(""); - new_kicktable.length = length; - - new_kicktable.x_min = x_min; - new_kicktable.x_max = x_max; - new_kicktable.x_nrpts = x_nrpts; - new_kicktable.x_kick = x_kick; + Kicktable new_kicktable = Kicktable( + x_min, x_max, x_nrpts, x_kick, + y_min, y_max, y_nrpts, y_kick, + length + ); + return Kicktable::add_kicktable(new_kicktable); +} - new_kicktable.y_min = y_min; - new_kicktable.y_max = y_max; - new_kicktable.y_nrpts = y_nrpts; - new_kicktable.y_kick = y_kick; - // looks through vector of kicktables... - for(unsigned int i=0; i& kicktable_list) { +void Kicktable::clear_kicktables() { kicktable_list.clear(); } @@ -181,5 +198,3 @@ bool Kicktable::operator==(const Kicktable& o) const { if (this->y_kick != o.y_kick) return false; return true; } - -template Status::type kicktable_getkicks(const int& kicktable_idx, const double& rx, const double& ry, double& hkick, double& vkick); diff --git a/src/tests.cpp b/src/tests.cpp index 8ed7ead..68c981a 100644 --- a/src/tests.cpp +++ b/src/tests.cpp @@ -323,8 +323,12 @@ int test_kicktable(Accelerator& accelerator) { Kicktable t; const Kicktable *ptrKicktable = nullptr; - add_kicktable("/home/fac_files/code/python/trackcpp/pytrack/id_kicktable.txt"); - add_kicktable("/home/fac_files/code/python/trackcpp/pytrack/id_kicktable2.txt"); + Kicktable::add_kicktable( + "/home/fac_files/code/python/trackcpp/pytrack/id_kicktable.txt" + ); + Kicktable::add_kicktable( + "/home/fac_files/code/python/trackcpp/pytrack/id_kicktable2.txt" + ); return 0; } From b60694b04e511fbe4276c46d34cf059a3af1c7fb Mon Sep 17 00:00:00 2001 From: Fernando Date: Fri, 23 Aug 2024 16:44:35 -0300 Subject: [PATCH 11/28] KICKTAB.ENH: Extend Kicktable to non-regular grids. --- include/trackcpp/kicktable.h | 100 ++++++++++++++---------------- src/kicktable.cpp | 114 ++++++++++++++--------------------- 2 files changed, 91 insertions(+), 123 deletions(-) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index e00d93e..2432575 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -18,8 +18,10 @@ #define _KICKTABLE_H #include "auxiliary.h" +// #include "../alglib/interpolation.h" #include #include +#include class Kicktable { @@ -30,20 +32,14 @@ class Kicktable { std::string filename; double length; - unsigned int x_nrpts, y_nrpts; - double x_min, x_max; - double y_min, y_max; + std::vector x_pos, y_pos; std::vector x_kick, y_kick; static std::vector kicktable_list; Kicktable( - const double x_min, - const double x_max, - const unsigned int x_nrpts, + const std::vector& x_pos, const std::vector& x_kick, - const double y_min, - const double y_max, - const unsigned int y_nrpts, + const std::vector& y_pos, const std::vector& y_kick, const double length = 1 ); @@ -55,30 +51,20 @@ class Kicktable { bool is_valid_kicktable() const; unsigned int get_idx(unsigned int ix, unsigned int iy) const { - return iy*x_nrpts+ix; + return iy * x_pos.size() + ix; } - double get_x(unsigned int ix) const + double get_x(unsigned int ix) const {return x_pos[ix];} + double get_y(unsigned int iy) const {return y_pos[iy];} + template unsigned int get_ix(const T& x) const { - return x_min + ix * (x_max - x_min) / (x_nrpts - 1.0); + return std::lower_bound(x_pos.begin(), x_pos.end(), x) - x_pos.begin(); } - double get_y(unsigned int iy) const + template unsigned int get_iy(const T& y) const { - return y_min + iy * (y_max - y_min) / (y_nrpts - 1.0); + return std::lower_bound(y_pos.begin(), y_pos.end(), y) - y_pos.begin(); } - template - unsigned int get_ix(const T& x) const - { - return (int) ((x - x_min) / ((x_max - x_min) / (x_nrpts - 1))); - } - template - unsigned int get_iy(const T& y) const - { - return (int) ((y - y_min) / ((y_max - y_min) / (y_nrpts - 1))); - } - - template - Status::type getkicks_bilinear( + template Status::type getkicks_bilinear( const T& rx, const T& ry, T& hkick, T& vkick ) const { @@ -96,55 +82,59 @@ class Kicktable { } // gets indices - const unsigned int ix = get_ix(rx); - const unsigned int iy = get_iy(ry); + unsigned int ix = get_ix(rx); + unsigned int ixp1, iyp1; + if (ix >= x_pos.size()-1){ixp1 = ix; --ix;} + else ixp1 = ix + 1; + + unsigned int iy = get_iy(ry); + if (iy >= y_pos.size()-1){iyp1 = iy; --iy;} + else iyp1 = iy + 1; /* coordinates */ const double x1 = get_x(ix); - const double x2 = get_x(ix+1); + const double x2 = get_x(ixp1); const double y1 = get_y(iy); - const double y2 = get_y(iy+1); + const double y2 = get_y(iyp1); + double denom = (x2 - x1) * (y2 - y1); // Bilinear interpolation // https://en.wikipedia.org/wiki/Bilinear_interpolation { /* hkick */ - const double fq11 = x_kick[get_idx(ix+0, iy+0)]; - const double fq12 = x_kick[get_idx(ix+0, iy+1)]; - const double fq21 = x_kick[get_idx(ix+1, iy+0)]; - const double fq22 = x_kick[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; + const double fq11 = x_kick[get_idx(ix, iy)]; + const double fq12 = x_kick[get_idx(ix, iyp1)]; + const double fq21 = x_kick[get_idx(ixp1, iy)]; + const double fq22 = x_kick[get_idx(ixp1, iyp1)]; + hkick = + f11 * (x2 - x) * (y2 - y) + f21 * (x - x1) * (y2 - y) + + f12 * (x2 - x) * (y - y1) + f22 * (x - x1) * (y - y1); + hkick /= denom; } { /* vkick */ - const double fq11 = y_kick[get_idx(ix+0, iy+0)]; - const double fq12 = y_kick[get_idx(ix+0, iy+1)]; - const double fq21 = y_kick[get_idx(ix+1, iy+0)]; - const double fq22 = y_kick[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; + const double fq11 = y_kick[get_idx(ix, iy)]; + const double fq12 = y_kick[get_idx(ix, iyp1)]; + const double fq21 = y_kick[get_idx(ixp1, iy)]; + const double fq22 = y_kick[get_idx(ixp1, iyp1)]; + vkick = + f11 * (x2 - x) * (y2 - y) + f21 * (x - x1) * (y2 - y) + + f12 * (x2 - x) * (y - y1) + f22 * (x - x1) * (y - y1); + vkick /= denom; } return Status::success; } - template - Status::type getkicks(const T& rx, const T& ry, T& hkick, T& vkick) const + template Status::type getkicks( + const T& rx, const T& ry, T& hkick, T& vkick + ) const { return getkicks_bilinear(rx, ry, hkick, vkick); } static int add_kicktable( - const double x_min, - const double x_max, - const unsigned int x_nrpts, + const std::vector& x_pos, const std::vector& x_kick, - const double y_min, - const double y_max, - const unsigned int y_nrpts, + const std::vector& y_pos, const std::vector& y_kick, const double length=1 ); diff --git a/src/kicktable.cpp b/src/kicktable.cpp index c510014..11bc3d2 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -24,10 +24,7 @@ std::vector Kicktable::kicktable_list; Kicktable::Kicktable(const std::string& filename_) : - filename(""), - x_nrpts(0), y_nrpts(0), - x_min(nan("")), x_max(nan("")), - y_min(nan("")), y_max(nan("")) + filename("") { if (filename_ != "") this->load_from_file(filename_); @@ -35,35 +32,32 @@ Kicktable::Kicktable(const std::string& filename_) : Kicktable::Kicktable( - const double x_min, - const double x_max, - const unsigned int x_nrpts, + const std::vector& x_pos, const std::vector& x_kick, - const double y_min, - const double y_max, - const unsigned int y_nrpts, + const std::vector& y_pos, const std::vector& y_kick, const double length ) : filename(""), - x_nrpts(x_nrpts), y_nrpts(y_nrpts), - x_min(x_min), x_max(x_max), - y_min(y_min), y_max(y_max), + x_pos(x_pos), y_pos(y_pos), x_kick(x_kick), y_kick(y_kick), length(length) {} bool Kicktable::is_valid_kicktable() const { - if (x_nrpts <= 1) + + if (x_pos.size() <= 1) return false; - if (y_nrpts <= 1) + if (y_pos.size() <= 1) return false; - if (x_kick.size() != x_nrpts * y_nrpts) + if (x_kick.size() != x_pos.size() * y_pos.size()) return false; - if (y_kick.size() != x_nrpts * y_nrpts) + if (y_kick.size() != x_kick.size()) return false; - if ((x_min >= x_max) | (y_min >= y_max)) + if (not std::is_sorted(x_pos.cbegin(), x_pos.cend())) + return false; + if (not std::is_sorted(y_pos.cbegin(), y_pos.cend())) return false; return true; } @@ -75,61 +69,57 @@ Status::type Kicktable::load_from_file(const std::string& filename_) { std::cout << "Could not find kicktable file " << filename_ << "!" << std::endl; return Status::file_not_found; } - this->filename = filename_; + filename = filename_; std::string str; + unsigned int x_nrpts, y_nrpts; // HEADER - getline(fp, str); // name of kicktable line - getline(fp, str); // author line - getline(fp, str); // label 'ID length[m]' - fp >> this->length; // length of element - getline(fp, str); // advances to new line - getline(fp, str); // label 'number of horizontal points' - fp >> this->x_nrpts; // number of horizontal points - getline(fp, str); // advances to new line - getline(fp, str); // label 'number of vertical points' - fp >> this->y_nrpts; // number of vertical points - getline(fp, str); // advances to new line - - this->x_kick.resize(x_nrpts * y_nrpts, 0); - this->y_kick.resize(x_nrpts * y_nrpts, 0); - - std::vector yvec(y_nrpts); // used to invert tables ordering, if necessary + getline(fp, str); // name of kicktable line + getline(fp, str); // author line + getline(fp, str); // label 'ID length[m]' + fp >> length; // length of element + getline(fp, str); // advances to new line + getline(fp, str); // label 'number of horizontal points' + fp >> x_nrpts; // number of horizontal points + getline(fp, str); // advances to new line + getline(fp, str); // label 'number of vertical points' + fp >> y_nrpts; // number of vertical points + getline(fp, str); // advances to new line + + x_kick.resize(x_nrpts * y_nrpts, 0); + y_kick.resize(x_nrpts * y_nrpts, 0); + x_pos.resize(x_nrpts, 0); + // used to invert tables ordering, if necessary + y_pos.resize(y_nrpts, 0); // HORIZONTAL KICK TABLE getline(fp, str); // label 'Horizontal KickTable in T^2.m^2' getline(fp, str); // label 'START' - for(unsigned int i=0; ix_nrpts; ++i) { - double posx; fp >> posx; - if (std::isnan(x_min) or posx < x_min) x_min = posx; - if (std::isnan(x_max) or posx > x_max) x_max = posx; - } + for(unsigned int i=0; i> x_pos[i];} for(int j=y_nrpts-1; j>=0; --j) { - double posy; fp >> posy; yvec.push_back(posy); - if (std::isnan(y_min) or posy < y_min) y_min = posy; - if (std::isnan(y_max) or posy > y_max) y_max = posy; + fp >> y_pos[j]; for(unsigned int i=0; i> x_kick[this->get_idx(i,j)]; + fp >> x_kick[get_idx(i, j)]; } getline(fp, str); // advances to new line // VERTICAL KICK TABLE getline(fp, str); // label 'Vertical KickTable in T^2.m^2' getline(fp, str); // label 'START' - for(unsigned int i=0; ix_nrpts; ++i) { double posx; fp >> posx; } + getline(fp, str); // horizontal position. Already set from x_kick for(int j=y_nrpts-1; j>=0; --j) { double posy; fp >> posy; for(unsigned int i=0; i> y_kick[this->get_idx(i,j)]; + fp >> y_kick[get_idx(i, j)]; } // invert tables, if necessary - if (yvec.size() > 1 and yvec[1] > yvec[0]) { - for(unsigned int i=0; ix_nrpts; ++i) { - for(unsigned int j=0; jy_nrpts/2; ++j) { - std::swap(x_kick[this->get_idx(i,j)], x_kick[this->get_idx(i,this->y_nrpts-j-1)]); - std::swap(y_kick[this->get_idx(i,j)], y_kick[this->get_idx(i,this->y_nrpts-j-1)]); + if (y_pos.size() > 1 and y_pos[1] > y_pos[0]) { + for(unsigned int i=0; i& x_pos, const std::vector& x_kick, - const double y_min, - const double y_max, - const unsigned int y_nrpts, + const std::vector& y_pos, const std::vector& y_kick, const double length ) { - Kicktable new_kicktable = Kicktable( - x_min, x_max, x_nrpts, x_kick, - y_min, y_max, y_nrpts, y_kick, - length - ); + Kicktable new_kicktable = Kicktable(x_pos, x_kick, y_pos, y_kick, length); return Kicktable::add_kicktable(new_kicktable); } @@ -188,12 +170,8 @@ void Kicktable::clear_kicktables() { bool Kicktable::operator==(const Kicktable& o) const { if (this == &o) return true; if (this->length != o.length) return false; - if (this->x_min != o.x_min) return false; - if (this->x_max != o.x_max) return false; - if (this->y_min != o.y_min) return false; - if (this->y_max != o.y_max) return false; - if (this->x_nrpts != o.x_nrpts) return false; - if (this->y_nrpts != o.y_nrpts) return false; + if (this->x_pos != o.x_pos) return false; + if (this->y_pos != o.y_pos) return false; if (this->x_kick != o.x_kick) return false; if (this->y_kick != o.y_kick) return false; return true; From cefadfaa8818ddabcf17f2d28dd43fb6a6ea628c Mon Sep 17 00:00:00 2001 From: Fernando Date: Fri, 23 Aug 2024 16:45:02 -0300 Subject: [PATCH 12/28] PASSM.DOC: Fix typo in documentation. thanks to @xresende. --- include/trackcpp/passmethods.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index c4a1c81..0763b54 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -113,7 +113,7 @@ Status::type kicktablethinkick( // 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 + // behavior of the kicks as function of energy, corroborating with no or // 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 @@ -168,7 +168,7 @@ void kickpolythinkick( // 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 + // behavior of the kicks as function of energy, corroborating with no or // 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 From e8f1f9fc59763719a5473b5a180e830f2e19ae50 Mon Sep 17 00:00:00 2001 From: Fernando Date: Mon, 26 Aug 2024 10:35:25 -0300 Subject: [PATCH 13/28] KICKTAB.BUG: Fix bugs of last commit. add rhs comparison of Tpsa objects with double. --- include/trackcpp/kicktable.h | 37 ++++++++++++----------------------- include/trackcpp/tpsa.h | 38 ++++++++++++++++++++++++++++++++++-- src/kicktable.cpp | 1 - 3 files changed, 48 insertions(+), 28 deletions(-) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index 2432575..8a34b55 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -68,19 +68,6 @@ class Kicktable { const T& rx, const T& ry, T& hkick, T& vkick ) const { - // checks x limits - if ((rx < x_min) or (rx > x_max)) { - //std::cout << "rx: " << double(rx) << std::endl; - hkick = nan(""); - return Status::kicktable_out_of_range; - } - // checks y limits - if ((ry < y_min) or (ry > y_max)) { - //std::cout << "ry: " << double(ry) << std::endl; - vkick = nan(""); - return Status::kicktable_out_of_range; - } - // gets indices unsigned int ix = get_ix(rx); unsigned int ixp1, iyp1; @@ -101,23 +88,23 @@ class Kicktable { // Bilinear interpolation // https://en.wikipedia.org/wiki/Bilinear_interpolation { /* hkick */ - const double fq11 = x_kick[get_idx(ix, iy)]; - const double fq12 = x_kick[get_idx(ix, iyp1)]; - const double fq21 = x_kick[get_idx(ixp1, iy)]; - const double fq22 = x_kick[get_idx(ixp1, iyp1)]; + const double f11 = x_kick[get_idx(ix, iy)]; + const double f12 = x_kick[get_idx(ix, iyp1)]; + const double f21 = x_kick[get_idx(ixp1, iy)]; + const double f22 = x_kick[get_idx(ixp1, iyp1)]; hkick = - f11 * (x2 - x) * (y2 - y) + f21 * (x - x1) * (y2 - y) + - f12 * (x2 - x) * (y - y1) + f22 * (x - x1) * (y - y1); + f11 * (x2 - rx) * (y2 - ry) + f21 * (rx - x1) * (y2 - ry) + + f12 * (x2 - rx) * (ry - y1) + f22 * (rx - x1) * (ry - y1); hkick /= denom; } { /* vkick */ - const double fq11 = y_kick[get_idx(ix, iy)]; - const double fq12 = y_kick[get_idx(ix, iyp1)]; - const double fq21 = y_kick[get_idx(ixp1, iy)]; - const double fq22 = y_kick[get_idx(ixp1, iyp1)]; + const double f11 = y_kick[get_idx(ix, iy)]; + const double f12 = y_kick[get_idx(ix, iyp1)]; + const double f21 = y_kick[get_idx(ixp1, iy)]; + const double f22 = y_kick[get_idx(ixp1, iyp1)]; vkick = - f11 * (x2 - x) * (y2 - y) + f21 * (x - x1) * (y2 - y) + - f12 * (x2 - x) * (y - y1) + f22 * (x - x1) * (y - y1); + f11 * (x2 - rx) * (y2 - ry) + f21 * (rx - x1) * (y2 - ry) + + f12 * (x2 - rx) * (ry - y1) + f22 * (rx - x1) * (ry - y1); vkick /= denom; } diff --git a/include/trackcpp/tpsa.h b/include/trackcpp/tpsa.h index 580a2cf..391baa2 100644 --- a/include/trackcpp/tpsa.h +++ b/include/trackcpp/tpsa.h @@ -46,8 +46,6 @@ // of using FFT are interesting since the Beam Dynamics community does not seem to have // realized that the method could speed up calculation of taylor maps... - - #ifndef TPSA_H #define TPSA_H @@ -375,6 +373,42 @@ bool Tpsa::operator >= (const TYPE& o_) const { return false; } +template +inline bool operator == (const TYPE& lhs, const Tpsa& rhs) +{ + return rhs == lhs; +} + +template +inline bool operator != (const TYPE& lhs, const Tpsa& rhs) +{ + return rhs != lhs; +} + +template +inline bool operator < (const TYPE& lhs, const Tpsa& rhs) +{ + return rhs > lhs; +} + +template +inline bool operator <= (const TYPE& lhs, const Tpsa& rhs) +{ + return rhs >= lhs; +} + +template +inline bool operator > (const TYPE& lhs, const Tpsa& rhs) +{ + return rhs < lhs; +} + +template +inline bool operator >= (const TYPE& lhs, const Tpsa& rhs) +{ + return rhs <= lhs; +} + template bool Tpsa::operator == (const Tpsa& o_) const { return (*this - o_) == (TYPE) 0; diff --git a/src/kicktable.cpp b/src/kicktable.cpp index 11bc3d2..ae2a0a7 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -127,7 +127,6 @@ Status::type Kicktable::load_from_file(const std::string& filename_) { } - int Kicktable::add_kicktable( const std::vector& x_pos, const std::vector& x_kick, From bae18911b46051f81e4c5e8d8d0f13decd9e9593 Mon Sep 17 00:00:00 2001 From: Fernando Date: Mon, 26 Aug 2024 10:36:33 -0300 Subject: [PATCH 14/28] TPSA.MNT: Remove duplicate include of math module. --- include/trackcpp/tpsa.h | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/include/trackcpp/tpsa.h b/include/trackcpp/tpsa.h index 391baa2..bf141c7 100644 --- a/include/trackcpp/tpsa.h +++ b/include/trackcpp/tpsa.h @@ -52,6 +52,7 @@ #include #include #include +#include // Expression Templates: IMPLEMENTATION OF BINOMIALS COEFFICIENTS AND RELATED RELEVANT EXPRESSIONS @@ -547,7 +548,6 @@ Tpsa abs(const Tpsa& a_) { if (a_ >= 0) return a_; else return -a_; } -#include template Tpsa sqrt(const Tpsa& a_) { Tpsa r; @@ -563,7 +563,6 @@ Tpsa sqrt(const Tpsa& a_) { return r; } -#include template Tpsa log(const Tpsa& a_) { Tpsa r; @@ -577,7 +576,6 @@ Tpsa log(const Tpsa& a_) { return r; } -#include template Tpsa cos(const Tpsa& a_) { Tpsa rc(1), rs; @@ -592,7 +590,6 @@ Tpsa cos(const Tpsa& a_) { return cos(a_.c[0]) * rc - sin(a_.c[0]) * rs; } -#include template Tpsa sin(const Tpsa& a_) { Tpsa rc(1), rs; @@ -607,13 +604,11 @@ Tpsa sin(const Tpsa& a_) { return sin(a_.c[0]) * rc + cos(a_.c[0]) * rs; } -#include template Tpsa tan(const Tpsa& a_) { return sin(a_)/cos(a_); } -#include template Tpsa cosh(const Tpsa& a_) { Tpsa rc(1), rs; @@ -628,7 +623,6 @@ Tpsa cosh(const Tpsa& a_) { return cosh(a_.c[0]) * rc + sinh(a_.c[0]) * rs; } -#include template Tpsa sinh(const Tpsa& a_) { Tpsa rc(1), rs; @@ -643,8 +637,6 @@ Tpsa sinh(const Tpsa& a_) { return sinh(a_.c[0]) * rc + cosh(a_.c[0]) * rs; } - -#include template Tpsa atan(const Tpsa& a_) { Tpsa r(atan(a_.c[0])); @@ -674,7 +666,6 @@ Tpsa atan(const Tpsa& a_) { return r; } -#include template Tpsa asin(const Tpsa& a_) { @@ -689,13 +680,11 @@ Tpsa asin(const Tpsa& a_) { } -#include template Tpsa acos(const Tpsa& a_) { return M_PI/2 - asin(a_); } -#include template Tpsa atanh(const Tpsa& a_) { Tpsa r(atan(a_.c[0])); From 1ae91ba998bded89408e7857c4365ab263d0f7ae Mon Sep 17 00:00:00 2001 From: Fernando Date: Mon, 26 Aug 2024 16:33:11 -0300 Subject: [PATCH 15/28] FLATFILE.ENH: Add kicktables to flat files. --- src/flat_file.cpp | 155 +++++++++++++++++++++++++++------------------- src/kicktable.cpp | 6 +- 2 files changed, 94 insertions(+), 67 deletions(-) diff --git a/src/flat_file.cpp b/src/flat_file.cpp index c0e3853..42cbc3a 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -21,6 +21,7 @@ #include #include #include +#include static const int hw = 18; // header field width static const int pw = 16; // parameter field width @@ -32,8 +33,8 @@ static bool has_t_vector(const double* t); static bool has_r_matrix(const double* r); static bool has_matrix66(const Matrix& r); static bool has_polynom(const std::vector& p); -static void write_6d_vector(std::ostream& fp, const std::string& label, const double* t); -static void write_6d_vector(std::ostream& fp, const std::string& label, const std::vector& t); +static void write_vector(std::ostream& fp, const std::string& label, const double* t, const std::size_t& siz); +static void write_vector(std::ostream& fp, const std::string& label, const std::vector& t); static void write_polynom(std::ostream& fp, const std::string& label, const std::vector& p); static void synchronize_polynomials(Element& e); static void synchronize_polynomials_kick(Element& e); @@ -103,8 +104,17 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) fp << std::setw(pw) << "pass_method" << pm_dict[e.pass_method] << '\n'; if (pm_dict[e.pass_method].compare("kicktable_pass") == 0) { unsigned int idx = e.kicktable_idx; - const Kicktable& ktable = Kicktable::get_kicktable(idx); - fp << std::setw(pw) << "kicktable_fname" << ktable.filename << '\n'; + const Kicktable& ktab = Kicktable::get_kicktable(idx); + if (ktab.filename.empty()) + { + write_vector(fp, "kicktable_x_pos", ktab.x_pos); + write_vector(fp, "kicktable_y_pos", ktab.y_pos); + write_vector(fp, "kicktable_x_kick", ktab.x_kick); + write_vector(fp, "kicktable_y_kick", ktab.y_kick); + fp << std::setw(pw) << "kicktable_length" << ktab.length << '\n'; + } + else + fp << std::setw(pw) << "kicktable_fname" << ktab.filename << '\n'; } if (e.nr_steps != 1) { fp.unsetf(std::ios_base::showpos); @@ -132,31 +142,31 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) if (e.angle_in != 0) { fp << std::setw(pw) << "angle_in" << e.angle_in << '\n'; } if (e.angle_out != 0) { fp << std::setw(pw) << "angle_out" << e.angle_out << '\n'; } if (e.rescale_kicks != 1.0) { fp << std::setw(pw) << "rescale_kicks" << e.rescale_kicks << '\n'; } - if (has_t_vector(e.t_in)) write_6d_vector(fp, "t_in", e.t_in); - if (has_t_vector(e.t_out)) write_6d_vector(fp, "t_out", e.t_out); + if (has_t_vector(e.t_in)) write_vector(fp, "t_in", e.t_in, 6); + if (has_t_vector(e.t_out)) write_vector(fp, "t_out", e.t_out, 6); if (has_r_matrix(e.r_in)) { - write_6d_vector(fp, "rx|r_in", &e.r_in[6*0]); - write_6d_vector(fp, "px|r_in", &e.r_in[6*1]); - write_6d_vector(fp, "ry|r_in", &e.r_in[6*2]); - write_6d_vector(fp, "py|r_in", &e.r_in[6*3]); - write_6d_vector(fp, "de|r_in", &e.r_in[6*4]); - write_6d_vector(fp, "dl|r_in", &e.r_in[6*5]); + write_vector(fp, "rx|r_in", &e.r_in[6*0], 6); + write_vector(fp, "px|r_in", &e.r_in[6*1], 6); + write_vector(fp, "ry|r_in", &e.r_in[6*2], 6); + write_vector(fp, "py|r_in", &e.r_in[6*3], 6); + write_vector(fp, "de|r_in", &e.r_in[6*4], 6); + write_vector(fp, "dl|r_in", &e.r_in[6*5], 6); } if (has_r_matrix(e.r_out)) { - write_6d_vector(fp, "rx|r_out", &e.r_out[6*0]); - write_6d_vector(fp, "px|r_out", &e.r_out[6*1]); - write_6d_vector(fp, "ry|r_out", &e.r_out[6*2]); - write_6d_vector(fp, "py|r_out", &e.r_out[6*3]); - write_6d_vector(fp, "de|r_out", &e.r_out[6*4]); - write_6d_vector(fp, "dl|r_out", &e.r_out[6*5]); + write_vector(fp, "rx|r_out", &e.r_out[6*0], 6); + write_vector(fp, "px|r_out", &e.r_out[6*1], 6); + write_vector(fp, "ry|r_out", &e.r_out[6*2], 6); + write_vector(fp, "py|r_out", &e.r_out[6*3], 6); + write_vector(fp, "de|r_out", &e.r_out[6*4], 6); + write_vector(fp, "dl|r_out", &e.r_out[6*5], 6); } if (has_matrix66(e.matrix66)) { - write_6d_vector(fp, "rx|matrix66", e.matrix66[0]); - write_6d_vector(fp, "px|matrix66", e.matrix66[1]); - write_6d_vector(fp, "ry|matrix66", e.matrix66[2]); - write_6d_vector(fp, "py|matrix66", e.matrix66[3]); - write_6d_vector(fp, "de|matrix66", e.matrix66[4]); - write_6d_vector(fp, "dl|matrix66", e.matrix66[5]); + write_vector(fp, "rx|matrix66", e.matrix66[0]); + write_vector(fp, "px|matrix66", e.matrix66[1]); + write_vector(fp, "ry|matrix66", e.matrix66[2]); + write_vector(fp, "py|matrix66", e.matrix66[3]); + write_vector(fp, "de|matrix66", e.matrix66[4]); + write_vector(fp, "dl|matrix66", e.matrix66[5]); } fp << '\n'; @@ -175,15 +185,16 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) bool found_hmin = false; bool found_vmin = false; unsigned int line_count = 0; - // once the counter below reaches 9, all info is available; - // x_max, x_min, x_nrpts, x_kick, y_max, y_min, y_nrpts, y_kick, length + // once the counter below reaches 5, all info is available; + // x_pos, y_pos, x_kick, y_kick, length // the filename property is not important: unsigned int kicktable_ready = 0; - double kicktable_x_max, kicktable_x_min; - double kicktable_y_max, kicktable_y_min, kicktable_length; - unsigned int x_nrpts, y_nrpts; - std::vector kicktable_x_kick, kicktable_y_kick; - std::string kicktable_filename; + double kicktable_length; + std::map> kicktable_info; + kicktable_info["x_pos"] = {}; + kicktable_info["y_pos"] = {}; + kicktable_info["x_kick"] = {}; + kicktable_info["y_kick"] = {}; while (std::getline(fp, line)) { //std::cout << line << std::endl; line_count++; @@ -283,31 +294,39 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) continue; } } - // if (cmd.compare("kicktable_x_kick") == 0) { - // unsigned int size = 0; - // while (not ss.eof()) { - // double m; - // ss >> m; - // if (ss.eof()) break; - // kicktable_x_kick.push_back(m); - // } - // ++kicktable_ready; - // if (kicktable_ready >= 9) { - // int idx = Kicktable::add_kicktable( - // kicktable_x_min, - // kicktable_x_max, - // kicktable_x_nrpts, - // kicktable_x_kick, - // kicktable_y_min, - // kicktable_y_max, - // kicktable_y_nrpts, - // kicktable_y_kick, - // kicktable_length - // ); - // if (idx >= 0) - // e.kicktable_idx = idx; - // } - // } + if (cmd.compare(0, 10, "kicktable_") == 0) { + if ( + auto item = kicktable_info.find(cmd.substr(10)); + item != kicktable_info.end() + ) + while (not ss.eof()) { + double m; ss >> m; + if (ss.eof()) break; + item->second.push_back(m); + } + else + ss >> kicktable_length; + + ++kicktable_ready; + if (kicktable_ready >= 5) { + int idx = Kicktable::add_kicktable( + kicktable_info["x_pos"], + kicktable_info["y_pos"], + kicktable_info["x_kick"], + kicktable_info["y_kick"], + kicktable_length + ); + if (idx >= 0) + e.kicktable_idx = idx; + else + return Status::kicktable_not_defined; + kicktable_ready = 0; + for (auto& n : kicktable_info) + // n.first: key, n.second: value + n.second = {}; + } + continue; + } if (cmd.compare("polynom_a") == 0) { std::vector order; std::vector multipole; @@ -599,18 +618,26 @@ static bool has_polynom(const std::vector& p) { return false; } -static void write_6d_vector(std::ostream& fp, const std::string& label, const double* t) { +static void write_vector( + std::ostream& fp, + const std::string& label, + const double* t, + const std::size_t& siz +) +{ fp << std::setw(pw) << label; - for (int i=0; i<6; ++i) - fp << t[i] << " "; + for (auto i=0; i& t) { - fp << std::setw(pw) << label; - for (int i=0; i<6; ++i) - fp << t[i] << " "; - fp << '\n'; +static void write_vector( + std::ostream& fp, + const std::string& label, + const std::vector& t +) +{ + write_vector(fp, label, t.data(), t.size()); } static void write_polynom(std::ostream& fp, const std::string& label, const std::vector& p) { diff --git a/src/kicktable.cpp b/src/kicktable.cpp index ae2a0a7..b339360 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -33,8 +33,8 @@ Kicktable::Kicktable(const std::string& filename_) : Kicktable::Kicktable( const std::vector& x_pos, - const std::vector& x_kick, const std::vector& y_pos, + const std::vector& x_kick, const std::vector& y_kick, const double length ) : @@ -129,13 +129,13 @@ Status::type Kicktable::load_from_file(const std::string& filename_) { int Kicktable::add_kicktable( const std::vector& x_pos, - const std::vector& x_kick, const std::vector& y_pos, + const std::vector& x_kick, const std::vector& y_kick, const double length ) { - Kicktable new_kicktable = Kicktable(x_pos, x_kick, y_pos, y_kick, length); + Kicktable new_kicktable = Kicktable(x_pos, y_pos, x_kick, y_kick, length); return Kicktable::add_kicktable(new_kicktable); } From e8eed563c6f109f6eaebfb07c10a063cfa04709f Mon Sep 17 00:00:00 2001 From: Fernando Date: Tue, 27 Aug 2024 10:29:08 -0300 Subject: [PATCH 16/28] MNT: change standard to C++20. --- INSTALL.md | 4 ++-- Makefile | 6 +++--- python_package/Makefile | 2 +- tests/Makefile | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/INSTALL.md b/INSTALL.md index 7b14689..42189a8 100644 --- a/INSTALL.md +++ b/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++20 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++20 support (tested with GCC 4.8) - blas - gsl (GNU Scientific Library) - SWIG 2 (tested with SWIG 2.0.11) diff --git a/Makefile b/Makefile index d573057..45282f6 100644 --- a/Makefile +++ b/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++20 -fPIC +DBG_FLAG = -O0 -g3 -std=c++20 -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 diff --git a/python_package/Makefile b/python_package/Makefile index 06edce0..824f6fc 100644 --- a/python_package/Makefile +++ b/python_package/Makefile @@ -25,7 +25,7 @@ NUMPY_VERSION := $(shell $(PYTHON) -c "import numpy; print(numpy.__version__)") OBJECTS = $(WRAPPEROBJS) $(INTERFACEOBJS) -CPPFLAGS = -std=c++11 -fPIC $(OPT_FLAG) +CPPFLAGS = -std=c++20 -fPIC $(OPT_FLAG) LIBS = $(shell gsl-config --libs) LIBS += -L../build -ltrackcpp INC = $(shell gsl-config --cflags) diff --git a/tests/Makefile b/tests/Makefile index 3375bbb..c5f8e6f 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -31,8 +31,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++20 -fPIC +DBG_FLAG = -O0 -g3 -std=c++20 -fPIC ARFLAGS = rcs DFLAGS = -DVERSION=$(VERSION) From a6bfc3bbe269ecced15454c6ffa05295ca1aaf47 Mon Sep 17 00:00:00 2001 From: Fernando Date: Tue, 27 Aug 2024 12:18:38 -0300 Subject: [PATCH 17/28] KICKTAB.BUG: Fix get_idx. index next to correct was being returned. --- include/trackcpp/kicktable.h | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index 8a34b55..f9fa8c7 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -57,11 +57,15 @@ class Kicktable { double get_y(unsigned int iy) const {return y_pos[iy];} template unsigned int get_ix(const T& x) const { - return std::lower_bound(x_pos.begin(), x_pos.end(), x) - x_pos.begin(); + 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 unsigned int get_iy(const T& y) const { - return std::lower_bound(y_pos.begin(), y_pos.end(), y) - y_pos.begin(); + int idx; + idx = std::lower_bound(y_pos.begin(), y_pos.end(), y) - y_pos.begin() - 1; + return idx >= 0 ? (unsigned int) idx : 0u; } template Status::type getkicks_bilinear( @@ -70,12 +74,12 @@ class Kicktable { { // 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; --ix;} - else ixp1 = ix + 1; - unsigned int iy = get_iy(ry); - if (iy >= y_pos.size()-1){iyp1 = iy; --iy;} + if (ix >= x_pos.size()-1) ixp1 = ix--; + else ixp1 = ix + 1; + if (iy >= y_pos.size()-1) iyp1 = iy--; else iyp1 = iy + 1; /* coordinates */ From 4cb2853b50f9720c6a09b5ed4e1e0861013de267 Mon Sep 17 00:00:00 2001 From: Fernando Date: Tue, 27 Aug 2024 12:28:20 -0300 Subject: [PATCH 18/28] KICKTAB.ENH: Add method to save Kicktable. --- include/trackcpp/kicktable.h | 5 +++ src/kicktable.cpp | 78 ++++++++++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index f9fa8c7..27561fd 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -47,6 +47,11 @@ class Kicktable { Kicktable(const Kicktable &) = default; Status::type load_from_file(const std::string& filename_); + Status::type save_to_file( + std::string& filename_, + const std::string author_name = " ", + bool file_flag = true + ); bool is_valid_kicktable() const; unsigned int get_idx(unsigned int ix, unsigned int iy) const diff --git a/src/kicktable.cpp b/src/kicktable.cpp index b339360..49265a6 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -124,6 +124,84 @@ Status::type Kicktable::load_from_file(const std::string& filename_) { } } return Status::success; +} + + +static const int hw = 18; // header field width +static const int pw = 16; // parameter field width +static const int np = 17; // number precision +Status::type Kicktable::save_to_file( + std::string& filename_, + const std::string author_name, + const bool file_flag +) +{ + // done with the help of chatgpt: + std::unique_ptr fp; + if (file_flag) + { + fp = std::make_unique(filename.c_str()); + if (!fp->good()) + return Status::file_not_found; + } + else + fp = std::make_unique(); + + fp->setf( + std::ios_base::left | + std::ios_base::scientific | + std::ios_base::uppercase + ); + fp->precision(np); + + // HEADER + *fp << "# Author: " << author_name << '\n'; + *fp << "#" << '\n'; + *fp << "# Total Length of Longitudinal Interval [m]" << '\n'; + *fp << length << '\n'; + *fp << "# Number of Horizontal Points" << '\n'; + *fp << x_pos.size() << '\n'; + *fp << "# Number of Vertical Points" << '\n'; + *fp << y_pos.size() << '\n'; + + // x_kick: + *fp << "# Total Horizontal 2nd Order Kick [T2m2]" << '\n'; + *fp << "START" << '\n'; + *fp << std::setw(pw) << ' '; + for (auto xi: x_pos) + *fp << std::setw(pw) << xi << ' '; + *fp << '\n'; + for (auto j=y_pos.size()-1; j>=0; ++j) + { + *fp << std::setw(pw) << y_pos[j] << ' '; + for (auto i=0; i=0; ++j) + { + *fp << std::setw(pw) << y_pos[j] << ' '; + for (auto i=0; i(fp.get())->str()); + return Status::success; } From 997d07d3511436ffc000229d8ecfef80c75e5f21 Mon Sep 17 00:00:00 2001 From: Fernando Date: Tue, 27 Aug 2024 12:29:51 -0300 Subject: [PATCH 19/28] KICKTAB.ENH: Add new methods to allow loading kicktable from string. --- include/trackcpp/kicktable.h | 4 +- src/kicktable.cpp | 89 +++++++++++++++++++++++------------- 2 files changed, 59 insertions(+), 34 deletions(-) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index 27561fd..6a1de4c 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -26,8 +26,6 @@ class Kicktable { - // Kicktable: assumes x_kick and y_kick tables sampled on the same regular (x,y)-point grid. - public: std::string filename; @@ -47,6 +45,7 @@ class Kicktable { Kicktable(const Kicktable &) = default; Status::type load_from_file(const std::string& filename_); + Status::type load_from_file(std::string& filename_, bool file_flag=true); Status::type save_to_file( std::string& filename_, const std::string author_name = " ", @@ -135,6 +134,7 @@ class Kicktable { const double length=1 ); static int add_kicktable(const std::string& filename); + static int add_kicktable(std::string& filename, bool file_flag=true); static int add_kicktable(const Kicktable &new_kicktable); static Kicktable& get_kicktable(const int& kicktable_idx) { diff --git a/src/kicktable.cpp b/src/kicktable.cpp index 49265a6..2b9e8e5 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -16,9 +16,11 @@ #include #include -#include +#include #include +#include #include +#include std::vector Kicktable::kicktable_list; @@ -30,7 +32,6 @@ Kicktable::Kicktable(const std::string& filename_) : this->load_from_file(filename_); } - Kicktable::Kicktable( const std::vector& x_pos, const std::vector& y_pos, @@ -62,30 +63,49 @@ bool Kicktable::is_valid_kicktable() const return true; } -Status::type Kicktable::load_from_file(const std::string& filename_) { +Status::type Kicktable::load_from_file(const std::string& filename_) +{ + std::string fname = filename_; + return load_from_file(fname, true); +} - std::ifstream fp(filename_); - if (fp.fail()) { - std::cout << "Could not find kicktable file " << filename_ << "!" << std::endl; - return Status::file_not_found; +Status::type Kicktable::load_from_file( + std::string& filename_, + const bool file_flag +) +{ + // done with the help of chatgpt: + std::unique_ptr fp; + if (file_flag) + { + fp = std::make_unique(filename_); + if (!fp->good()) + { + std::cout << "Could not find kicktable file "; + std::cout << filename_ << "!" << std::endl; + return Status::file_not_found; + } + filename = filename_; } - filename = filename_; + else + fp = std::make_unique(filename_); + std::string str; unsigned int x_nrpts, y_nrpts; // HEADER - getline(fp, str); // name of kicktable line - getline(fp, str); // author line - getline(fp, str); // label 'ID length[m]' - fp >> length; // length of element - getline(fp, str); // advances to new line - getline(fp, str); // label 'number of horizontal points' - fp >> x_nrpts; // number of horizontal points - getline(fp, str); // advances to new line - getline(fp, str); // label 'number of vertical points' - fp >> y_nrpts; // number of vertical points - getline(fp, str); // advances to new line + getline(*fp, str); // author line + getline(*fp, str); // empty line + getline(*fp, str); // label 'ID length[m]' + *fp >> length; // length of element + getline(*fp, str); // advances to new line + getline(*fp, str); // label 'number of horizontal points' + *fp >> x_nrpts; // number of horizontal points + getline(*fp, str); // advances to new line + getline(*fp, str); // label 'number of vertical points' + *fp >> y_nrpts; // number of vertical points + getline(*fp, str); // advances to new line x_kick.resize(x_nrpts * y_nrpts, 0); y_kick.resize(x_nrpts * y_nrpts, 0); @@ -94,28 +114,28 @@ Status::type Kicktable::load_from_file(const std::string& filename_) { y_pos.resize(y_nrpts, 0); // HORIZONTAL KICK TABLE - getline(fp, str); // label 'Horizontal KickTable in T^2.m^2' - getline(fp, str); // label 'START' - for(unsigned int i=0; i> x_pos[i];} + getline(*fp, str); // label 'Horizontal KickTable in T^2.m^2' + getline(*fp, str); // label 'START' + for(unsigned int i=0; i> x_pos[i];} for(int j=y_nrpts-1; j>=0; --j) { - fp >> y_pos[j]; + *fp >> y_pos[j]; for(unsigned int i=0; i> x_kick[get_idx(i, j)]; + *fp >> x_kick[get_idx(i, j)]; } - getline(fp, str); // advances to new line + getline(*fp, str); // advances to new line // VERTICAL KICK TABLE - getline(fp, str); // label 'Vertical KickTable in T^2.m^2' - getline(fp, str); // label 'START' - getline(fp, str); // horizontal position. Already set from x_kick + getline(*fp, str); // label 'Vertical KickTable in T^2.m^2' + getline(*fp, str); // label 'START' + getline(*fp, str); // horizontal position. Already set from x_kick for(int j=y_nrpts-1; j>=0; --j) { - double posy; fp >> posy; + double posy; *fp >> posy; for(unsigned int i=0; i> y_kick[get_idx(i, j)]; + *fp >> y_kick[get_idx(i, j)]; } // invert tables, if necessary - if (y_pos.size() > 1 and y_pos[1] > y_pos[0]) { + if (y_pos.size() > 1 && y_pos[1] > y_pos[0]) { for(unsigned int i=0; i Date: Tue, 27 Aug 2024 12:33:47 -0300 Subject: [PATCH 20/28] FLATFILE.MNT: Change style of some methods. --- src/flat_file.cpp | 40 ++++++++++++++++++++++------------------ 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/src/flat_file.cpp b/src/flat_file.cpp index 42cbc3a..9dbbd08 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -45,38 +45,42 @@ static Status::type read_flat_file_trackcpp(std::istream&, Accelerator& accelera // -- implementation of API -- -Status::type read_flat_file(std::string& filename, Accelerator& accelerator, bool file_flag) { - if (file_flag) { +Status::type read_flat_file(std::string& filename, Accelerator& accelerator, bool file_flag) +{ + if (file_flag) + { std::ifstream fp(filename.c_str()); - if (fp.fail()) return Status::file_not_found; + if (!fp.good()) + return Status::file_not_found; read_flat_file_trackcpp(fp, accelerator); - fp.close(); - return Status::success; - } else { + } + else + { std::stringstream fp; fp.str(filename); read_flat_file_trackcpp(fp, accelerator); - //filename = fp.str(); - return Status::success; } - + return Status::success; } -Status::type write_flat_file(std::string& filename, const Accelerator& accelerator, bool file_flag) { - - if (file_flag) { +Status::type write_flat_file( + std::string& filename, const Accelerator& accelerator, bool file_flag +) +{ + if (file_flag) + { std::ofstream fp(filename.c_str()); - if (fp.fail()) return Status::file_not_found; + if (!fp.good()) + return Status::file_not_found; write_flat_file_trackcpp(fp, accelerator); - fp.close(); - return Status::success; - } else { + } + else + { std::stringstream fp; write_flat_file_trackcpp(fp, accelerator); filename = fp.str(); // very costly. no other way?! - return Status::success; } - + return Status::success; } void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) { From 8b74b81209868315418a3c7b1c2b7e6cfb4d2d59 Mon Sep 17 00:00:00 2001 From: Fernando Date: Tue, 27 Aug 2024 13:36:18 -0300 Subject: [PATCH 21/28] MNT: change back to C++17. --- INSTALL.md | 4 ++-- Makefile | 4 ++-- python_package/Makefile | 2 +- tests/Makefile | 4 ++-- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/INSTALL.md b/INSTALL.md index 42189a8..aac7c55 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -19,7 +19,7 @@ INSTRUCTIONS REQUIREMENTS -- a C++ compiler with C++20 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++20 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) diff --git a/Makefile b/Makefile index 45282f6..69c97ad 100644 --- a/Makefile +++ b/Makefile @@ -36,8 +36,8 @@ CC = gcc CXX = g++ AR = ar MACHINE = -m64 -OPT_FLAG = -O3 -std=c++20 -fPIC -DBG_FLAG = -O0 -g3 -std=c++20 -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 \ diff --git a/python_package/Makefile b/python_package/Makefile index 824f6fc..fa3a52f 100644 --- a/python_package/Makefile +++ b/python_package/Makefile @@ -25,7 +25,7 @@ NUMPY_VERSION := $(shell $(PYTHON) -c "import numpy; print(numpy.__version__)") OBJECTS = $(WRAPPEROBJS) $(INTERFACEOBJS) -CPPFLAGS = -std=c++20 -fPIC $(OPT_FLAG) +CPPFLAGS = -std=c++17 -fPIC $(OPT_FLAG) LIBS = $(shell gsl-config --libs) LIBS += -L../build -ltrackcpp INC = $(shell gsl-config --cflags) diff --git a/tests/Makefile b/tests/Makefile index c5f8e6f..98792bd 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -31,8 +31,8 @@ CC = gcc CXX = g++ AR = ar MACHINE = -m64 -OPT_FLAG = -O3 -std=c++20 -fPIC -DBG_FLAG = -O0 -g3 -std=c++20 -fPIC +OPT_FLAG = -O3 -std=c++17 -fPIC +DBG_FLAG = -O0 -g3 -std=c++17 -fPIC ARFLAGS = rcs DFLAGS = -DVERSION=$(VERSION) From 4212582f35cb03e5f0b24efd319857b8a4adbd9e Mon Sep 17 00:00:00 2001 From: Fernando Date: Tue, 27 Aug 2024 13:43:33 -0300 Subject: [PATCH 22/28] PYTHON.MNT: Remove unnecessary wrappers. --- python_package/interface.cpp | 122 +++++++++-------------------------- python_package/interface.h | 25 ------- 2 files changed, 30 insertions(+), 117 deletions(-) diff --git a/python_package/interface.cpp b/python_package/interface.cpp index caa18db..c3468b4 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -22,8 +22,8 @@ Status::type track_elementpass_wrapper ( const Accelerator& accelerator, const Element& el, double *pos, int n1, int n2 -) { - +) +{ std::vector > post; post.reserve(n2); @@ -48,8 +48,8 @@ Status::type track_linepass_wrapper( double *orig_pos, int ni1, int ni2, double *pos, int n1, int n2, LinePassArgs& args -) { - +) +{ std::vector > post; std::vector > orig_post; @@ -97,8 +97,8 @@ Status::type track_ringpass_wrapper ( double *orig_pos, int ni1, int ni2, double *pos, int n1, int n2, RingPassArgs& args -) { - +) +{ std::vector > post; std::vector > orig_post; @@ -143,12 +143,13 @@ Status::type track_ringpass_wrapper ( } Status::type calc_twiss_wrapper ( - Accelerator& accelerator, - const Pos& fixed_point, - Matrix& m66, - double *twiss, int n1, int n2, - Twiss twiss0) { - + Accelerator& accelerator, + const Pos& fixed_point, + Matrix& m66, + double *twiss, int n1, int n2, + Twiss twiss0 +) +{ std::vector > post; std::vector > orig_post; @@ -181,96 +182,31 @@ Status::type calc_twiss_wrapper ( return status; } -Element marker_wrapper(const std::string &fam_name_) { - return Element::marker(fam_name_); -} - -Element bpm_wrapper(const std::string &fam_name_) { - return Element::bpm(fam_name_); -} - -Element hcorrector_wrapper(const std::string &fam_name_, const double &length_, const double &hkick_) { - return Element::hcorrector(fam_name_, length_, hkick_); -} - -Element vcorrector_wrapper(const std::string &fam_name_, const double &length_, const double &vkick_) { - return Element::vcorrector(fam_name_, length_, vkick_); -} - -Element corrector_wrapper(const std::string &fam_name_, const double &length_, const double &hkick_, const double &vkick_) { - return Element::corrector(fam_name_, length_, hkick_, vkick_); -} - -Element drift_wrapper(const std::string &fam_name_, const double &length_) { - return Element::drift(fam_name_, length_); -} - -Element drift_g2l_wrapper(const std::string &fam_name_, const double &length_) { - return Element::drift_g2l(fam_name_, length_); -} - -Element matrix_wrapper(const std::string &fam_name_, const double &length_) { - return Element::matrix(fam_name_, length_); -} -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_, - const std::vector& polynom_a_, const std::vector& polynom_b_, - const double& K_, const double& S_) { - return Element::rbend(fam_name_, length_, angle_, angle_in_, angle_out_, gap_, fint_in_, fint_out_, polynom_a_, polynom_b_, K_, S_); -} - -Element quadrupole_wrapper(const std::string &fam_name_, const double &length_, const double &K_, const int nr_steps_) { - return Element::quadrupole(fam_name_, length_, K_, nr_steps_); -} - -Element sextupole_wrapper(const std::string &fam_name_, const double &length_, const double &S_, const int nr_steps_) { - return Element::sextupole(fam_name_, length_, S_, nr_steps_); -} - -Element rfcavity_wrapper(const std::string &fam_name_, const double &length_, const double &frequency_, const double &voltage_, const double &phase_lag_) { - return Element::rfcavity(fam_name_, length_, frequency_, voltage_, phase_lag_); -} - -Element kickmap_wrapper(const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_, const double& rescale_length_, const double& rescale_kicks_) { - 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) { +Status::type read_flat_file_wrapper( + String& fname, Accelerator& accelerator, bool file_flag +) +{ return read_flat_file(fname.data, accelerator, file_flag); } -Status::type write_flat_file_wrapper(String& fname, const Accelerator& accelerator, bool file_flag) { - return write_flat_file(fname.data, accelerator, file_flag); -} - -Status::type kicktable_getkicks_wrapper( - const int& kicktable_idx, - const double& rx, - const double& ry, - double& hkick__, - double& vkick__ +Status::type write_flat_file_wrapper( + String& fname, const Accelerator& accelerator, bool file_flag ) { - if (not Kicktable::is_valid_kicktable_index(kicktable_idx)) - return Status::kicktable_not_defined; - const Kicktable &kicktable = Kicktable::get_kicktable(kicktable_idx); - return kicktable.getkicks(rx, ry, hkick__, vkick__); + return write_flat_file(fname.data, accelerator, file_flag); } + Status::type track_findm66_wrapper( Accelerator& accelerator, const Pos& fixed_point, double *cumul_tm, int n1_tm, int n2_tm, int n3_tm, double *m66, int n1_m66, int n2_m66, Pos& v0, - std::vector& indices) { - + std::vector& indices +) +{ std::vector vec_tm; Matrix vec_m66; @@ -294,8 +230,9 @@ Status::type track_diffusionmatrix_wrapper( const Accelerator& accelerator, const Pos& fixed_point, double *cumul_tm, int n1_tm, int n2_tm, int n3_tm, - double *bdiffmats, int n1_bd, int n2_bd, int n3_bd){ - + double *bdiffmats, int n1_bd, int n2_bd, int n3_bd +) +{ std::vector vec_tm; std::vector vec_bd; @@ -326,8 +263,9 @@ void naff_general_wrapper( bool is_real, int nr_ff, int win, double *ff_out, int n1_ff_out, int n2_ff_out, double *re_out, int n1_re_out, int n2_re_out, - double *im_out, int n1_im_out, int n2_im_out){ - + double *im_out, int n1_im_out, int n2_im_out +) +{ std::vector vec_re_in; std::vector vec_im_in; std::vector vec_ff_out; diff --git a/python_package/interface.h b/python_package/interface.h index 32f3a2b..57acf5c 100644 --- a/python_package/interface.h +++ b/python_package/interface.h @@ -79,33 +79,8 @@ Status::type calc_twiss_wrapper ( Twiss twiss0 ); -Element marker_wrapper(const std::string& fam_name_); -Element bpm_wrapper(const std::string& fam_name_); -Element hcorrector_wrapper(const std::string& fam_name_, const double& length_, const double& hkick_); -Element vcorrector_wrapper(const std::string& fam_name_, const double& length_, const double& vkick_); -Element corrector_wrapper(const std::string& fam_name_, const double& length_, const double& hkick_, const double& vkick_); -Element drift_wrapper(const std::string& fam_name_, const double& length_); -Element drift_g2l_wrapper(const std::string& fam_name_, const double& length_); -Element matrix_wrapper(const std::string& fam_name_, const double& length_); -Element quadrupole_wrapper(const std::string& fam_name_, const double& length_, const double& K_, const int nr_steps_ = 10); -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_, - const std::vector& polynom_a_, const std::vector& polynom_b_, - const double& K_, const double& S_); Status::type write_flat_file_wrapper(String& fname, const Accelerator& accelerator, bool file_flag = true); Status::type read_flat_file_wrapper(String& fname, Accelerator& accelerator, bool file_flag = true); -Status::type kicktable_getkicks_wrapper( - const int& kicktable_idx, - const double& rx, - const double& ry, - double& hkick__, - double& vkick__ -); Status::type track_findm66_wrapper( Accelerator& accelerator, From 5fc73a448ec216795899db46baa5528e47dfb8f5 Mon Sep 17 00:00:00 2001 From: Fernando Date: Tue, 27 Aug 2024 13:44:11 -0300 Subject: [PATCH 23/28] PYTHON.ENH: Add wrapper for Kicktable::getkicks. --- include/trackcpp/kicktable.h | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index 6a1de4c..bc41a9c 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -126,6 +126,13 @@ class Kicktable { return getkicks_bilinear(rx, ry, hkick, vkick); } + Status::type getkicks( + const double& rx, const double& ry, double& hkick__, double& vkick__ + ) const + { + return getkicks(rx, ry, hkick__, vkick__); + } + static int add_kicktable( const std::vector& x_pos, const std::vector& x_kick, From a784e0339931cb4a5a8b37f9f09c8dacc860bf1c Mon Sep 17 00:00:00 2001 From: Fernando Date: Tue, 27 Aug 2024 13:45:43 -0300 Subject: [PATCH 24/28] PYTHON.DOC: Update comment on getkicks. --- python_package/trackcpp.i | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/python_package/trackcpp.i b/python_package/trackcpp.i index 3c2889f..5c447c4 100644 --- a/python_package/trackcpp.i +++ b/python_package/trackcpp.i @@ -66,8 +66,10 @@ double get_double_max() { %} -%extend std::vector { - double* data_() { +%extend std::vector +{ + double* data_() + { return $self->data(); } } @@ -110,7 +112,7 @@ double get_double_max() { %apply (double* INPLACE_ARRAY2, int DIM1, int DIM2 ) { (double *im_out, int n1_im_out, int n2_im_out)} -// For kicktable_getkicks_wrapper +// For Kicktable::getkicks %apply double *OUTPUT { double &hkick__, double &vkick__ }; From f62cb4929899c2582b09fe8190e2ba24e81ba8be Mon Sep 17 00:00:00 2001 From: Fernando Date: Tue, 27 Aug 2024 17:00:03 -0300 Subject: [PATCH 25/28] KICKTAB.ENH: Fix some bugs and make interface more suitable for python package. basic functions were tested. --- include/trackcpp/kicktable.h | 18 +++-- python_package/trackcpp.i | 1 - src/kicktable.cpp | 125 ++++++++++++++++------------------- 3 files changed, 68 insertions(+), 76 deletions(-) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index bc41a9c..e3876a9 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -44,13 +44,16 @@ class Kicktable { Kicktable(const std::string& filename_ = ""); Kicktable(const Kicktable &) = default; - Status::type load_from_file(const std::string& filename_); - Status::type load_from_file(std::string& filename_, bool file_flag=true); + 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( - std::string& filename_, - const std::string author_name = " ", - bool file_flag = true + 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 @@ -71,6 +74,8 @@ class Kicktable { 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(x);} + unsigned int get_iy(const double& y) const {return get_iy(y);} template Status::type getkicks_bilinear( const T& rx, const T& ry, T& hkick, T& vkick @@ -140,8 +145,7 @@ class Kicktable { const std::vector& y_kick, const double length=1 ); - static int add_kicktable(const std::string& filename); - static int add_kicktable(std::string& filename, bool file_flag=true); + static int add_kicktable(const std::string& filename, bool file_flag=true); static int add_kicktable(const Kicktable &new_kicktable); static Kicktable& get_kicktable(const int& kicktable_idx) { diff --git a/python_package/trackcpp.i b/python_package/trackcpp.i index 5c447c4..851c96d 100644 --- a/python_package/trackcpp.i +++ b/python_package/trackcpp.i @@ -115,7 +115,6 @@ double get_double_max() { // For Kicktable::getkicks %apply double *OUTPUT { double &hkick__, double &vkick__ }; - %include "../include/trackcpp/elements.h" %include "../include/trackcpp/kicktable.h" %include "../include/trackcpp/auxiliary.h" diff --git a/src/kicktable.cpp b/src/kicktable.cpp index 2b9e8e5..81c6ce9 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -63,14 +63,8 @@ bool Kicktable::is_valid_kicktable() const return true; } -Status::type Kicktable::load_from_file(const std::string& filename_) -{ - std::string fname = filename_; - return load_from_file(fname, true); -} - Status::type Kicktable::load_from_file( - std::string& filename_, + const std::string& filename_, const bool file_flag ) { @@ -90,7 +84,6 @@ Status::type Kicktable::load_from_file( else fp = std::make_unique(filename_); - std::string str; unsigned int x_nrpts, y_nrpts; @@ -135,7 +128,7 @@ Status::type Kicktable::load_from_file( } // invert tables, if necessary - if (y_pos.size() > 1 && y_pos[1] > y_pos[0]) { + if (y_pos.size() > 1 && y_pos[1] < y_pos[0]) { for(unsigned int i=0; i fp; - if (file_flag) - { - fp = std::make_unique(filename.c_str()); - if (!fp->good()) - return Status::file_not_found; - } - else - fp = std::make_unique(); - - fp->setf( + fp.setf( std::ios_base::left | std::ios_base::scientific | - std::ios_base::uppercase + std::ios_base::uppercase | + std::ios_base::showpos ); - fp->precision(np); + fp.precision(np); // HEADER - *fp << "# Author: " << author_name << '\n'; - *fp << "#" << '\n'; - *fp << "# Total Length of Longitudinal Interval [m]" << '\n'; - *fp << length << '\n'; - *fp << "# Number of Horizontal Points" << '\n'; - *fp << x_pos.size() << '\n'; - *fp << "# Number of Vertical Points" << '\n'; - *fp << y_pos.size() << '\n'; + fp << "# Author: " << author_name << '\n'; + fp << "#" << '\n'; + fp << "# Total Length of Longitudinal Interval [m]" << '\n'; + fp << length << '\n'; + fp << "# Number of Horizontal Points" << '\n'; + fp << x_pos.size() << '\n'; + fp << "# Number of Vertical Points" << '\n'; + fp << y_pos.size() << '\n'; // x_kick: - *fp << "# Total Horizontal 2nd Order Kick [T2m2]" << '\n'; - *fp << "START" << '\n'; - *fp << std::setw(pw) << ' '; + fp << "# Total Horizontal 2nd Order Kick [T2m2]" << '\n'; + fp << "START" << '\n'; + fp << std::setw(pw) << ' ' << ' '; for (auto xi: x_pos) - *fp << std::setw(pw) << xi << ' '; - *fp << '\n'; - for (auto j=y_pos.size()-1; j>=0; ++j) + fp << std::setw(pw) << xi << ' '; + fp << '\n'; + for (int j=y_pos.size()-1; j>=0; --j) { - *fp << std::setw(pw) << y_pos[j] << ' '; + fp << std::setw(pw) << y_pos[j] << ' '; for (auto i=0; i=0; ++j) + fp << std::setw(pw) << xi << ' '; + fp << '\n'; + for (int j=y_pos.size()-1; j>=0; --j) { - *fp << std::setw(pw) << y_pos[j] << ' '; + fp << std::setw(pw) << y_pos[j] << ' '; for (auto i=0; i(fp.get())->str()); return Status::success; +} +Status::type Kicktable::save_to_file( + const std::string filename_, + const std::string author_name +) +{ + std::ofstream fp(filename_.c_str()); + if (!fp.good()) + return Status::file_not_found; + return _dump_to_stream(fp, author_name); +} + +std::string Kicktable::save_to_string(const std::string author_name) +{ + std::stringstream fp; + _dump_to_stream(fp, author_name); + std::string stg = std::move(fp.str()); + return stg; } int Kicktable::add_kicktable( @@ -237,14 +233,7 @@ int Kicktable::add_kicktable( return Kicktable::add_kicktable(new_kicktable); } - -int Kicktable::add_kicktable(const std::string& filename) -{ - // loads a new kicktable from file and inserts it into vector of kicktables - Kicktable new_kicktable(filename); - return Kicktable::add_kicktable(new_kicktable); -} -int Kicktable::add_kicktable(std::string& filename, bool file_flag) +int Kicktable::add_kicktable(const std::string& filename, bool file_flag) { Kicktable new_kicktable; new_kicktable.load_from_file(filename, file_flag); From df77a8979b095c2637f1a6099f6fb7b9f8738925 Mon Sep 17 00:00:00 2001 From: Fernando Date: Wed, 28 Aug 2024 17:24:09 -0300 Subject: [PATCH 26/28] FLATFILE.BUG: increase width of parameters field names and add space to end of name to ensure values are not glued to them. --- src/flat_file.cpp | 144 ++++++++++++++++++++++++++-------------------- 1 file changed, 81 insertions(+), 63 deletions(-) diff --git a/src/flat_file.cpp b/src/flat_file.cpp index 9dbbd08..44ba751 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -24,7 +24,7 @@ #include static const int hw = 18; // header field width -static const int pw = 16; // parameter field width +static const int pw = 17; // parameter field width static const int np = 17; // number precision static bool read_boolean_string(std::istringstream& ss); static int process_rad_property(std::istringstream& ss); @@ -88,12 +88,12 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) fp.setf(std::ios_base::left | std::ios_base::scientific | std::ios_base::uppercase); fp.precision(np); - fp << std::setw(hw) << "% energy" << accelerator.energy << " eV\n"; - fp << std::setw(hw) << "% harmonic_number" << accelerator.harmonic_number << "\n"; - fp << std::setw(hw) << "% cavity_on" << get_boolean_string(accelerator.cavity_on) << "\n"; - fp << std::setw(hw) << "% radiation_on" << accelerator.radiation_on << "\n"; - fp << std::setw(hw) << "% vchamber_on" << get_boolean_string(accelerator.vchamber_on) << "\n"; - fp << std::setw(hw) << "% lattice_version" << accelerator.lattice_version << "\n"; + fp << std::setw(hw) << "% energy " << accelerator.energy << " eV\n"; + fp << std::setw(hw) << "% harmonic_number " << accelerator.harmonic_number << "\n"; + fp << std::setw(hw) << "% cavity_on " << get_boolean_string(accelerator.cavity_on) << "\n"; + fp << std::setw(hw) << "% radiation_on " << accelerator.radiation_on << "\n"; + fp << std::setw(hw) << "% vchamber_on " << get_boolean_string(accelerator.vchamber_on) << "\n"; + fp << std::setw(hw) << "% lattice_version " << accelerator.lattice_version << "\n"; fp << '\n'; fp.setf(std::ios_base::showpos); @@ -101,78 +101,96 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) for (auto i=0; i Date: Tue, 10 Sep 2024 16:06:18 -0300 Subject: [PATCH 27/28] BUG: several bug fixes. --- include/trackcpp/kicktable.h | 45 ++++++++------ include/trackcpp/passmethods.hpp | 32 ++++++---- python_package/trackcpp.i | 103 ++++++++++++++++++++++--------- src/kicktable.cpp | 4 +- 4 files changed, 122 insertions(+), 62 deletions(-) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index e3876a9..6783f7d 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -36,8 +36,8 @@ class Kicktable { Kicktable( const std::vector& x_pos, - const std::vector& x_kick, const std::vector& y_pos, + const std::vector& x_kick, const std::vector& y_kick, const double length = 1 ); @@ -81,43 +81,47 @@ class Kicktable { 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); - double denom = (x2 - x1) * (y2 - y1); + 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[get_idx(ix, iy)]; - const double f12 = x_kick[get_idx(ix, iyp1)]; - const double f21 = x_kick[get_idx(ixp1, iy)]; - const double f22 = x_kick[get_idx(ixp1, iyp1)]; - hkick = - f11 * (x2 - rx) * (y2 - ry) + f21 * (rx - x1) * (y2 - ry) + - f12 * (x2 - rx) * (ry - y1) + f22 * (rx - x1) * (ry - y1); + 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[get_idx(ix, iy)]; - const double f12 = y_kick[get_idx(ix, iyp1)]; - const double f21 = y_kick[get_idx(ixp1, iy)]; - const double f22 = y_kick[get_idx(ixp1, iyp1)]; - vkick = - f11 * (x2 - rx) * (y2 - ry) + f21 * (rx - x1) * (y2 - ry) + - f12 * (x2 - rx) * (ry - y1) + f22 * (rx - x1) * (ry - y1); + 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; } @@ -140,12 +144,12 @@ class Kicktable { static int add_kicktable( const std::vector& x_pos, - const std::vector& x_kick, const std::vector& y_pos, + const std::vector& x_kick, const std::vector& y_kick, const double length=1 ); - static int add_kicktable(const std::string& filename, bool file_flag=true); + static int add_kicktable(const std::string filename); static int add_kicktable(const Kicktable &new_kicktable); static Kicktable& get_kicktable(const int& kicktable_idx) { @@ -155,9 +159,10 @@ class Kicktable { } static bool is_valid_kicktable_index(const int idx) { - return ((idx>0) & (idx < kicktable_list.size())); + return ((idx>=0) & (idx < kicktable_list.size())); } static void clear_kicktables(); + static size_t get_kicktable_list_size(){return kicktable_list.size();} bool operator==(const Kicktable& o) const; bool operator!=(const Kicktable& o) const { return !(*this == o); } diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 0763b54..11702a0 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -146,19 +146,27 @@ void kickpolythinkick( T hkick = 0; T vkick = 0; T rx_p = 1; - unsigned int k = 0; + unsigned int k0 = 0; + const size_t siz = polyx.size(); + // std::cout << "call" << std::endl; 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; + T ry_p = 1; + unsigned int k = k0; + for (auto j = i; j <= order; ++j) + { + if (k >= siz) break; + // std::cout << " o: " << order << " i: " << i << " j: " << j; + // std::cout << " k0: " << k0 << " k: " << k << " rx: " << rx_p; + // std::cout << " ry: " << ry_p << " px: " << polyx[k]; + // std::cout << " py: " << polyy[k] << std::endl; + hkick += polyx[k] * rx_p * ry_p; + vkick += polyy[k] * rx_p * ry_p; + ry_p *= pos.ry; + k += j + 2; + } + k0 += i + 1; + rx_p *= pos.rx; } // According to Ellaune's theory of kick maps: // https://accelconf.web.cern.ch/e92/PDF/EPAC1992_0661.PDF @@ -562,7 +570,7 @@ Status::type pm_kickpoly_pass( // 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; + const unsigned int order = (sqrt(1 + 8*elem.polynom_kickx.size()) - 1)/2; global_2_local(pos, elem); for(unsigned int i=0; i; @@ -54,16 +84,15 @@ namespace std { } %inline %{ -double c_array_get(double* v, int i) { - return v[i]; -} -void c_array_set(double* v, int i, double x) { - v[i] = x; -} -double get_double_max() { - return DBL_MAX; -} - + double c_array_get(double* v, int i) { + return v[i]; + } + void c_array_set(double* v, int i, double x) { + v[i] = x; + } + double get_double_max() { + return DBL_MAX; + } %} %extend std::vector @@ -91,26 +120,44 @@ double get_double_max() { (double* twiss, int n1, int n2)} //For find_m66 and diffusion_matrix -%apply (double* INPLACE_ARRAY3, int DIM1, int DIM2, int DIM3 ) { - (double *cumul_tm, int n1_tm, int n2_tm, int n3_tm)} -%apply (double* INPLACE_ARRAY3, int DIM1, int DIM2, int DIM3 ) { - (double *elem_tm, int n1_tm, int n2_tm, int n3_tm)} -%apply (double* INPLACE_ARRAY3, int DIM1, int DIM2, int DIM3 ) { - (double *bdiffmats, int n1_bd, int n2_bd, int n3_bd)} -%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2 ) { - (double *m66, int n1_m66, int n2_m66)} +%apply (double* INPLACE_ARRAY3, int DIM1, int DIM2, int DIM3 ) +{ + (double *cumul_tm, int n1_tm, int n2_tm, int n3_tm) +} +%apply (double* INPLACE_ARRAY3, int DIM1, int DIM2, int DIM3 ) +{ + (double *elem_tm, int n1_tm, int n2_tm, int n3_tm) +} +%apply (double* INPLACE_ARRAY3, int DIM1, int DIM2, int DIM3 ) +{ + (double *bdiffmats, int n1_bd, int n2_bd, int n3_bd) +} +%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2 ) +{ + (double *m66, int n1_m66, int n2_m66) +} // For naff_general -%apply (double* IN_ARRAY2, int DIM1, int DIM2 ){ - (double* re_in, int n1_re_in, int n2_re_in)} -%apply (double* IN_ARRAY2, int DIM1, int DIM2 ){ - (double* im_in, int n1_im_in, int n2_im_in)} -%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2 ) { - (double *ff_out, int n1_ff_out, int n2_ff_out)} -%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2 ) { - (double *re_out, int n1_re_out, int n2_re_out)} -%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2 ) { - (double *im_out, int n1_im_out, int n2_im_out)} +%apply (double* IN_ARRAY2, int DIM1, int DIM2 ) +{ + (double* re_in, int n1_re_in, int n2_re_in) +} +%apply (double* IN_ARRAY2, int DIM1, int DIM2 ) +{ + (double* im_in, int n1_im_in, int n2_im_in) +} +%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2 ) +{ + (double *ff_out, int n1_ff_out, int n2_ff_out) +} +%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2 ) +{ + (double *re_out, int n1_re_out, int n2_re_out) +} +%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2 ) +{ + (double *im_out, int n1_im_out, int n2_im_out) +} // For Kicktable::getkicks %apply double *OUTPUT { double &hkick__, double &vkick__ }; diff --git a/src/kicktable.cpp b/src/kicktable.cpp index 81c6ce9..d0bc095 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -233,10 +233,10 @@ int Kicktable::add_kicktable( return Kicktable::add_kicktable(new_kicktable); } -int Kicktable::add_kicktable(const std::string& filename, bool file_flag) +int Kicktable::add_kicktable(const std::string filename) { Kicktable new_kicktable; - new_kicktable.load_from_file(filename, file_flag); + new_kicktable.load_from_file(filename, true); return Kicktable::add_kicktable(new_kicktable); } From 66193280ca24dad07b80b270e37a03360fb59785 Mon Sep 17 00:00:00 2001 From: Fernando Date: Wed, 11 Sep 2024 17:25:44 -0300 Subject: [PATCH 28/28] KICKTAB.BUG: Use std::list instead of std::vector as container for kicktable_list to fix broken references. --- include/trackcpp/kicktable.h | 10 ++++++--- python_package/trackcpp.i | 3 ++- src/kicktable.cpp | 40 ++++++++++++++++++++---------------- 3 files changed, 31 insertions(+), 22 deletions(-) diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index 6783f7d..c334330 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -21,6 +21,7 @@ // #include "../alglib/interpolation.h" #include #include +#include #include @@ -32,7 +33,7 @@ class Kicktable { double length; std::vector x_pos, y_pos; std::vector x_kick, y_kick; - static std::vector kicktable_list; + static std::list kicktable_list; Kicktable( const std::vector& x_pos, @@ -150,12 +151,15 @@ class Kicktable { const double length=1 ); static int add_kicktable(const std::string filename); - static int add_kicktable(const Kicktable &new_kicktable); + 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."); - return kicktable_list[kicktable_idx]; + + auto it = kicktable_list.begin(); + std::advance(it, kicktable_idx); + return *it; } static bool is_valid_kicktable_index(const int idx) { diff --git a/python_package/trackcpp.i b/python_package/trackcpp.i index 3a5989b..d490d8f 100644 --- a/python_package/trackcpp.i +++ b/python_package/trackcpp.i @@ -35,6 +35,7 @@ %include "carrays.i" %include "std_string.i" %include "std_vector.i" +%include "std_list.i" %include "stl.i" %include "typemaps.i" %include "std_except.i" @@ -80,7 +81,7 @@ namespace std { %template(CppDoubleMatrixVector) vector< vector< vector > >; %template(CppTwissVector) vector; %template(CppMatrixVector) vector< Matrix >; - %template(CppKicktableVector) vector< Kicktable >; + %template(CppKicktableVector) list< Kicktable >; } %inline %{ diff --git a/src/kicktable.cpp b/src/kicktable.cpp index d0bc095..8e324f3 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -23,7 +23,7 @@ #include -std::vector Kicktable::kicktable_list; +std::list Kicktable::kicktable_list; Kicktable::Kicktable(const std::string& filename_) : filename("") @@ -229,29 +229,33 @@ int Kicktable::add_kicktable( const double length ) { - Kicktable new_kicktable = Kicktable(x_pos, y_pos, x_kick, y_kick, length); - return Kicktable::add_kicktable(new_kicktable); + Kicktable new_ktab = Kicktable(x_pos, y_pos, x_kick, y_kick, length); + return Kicktable::add_kicktable(new_ktab); } int Kicktable::add_kicktable(const std::string filename) { - Kicktable new_kicktable; - new_kicktable.load_from_file(filename, true); - return Kicktable::add_kicktable(new_kicktable); + Kicktable new_ktab; + new_ktab.load_from_file(filename, true); + return Kicktable::add_kicktable(new_ktab); } -int Kicktable::add_kicktable(const Kicktable &new_kicktable) +int Kicktable::add_kicktable(const Kicktable &new_ktab) { - if (not new_kicktable.is_valid_kicktable()) + if (not new_ktab.is_valid_kicktable()) return -1; // looks through vector of kicktables... - for(int i=0; ilength != o.length) return false; - if (this->x_pos != o.x_pos) return false; - if (this->y_pos != o.y_pos) return false; - if (this->x_kick != o.x_kick) return false; - if (this->y_kick != o.y_kick) return false; + if (length != o.length) return false; + if (x_pos != o.x_pos) return false; + if (y_pos != o.y_pos) return false; + if (x_kick != o.x_kick) return false; + if (y_kick != o.y_kick) return false; return true; }