diff --git a/INSTALL.md b/INSTALL.md index 7b14689..aac7c55 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++17 support (tested with GCC 4.8) - pthread - blas - gsl (GNU Scientific Library) @@ -49,7 +49,7 @@ INSTRUCTIONS REQUIREMENTS -- a C++ compiler with C++11 support (tested with GCC 4.8) +- a C++ compiler with C++17 support (tested with GCC 4.8) - blas - gsl (GNU Scientific Library) - SWIG 2 (tested with SWIG 2.0.11) diff --git a/Makefile b/Makefile index d573057..69c97ad 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++17 -fPIC +DBG_FLAG = -O0 -g3 -std=c++17 -fPIC ARFLAGS = rcs DFLAGS = -DVERSION=$(VERSION) LIBSOURCES_CPP = lattice.cpp \ @@ -58,7 +58,7 @@ LIBSOURCES_CPP = lattice.cpp \ auxiliary.cpp BINSOURCES_CPP = exec.cpp \ tests.cpp \ - commands.cpp \ + commands.cpp AUXFILES = VERSION 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/elements.h b/include/trackcpp/elements.h index 94942d9..97f50ce 100644 --- a/include/trackcpp/elements.h +++ b/include/trackcpp/elements.h @@ -61,11 +61,13 @@ class Element { double frequency = 0; // [Hz] double voltage = 0; // [V] double phase_lag = 0; // [rad] - int kicktable_idx = -1; // index of kickmap object in kicktable_list + int kicktable_idx = -1; //index of kickmap object in kicktable_list double rescale_kicks = 1.0; // for kickmaps std::vector 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/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); diff --git a/include/trackcpp/kicktable.h b/include/trackcpp/kicktable.h index 9457287..c334330 100644 --- a/include/trackcpp/kicktable.h +++ b/include/trackcpp/kicktable.h @@ -18,108 +18,159 @@ #define _KICKTABLE_H #include "auxiliary.h" +// #include "../alglib/interpolation.h" #include #include +#include +#include class Kicktable { - // Kicktable: assumes x_kick and y_kick tables sampled on the same regular (x,y)-point grid. - public: std::string filename; double length; - unsigned int x_nrpts, y_nrpts; - double x_min, x_max; - double y_min, y_max; + std::vector x_pos, y_pos; std::vector x_kick, y_kick; - + static std::list kicktable_list; + + Kicktable( + const std::vector& x_pos, + const std::vector& y_pos, + const std::vector& x_kick, + 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 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) { - - const Kicktable* kicktable = &kicktable_list[kicktable_idx]; + Status::type load_from_file( + const std::string& filename_, bool file_flag=true + ); + Status::type _dump_to_stream( + std::ostream& fp, const std::string author_name = " " + ); + Status::type save_to_file( + const std::string filename_, const std::string author_name = " " + ); + std::string save_to_string(const std::string author_name = " "); + + bool is_valid_kicktable() const; + unsigned int get_idx(unsigned int ix, unsigned int iy) const + { + return iy * x_pos.size() + ix; + } + double get_x(unsigned int ix) const {return x_pos[ix];} + double get_y(unsigned int iy) const {return y_pos[iy];} + template unsigned int get_ix(const T& x) const + { + int idx; + idx = std::lower_bound(x_pos.begin(), x_pos.end(), x) - x_pos.begin() - 1; + return idx >= 0 ? (unsigned int) idx : 0u; + } + template unsigned int get_iy(const T& y) const + { + int idx; + idx = std::lower_bound(y_pos.begin(), y_pos.end(), y) - y_pos.begin() - 1; + return idx >= 0 ? (unsigned int) idx : 0u; + } + unsigned int get_ix(const double& x) const {return get_ix(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 + ) const + { + + // gets indices + unsigned int ix = get_ix(rx); + unsigned int iy = get_iy(ry); + unsigned int ixp1, iyp1; + if (ix >= x_pos.size()-1) ixp1 = ix--; + else ixp1 = ix + 1; + if (iy >= y_pos.size()-1) iyp1 = iy--; + else iyp1 = iy + 1; + const unsigned int i00 = get_idx(ix, iy); + const unsigned int i01 = get_idx(ix, iyp1); + const unsigned int i10 = get_idx(ixp1, iy); + const unsigned int i11 = get_idx(ixp1, iyp1); + + /* coordinates */ + const double x1 = get_x(ix); + const double x2 = get_x(ixp1); + const double y1 = get_y(iy); + const double y2 = get_y(iyp1); + const double denom = (x2 - x1) * (y2 - y1); + const T dx1 = rx - x1; + const T dx2 = x2 - rx; + const T dy1 = ry - y1; + const T dy2 = y2 - ry; + + // Bilinear interpolation + // https://en.wikipedia.org/wiki/Bilinear_interpolation + { /* hkick */ + const double& f11 = x_kick[i00]; + const double& f12 = x_kick[i01]; + const double& f21 = x_kick[i10]; + const double& f22 = x_kick[i11]; + hkick = (f11 * dx2 + f21 * dx1) * dy2 + (f12 * dx2 + f22 * dx1) * dy1; + hkick /= denom; + } + { /* vkick */ + const double& f11 = y_kick[i00]; + const double& f12 = y_kick[i01]; + const double& f21 = y_kick[i10]; + const double& f22 = y_kick[i11]; + vkick = (f11 * dx2 + f21 * dx1) * dy2 + (f12 * dx2 + f22 * dx1) * dy1; + vkick /= denom; + } + + return Status::success; + } - // checks x limits - const double& xmin = kicktable->x_min; - const double& xmax = kicktable->x_max; - if ((rx < xmin) or (rx > xmax)) { - //std::cout << "rx: " << double(rx) << std::endl; - hkick = nan(""); - return Status::kicktable_out_of_range; + template Status::type getkicks( + const T& rx, const T& ry, T& hkick, T& vkick + ) const + { + return getkicks_bilinear(rx, ry, hkick, vkick); } - // checks y limits - const double& ymin = kicktable->y_min; - const double& ymax = kicktable->y_max; - if ((ry < ymin) or (ry > ymax)) { - //std::cout << "ry: " << double(ry) << std::endl; - vkick = nan(""); - return Status::kicktable_out_of_range; + Status::type getkicks( + const double& rx, const double& ry, double& hkick__, double& vkick__ + ) const + { + return getkicks(rx, ry, hkick__, vkick__); } - // gets indices - const unsigned int ix = kicktable->get_ix(rx); - const unsigned int iy = kicktable->get_iy(ry); - - /* coordinates */ - const double x1 = kicktable->get_x(ix); - const double x2 = kicktable->get_x(ix+1); - const double y1 = kicktable->get_y(iy); - const double y2 = kicktable->get_y(iy+1); - - // Bilinear interpolation - https://en.wikipedia.org/wiki/Bilinear_interpolation - - { /* hkick */ - const double fq11 = kicktable->x_kick[kicktable->get_idx(ix+0, iy+0)]; - const double fq12 = kicktable->x_kick[kicktable->get_idx(ix+0, iy+1)]; - const double fq21 = kicktable->x_kick[kicktable->get_idx(ix+1, iy+0)]; - const double fq22 = kicktable->x_kick[kicktable->get_idx(ix+1, iy+1)]; - const T fR1 = ((x2 - rx) * fq11 + (rx - x1) * fq21) / (x2 - x1); - const T fR2 = ((x2 - rx) * fq12 + (rx - x1) * fq22) / (x2 - x1); - const T fP = ((y2 - ry) * fR1 + (ry - y1) * fR2 ) / (y2 - y1); - hkick = fP; + static int add_kicktable( + const std::vector& x_pos, + 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); + static int add_kicktable(const Kicktable &new_ktab); + static Kicktable& get_kicktable(const int& kicktable_idx) + { + if (not is_valid_kicktable_index(kicktable_idx)) + throw std::out_of_range("kicktable_idx is out of range."); + + auto it = kicktable_list.begin(); + std::advance(it, kicktable_idx); + return *it; } - { /* vkick */ - const double fq11 = kicktable->y_kick[kicktable->get_idx(ix+0, iy+0)]; - const double fq12 = kicktable->y_kick[kicktable->get_idx(ix+0, iy+1)]; - const double fq21 = kicktable->y_kick[kicktable->get_idx(ix+1, iy+0)]; - const double fq22 = kicktable->y_kick[kicktable->get_idx(ix+1, iy+1)]; - const T fR1 = ((x2 - rx) * fq11 + (rx - x1) * fq21) / (x2 - x1); - const T fR2 = ((x2 - rx) * fq12 + (rx - x1) * fq22) / (x2 - x1); - const T fP = ((y2 - ry) * fR1 + (ry - y1) * fR2 ) / (y2 - y1); - vkick = fP; + static bool is_valid_kicktable_index(const int idx) + { + return ((idx>=0) & (idx < kicktable_list.size())); } + static void clear_kicktables(); + static size_t get_kicktable_list_size(){return kicktable_list.size();} - return Status::success; -} + bool operator==(const Kicktable& o) const; + bool operator!=(const Kicktable& o) const { return !(*this == o); } -template -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/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/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..11702a0 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 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); - pos.px += rescale_kicks * hkick / (brho * brho) / nr_steps; - pos.py += rescale_kicks * vkick / (brho * brho) / nr_steps; + 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 + // + // 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 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 + // 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,60 @@ 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 k0 = 0; + const size_t siz = polyx.size(); + // std::cout << "call" << std::endl; + for (auto i = 0; i <= order; ++i) + { + 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 + // + // 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 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 + // 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) { @@ -442,15 +515,21 @@ 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); - const double brho = get_magnetic_rigidity(accelerator.energy); + double sl = elem.length / float(elem.nr_steps); + double brho2 = get_magnetic_rigidity(accelerator.energy); + brho2 *= brho2; global_2_local(pos, elem); if (elem.kicktable_idx < 0) { @@ -458,7 +537,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, + kicktable, + brho2, + elem.nr_steps, + elem.rescale_kicks + ); if (status != Status::success) return status; drift(pos, sl / 2); } @@ -468,6 +553,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) { diff --git a/include/trackcpp/tpsa.h b/include/trackcpp/tpsa.h index 580a2cf..bf141c7 100644 --- a/include/trackcpp/tpsa.h +++ b/include/trackcpp/tpsa.h @@ -46,14 +46,13 @@ // 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 #include #include #include +#include // Expression Templates: IMPLEMENTATION OF BINOMIALS COEFFICIENTS AND RELATED RELEVANT EXPRESSIONS @@ -375,6 +374,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; @@ -513,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; @@ -529,7 +563,6 @@ Tpsa sqrt(const Tpsa& a_) { return r; } -#include template Tpsa log(const Tpsa& a_) { Tpsa r; @@ -543,7 +576,6 @@ Tpsa log(const Tpsa& a_) { return r; } -#include template Tpsa cos(const Tpsa& a_) { Tpsa rc(1), rs; @@ -558,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; @@ -573,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; @@ -594,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; @@ -609,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])); @@ -640,7 +666,6 @@ Tpsa atan(const Tpsa& a_) { return r; } -#include template Tpsa asin(const Tpsa& a_) { @@ -655,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])); 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/python_package/Makefile b/python_package/Makefile index 06edce0..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++11 -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/python_package/interface.cpp b/python_package/interface.cpp index 383c7fe..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,73 +182,21 @@ 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_); -} - -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) { +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__) { - return kicktable_getkicks(kicktable_idx, rx, ry, hkick__, vkick__); -} Status::type track_findm66_wrapper( Accelerator& accelerator, @@ -255,8 +204,9 @@ Status::type track_findm66_wrapper( 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; @@ -280,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; @@ -312,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 e35e092..57acf5c 100644 --- a/python_package/interface.h +++ b/python_package/interface.h @@ -79,26 +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 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, diff --git a/python_package/trackcpp.i b/python_package/trackcpp.i index 766a189..d490d8f 100644 --- a/python_package/trackcpp.i +++ b/python_package/trackcpp.i @@ -35,8 +35,39 @@ %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" + + +// Define a custom exception in Python +%pythoncode { + class TrackingError(Exception): + pass +} + +%exception { + try + { + $action + } + catch (const std::out_of_range &e) + { + PyObject *ex = PyObject_GetAttrString( + PyImport_AddModule("trackcpp"), "TrackingError" + ); + if (ex != NULL) PyErr_SetString(ex, e.what()); + else PyErr_SetString(PyExc_IndexError, e.what()); + SWIG_fail; + } + catch (...) + { + PyErr_SetString(PyExc_RuntimeError, "Unknown exception occurred"); + SWIG_fail; + } +} + namespace std { %template(CppStringVector) vector; @@ -50,24 +81,25 @@ namespace std { %template(CppDoubleMatrixVector) vector< vector< vector > >; %template(CppTwissVector) vector; %template(CppMatrixVector) vector< Matrix >; - %template(CppKicktableVector) vector< Kicktable >; + %template(CppKicktableVector) list< Kicktable >; } %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 { - double* data_() { +%extend std::vector +{ + double* data_() + { return $self->data(); } } @@ -89,31 +121,48 @@ 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 +// 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/elements.cpp b/src/elements.cpp index 5d96fb1..0f93d9b 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -148,14 +148,19 @@ 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_); + 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 #include #include +#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); @@ -32,10 +33,11 @@ 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); 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); @@ -43,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) { @@ -82,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); @@ -95,67 +101,96 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) for (auto i=0; i> 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++; @@ -263,7 +308,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 { @@ -271,6 +316,39 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator) continue; } } + 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; @@ -305,6 +383,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> 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; @@ -433,6 +545,15 @@ static void synchronize_polynomials(Element& e) { } +static void synchronize_polynomials_kick(Element& e) { + + unsigned int size = (e.polynom_kickx.size() > e.polynom_kicky.size()) ? e.polynom_kickx.size() : e.polynom_kicky.size(); + e.polynom_kickx.resize(size, 0); + e.polynom_kicky.resize(size, 0); + +} + +//This method is only used to read tracy3 files. static void read_polynomials(std::ifstream& fp, Element& e) { unsigned int nr_monomials, n_design, order; e.polynom_a = std::vector(Element::default_polynom); @@ -519,18 +640,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 aace606..8e324f3 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -16,128 +16,258 @@ #include #include -#include +#include #include +#include #include +#include + -std::vector kicktable_list; +std::list Kicktable::kicktable_list; Kicktable::Kicktable(const std::string& filename_) : + filename("") +{ + if (filename_ != "") + this->load_from_file(filename_); +} + +Kicktable::Kicktable( + const std::vector& x_pos, + const std::vector& y_pos, + const std::vector& x_kick, + const std::vector& y_kick, + const double length +) : filename(""), - x_nrpts(0), y_nrpts(0), - x_min(nan("")), x_max(nan("")), - y_min(nan("")), y_max(nan("")) { + x_pos(x_pos), y_pos(y_pos), + x_kick(x_kick), y_kick(y_kick), + length(length) +{} - if (filename_ != "") { - this->load_from_file(filename_); - } +bool Kicktable::is_valid_kicktable() const +{ + if (x_pos.size() <= 1) + return false; + if (y_pos.size() <= 1) + return false; + if (x_kick.size() != x_pos.size() * y_pos.size()) + return false; + if (y_kick.size() != x_kick.size()) + return false; + 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; } -Status::type Kicktable::load_from_file(const std::string& filename_) { - - 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( + const 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_; } - this->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 >> 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); // 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); + 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; - } + 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) { - 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 + 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); // 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[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 && y_pos[1] < y_pos[0]) { + for(unsigned int i=0; i=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& kicktable_list) { +int Kicktable::add_kicktable( + const std::vector& x_pos, + const std::vector& y_pos, + const std::vector& x_kick, + const std::vector& y_kick, + const double length +) +{ + 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_ktab; + new_ktab.load_from_file(filename, true); + return Kicktable::add_kicktable(new_ktab); +} + +int Kicktable::add_kicktable(const Kicktable &new_ktab) +{ + if (not new_ktab.is_valid_kicktable()) + return -1; + + // looks through vector of kicktables... + const auto it = std::find( + kicktable_list.begin(), kicktable_list.end(), new_ktab); + + if (it == kicktable_list.end()) + { + kicktable_list.push_back(new_ktab); + return (int) kicktable_list.size() - 1; + } + + return std::distance(kicktable_list.begin(), it); +} + +void Kicktable::clear_kicktables() { kicktable_list.clear(); } 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_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; } - -template Status::type kicktable_getkicks(const int& kicktable_idx, const double& rx, const double& ry, double& hkick, double& vkick); 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