From 6784973c2101e3a95daeb80682fad70d2439d1d5 Mon Sep 17 00:00:00 2001 From: vitor Date: Mon, 3 Jun 2024 08:33:00 -0300 Subject: [PATCH 01/20] testing-different-cavity-pass-methods --- Makefile | 42 +++++------ include/trackcpp/accelerator.h | 1 + include/trackcpp/auxiliary.h | 45 ++++++++++-- include/trackcpp/passmethods.h | 5 +- include/trackcpp/passmethods.hpp | 121 ++++++++++++++++++++++++++++++- include/trackcpp/tracking.h | 39 +++++++--- python_package/interface.cpp | 7 +- python_package/interface.h | 3 +- src/accelerator.cpp | 31 ++++++++ src/elements.cpp | 2 +- src/flat_file.cpp | 2 +- src/lattice.cpp | 2 +- src/tracking.cpp | 89 ++++++++++++++++++++++- 13 files changed, 338 insertions(+), 51 deletions(-) diff --git a/Makefile b/Makefile index d573057..afbbfa0 100644 --- a/Makefile +++ b/Makefile @@ -85,27 +85,27 @@ PYTHON_PACKAGE_DIR = python_package $(shell touch $(SRCDIR)/output.cpp) # this is so that last compilation time always goes into executable -WARNINGS_CFLAGS += -Wall -WARNINGS_CFLAGS += -Wextra # reasonable and standard -WARNINGS_CFLAGS += -Wshadow # warn the user if a variable declaration shadows one from a parent context -WARNINGS_CFLAGS += -Wnon-virtual-dtor # warn the user if a class with virtual functions has a non-virtual destructor. - # This helps catch hard to track down memory errors -# WARNINGS_CFLAGS += -Wdouble-promotion # warn if float is implicit promoted to double -# WARNINGS_CFLAGS += -Wformat=2 # warn on security issues around functions that format output (ie printf) -# WARNINGS_CFLAGS += -Wold-style-cast # warn for c-style casts -# WARNINGS_CFLAGS += -Wsign-conversion # warn on sign conversions -WARNINGS_CFLAGS += -Wcast-align # warn for potential performance problem casts -WARNINGS_CFLAGS += -Wconversion # warn on type conversions that may lose data -# WARNINGS_CFLAGS += -Wduplicated-branches # warn if if / else branches have duplicated code // not available in g++ 6.3.0 -WARNINGS_CFLAGS += -Wduplicated-cond # warn if if / else chain has duplicated conditions -# WARNINGS_CFLAGS += -Wimplicit-fallthrough # warn on statements that fallthrough without an explicit annotation // not available in g++ 6.3.0 -WARNINGS_CFLAGS += -Wlogical-op # warn about logical operations being used where bitwise were probably wanted -WARNINGS_CFLAGS += -Wmisleading-indentation # warn if indentation implies blocks where blocks do not exist -WARNINGS_CFLAGS += -Wnull-dereference # warn if a null dereference is detected -WARNINGS_CFLAGS += -Woverloaded-virtual # warn if you overload (not override) a virtual function -WARNINGS_CFLAGS += -Wpedantic # warn if non-standard C++ is used -WARNINGS_CFLAGS += -Wunused # warn on anything being unused -WARNINGS_CFLAGS += -Wuseless-cast # warn if you perform a cast to the same type +# WARNINGS_CFLAGS += -Wall +# WARNINGS_CFLAGS += -Wextra # reasonable and standard +# WARNINGS_CFLAGS += -Wshadow # warn the user if a variable declaration shadows one from a parent context +# WARNINGS_CFLAGS += -Wnon-virtual-dtor # warn the user if a class with virtual functions has a non-virtual destructor. +# # This helps catch hard to track down memory errors +# # WARNINGS_CFLAGS += -Wdouble-promotion # warn if float is implicit promoted to double +# # WARNINGS_CFLAGS += -Wformat=2 # warn on security issues around functions that format output (ie printf) +# # WARNINGS_CFLAGS += -Wold-style-cast # warn for c-style casts +# # WARNINGS_CFLAGS += -Wsign-conversion # warn on sign conversions +# WARNINGS_CFLAGS += -Wcast-align # warn for potential performance problem casts +# WARNINGS_CFLAGS += -Wconversion # warn on type conversions that may lose data +# # WARNINGS_CFLAGS += -Wduplicated-branches # warn if if / else branches have duplicated code // not available in g++ 6.3.0 +# WARNINGS_CFLAGS += -Wduplicated-cond # warn if if / else chain has duplicated conditions +# # WARNINGS_CFLAGS += -Wimplicit-fallthrough # warn on statements that fallthrough without an explicit annotation // not available in g++ 6.3.0 +# WARNINGS_CFLAGS += -Wlogical-op # warn about logical operations being used where bitwise were probably wanted +# WARNINGS_CFLAGS += -Wmisleading-indentation # warn if indentation implies blocks where blocks do not exist +# WARNINGS_CFLAGS += -Wnull-dereference # warn if a null dereference is detected +# WARNINGS_CFLAGS += -Woverloaded-virtual # warn if you overload (not override) a virtual function +# WARNINGS_CFLAGS += -Wpedantic # warn if non-standard C++ is used +# WARNINGS_CFLAGS += -Wunused # warn on anything being unused +# WARNINGS_CFLAGS += -Wuseless-cast # warn if you perform a cast to the same type ifeq ($(MAKECMDGOALS),trackcpp-debug) CFLAGS = $(MACHINE) $(DBG_FLAG) $(DFLAGS) -pthread diff --git a/include/trackcpp/accelerator.h b/include/trackcpp/accelerator.h index 0653507..98c424f 100644 --- a/include/trackcpp/accelerator.h +++ b/include/trackcpp/accelerator.h @@ -39,6 +39,7 @@ class Accelerator { bool operator!=(const Accelerator& o) const { return !(*this == o); }; bool isequal(const Accelerator& a) const { return *this == a; } // necessary for python_package double get_length() const; + double get_time_aware_frac() const; friend std::ostream& operator<< (std::ostream &out, const Accelerator& a); }; diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 65cfca5..d92679e 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -31,28 +31,55 @@ class PassMethodsClass { static const int pm_str_mpole_symplectic4_pass = 2; static const int pm_bnd_mpole_symplectic4_pass = 3; static const int pm_corrector_pass = 4; - static const int pm_cavity_pass = 5; + static const int pm_cavity_0_pass = 5; static const int pm_thinquad_pass = 6; static const int pm_thinsext_pass = 7; 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_cavity_1comp_pass = 11; + static const int pm_cavity_1frac_pass = 12; + static const int pm_cavity_2_pass = 13; + static const int pm_nr_pms = 14; // counter for number of passmethods PassMethodsClass() { passmethods.push_back("identity_pass"); passmethods.push_back("drift_pass"); passmethods.push_back("str_mpole_symplectic4_pass"); passmethods.push_back("bnd_mpole_symplectic4_pass"); passmethods.push_back("corrector_pass"); - passmethods.push_back("cavity_pass"); + passmethods.push_back("cavity_0_pass"); passmethods.push_back("thinquad_pass"); passmethods.push_back("thinsext_pass"); passmethods.push_back("kicktable_pass"); passmethods.push_back("matrix_pass"); passmethods.push_back("drift_g2l_pass"); + passmethods.push_back("cavity_1comp_pass"); + passmethods.push_back("cavity_1frac_pass"); + passmethods.push_back("cavity_2_pass"); } int size() const { return passmethods.size(); } std::string operator[](const int i) const { return passmethods[i]; } + bool is_time_aware_pm(const int i) const { + bool flag = false; + switch (i) + { + case pm_cavity_0_pass: + flag = true; + break; + case pm_cavity_1comp_pass: + flag = true; + break; + case pm_cavity_1frac_pass: + flag = true; + break; + case pm_cavity_2_pass: + flag = true; + break; + default: + break; + } + return flag; + }; private: std::vector passmethods; }; @@ -65,13 +92,16 @@ struct PassMethod { pm_str_mpole_symplectic4_pass = 2, pm_bnd_mpole_symplectic4_pass = 3, pm_corrector_pass = 4, - pm_cavity_pass = 5, + pm_cavity_0_pass = 5, pm_thinquad_pass = 6, pm_thinsext_pass = 7, pm_kickmap_pass = 8, pm_matrix_pass = 9, pm_drift_g2l_pass = 10, - pm_nr_pms = 11, + pm_cavity_1comp_pass = 11, + pm_cavity_1frac_pass = 12, + pm_cavity_2_pass = 13, + pm_nr_pms = 14, }; }; @@ -83,12 +113,15 @@ const std::vector pm_dict = { "str_mpole_symplectic4_pass", "bnd_mpole_symplectic4_pass", "corrector_pass", - "cavity_pass", + "cavity_0_pass", "thinquad_pass", "thinsext_pass", "kicktable_pass", "matrix_pass", "drift_g2l_pass", + "cavity_1comp_pass", + "cavity_1frac_pass", + "cavity_2_pass", }; struct RadiationState { diff --git a/include/trackcpp/passmethods.h b/include/trackcpp/passmethods.h index 51a38aa..1ad9bbe 100644 --- a/include/trackcpp/passmethods.h +++ b/include/trackcpp/passmethods.h @@ -60,7 +60,10 @@ template Status::type pm_drift_pass (Pos &pos, c template Status::type pm_str_mpole_symplectic4_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_bnd_mpole_symplectic4_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_corrector_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); -template Status::type pm_cavity_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_cavity_0_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_cavity_1comp_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_cavity_1frac_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_cavity_2_pass (Pos &pos, const Element &elem, const Accelerator& accelerator, const double nturn = 0.0); template Status::type pm_thinquad_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_thinsext_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_kickmap_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index a42e10e..b918c55 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -387,7 +387,7 @@ Status::type pm_corrector_pass(Pos &pos, const Element &elem, } template -Status::type pm_cavity_pass(Pos &pos, const Element &elem, +Status::type pm_cavity_0_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); @@ -420,6 +420,125 @@ Status::type pm_cavity_pass(Pos &pos, const Element &elem, return Status::success; } +template +Status::type pm_cavity_1comp_pass(Pos &pos, const Element &elem, + const Accelerator& accelerator) { + + if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); + + global_2_local(pos, elem); + double nv = elem.voltage / accelerator.energy; + double frf = elem.frequency; + double harmonic_number = accelerator.harmonic_number; + double L0 = accelerator.get_length(); + double ddl = (light_speed*harmonic_number/frf - L0); + if (elem.length == 0) { + T &de = pos.de, &dl = pos.dl; + de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); + dl -= ddl; + } else { + T &rx = pos.rx, &px = pos.px; + T &ry = pos.ry, &py = pos.py; + T &de = pos.de, &dl = pos.dl; + // drift half length + T pnorm = 1 / (1 + de); + T norml = (0.5 * elem.length) * pnorm; + rx += norml * px; + ry += norml * py; + dl += 0.5 * norml * pnorm * (px*px + py*py); + // longitudinal momentum kick + de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); + dl -= ddl; + // drift half length + pnorm = 1.0 / (1.0 + de); + norml = (0.5 * elem.length) * pnorm; + rx += norml * px; + ry += norml * py; + dl += 0.5 * norml * pnorm * (px*px + py*py); + } + local_2_global(pos, elem); + return Status::success; +} + +template +Status::type pm_cavity_1frac_pass(Pos &pos, const Element &elem, + const Accelerator& accelerator) { + + if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); + + global_2_local(pos, elem); + double nv = elem.voltage / accelerator.energy; + double frf = elem.frequency; + double harmonic_number = accelerator.harmonic_number; + double L0 = accelerator.get_length(); + double time_aware_frac = accelerator.get_time_aware_frac(); + double ddl = (light_speed*harmonic_number/frf - L0)*time_aware_frac; + if (elem.length == 0) { + T &de = pos.de, &dl = pos.dl; + de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); + dl -= ddl; + } else { + T &rx = pos.rx, &px = pos.px; + T &ry = pos.ry, &py = pos.py; + T &de = pos.de, &dl = pos.dl; + // drift half length + T pnorm = 1 / (1 + de); + T norml = (0.5 * elem.length) * pnorm; + rx += norml * px; + ry += norml * py; + dl += 0.5 * norml * pnorm * (px*px + py*py); + // longitudinal momentum kick + de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); + dl -= ddl; + // drift half length + pnorm = 1.0 / (1.0 + de); + norml = (0.5 * elem.length) * pnorm; + rx += norml * px; + ry += norml * py; + dl += 0.5 * norml * pnorm * (px*px + py*py); + } + local_2_global(pos, elem); + return Status::success; +} + +template +Status::type pm_cavity_2_pass(Pos &pos, const Element &elem, + const Accelerator& accelerator, const double nturn) { + + if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); + + global_2_local(pos, elem); + double nv = elem.voltage / accelerator.energy; + double frf = elem.frequency; + double harmonic_number = accelerator.harmonic_number; + double L0 = accelerator.get_length(); + double ddl = nturn * (light_speed*harmonic_number/frf - L0); + if (elem.length == 0) { + T &de = pos.de, &dl = pos.dl; + de += -nv * sin(TWOPI*frf *(dl-ddl)/ light_speed - elem.phase_lag); + } else { + T &rx = pos.rx, &px = pos.px; + T &ry = pos.ry, &py = pos.py; + T &de = pos.de, &dl = pos.dl; + // drift half length + T pnorm = 1 / (1 + de); + T norml = (0.5 * elem.length) * pnorm; + rx += norml * px; + ry += norml * py; + dl += 0.5 * norml * pnorm * (px*px + py*py); + // longitudinal momentum kick + de += -nv * sin(TWOPI*frf *(dl-ddl)/ light_speed - elem.phase_lag); + // drift half length + pnorm = 1.0 / (1.0 + de); + norml = (0.5 * elem.length) * pnorm; + rx += norml * px; + ry += norml * py; + dl += 0.5 * norml * pnorm * (px*px + py*py); + } + local_2_global(pos, elem); + return Status::success; +} + template Status::type pm_thinquad_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index 1e48205..3f36dd1 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -44,6 +44,10 @@ Status::type track_findorbit6( Accelerator& accelerator, std::vector >& closed_orbit, const Pos& fixed_point_guess = Pos(0)); +Status::type track_findorbit6_dct( + Accelerator& accelerator, std::vector >& closed_orbit, + const Pos& fixed_point_guess = Pos(0), const double dct = 0.0); + Pos linalg_solve4_posvec( const std::vector >& M, const Pos& b); @@ -55,7 +59,7 @@ template Status::type track_elementpass ( const Element& el, // element through which to track particle Pos &orig_pos, // initial electron coordinates - const Accelerator& accelerator) { + const Accelerator& accelerator, const double nturn = 0.0) { Status::type status = Status::success; @@ -75,8 +79,8 @@ Status::type track_elementpass ( case PassMethod::pm_corrector_pass: if ((status = pm_corrector_pass(orig_pos, el, accelerator)) != Status::success) return status; break; - case PassMethod::pm_cavity_pass: - if ((status = pm_cavity_pass(orig_pos, el, accelerator)) != Status::success) return status; + case PassMethod::pm_cavity_0_pass: + if ((status = pm_cavity_0_pass(orig_pos, el, accelerator)) != Status::success) return status; break; case PassMethod::pm_thinquad_pass: if ((status = pm_thinquad_pass(orig_pos, el, accelerator)) != Status::success) return status; @@ -93,6 +97,15 @@ 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_cavity_1comp_pass: + if ((status = pm_cavity_1comp_pass(orig_pos, el, accelerator)) != Status::success) return status; + break; + case PassMethod::pm_cavity_1frac_pass: + if ((status = pm_cavity_1frac_pass(orig_pos, el, accelerator)) != Status::success) return status; + break; + case PassMethod::pm_cavity_2_pass: + if ((status = pm_cavity_2_pass(orig_pos, el, accelerator, nturn)) != Status::success) return status; + break; default: return Status::passmethod_not_defined; } @@ -105,12 +118,12 @@ template Status::type track_elementpass ( const Element& el, // element through which to track particle std::vector >& orig_pos, // initial electron coordinates - const Accelerator& accelerator) { + const Accelerator& accelerator, const double nturn) { Status::type status = Status::success; for(auto&& pos: orig_pos) { - Status::type status2 = track_elementpass(el, pos, accelerator); + Status::type status2 = track_elementpass(el, pos, accelerator, nturn); if (status2 != Status::success) status = status2; } return status; @@ -140,7 +153,7 @@ Status::type track_linepass ( std::vector >& pos, // vector with tracked electron coordinates at start of every element and at the end of last one. unsigned int& element_offset, // index of starting element for tracking Plane::type& lost_plane, // return plane in which particle was lost, if the case. - std::vector& indices) {// indices to return; + std::vector& indices, double nturn = 0.0) {// indices to return; Status::type status = Status::success; lost_plane = Plane::no_plane; @@ -165,7 +178,7 @@ Status::type track_linepass ( // stores trajectory at entrance of each element if (indcs[i]) pos.push_back(orig_pos); - status = track_elementpass (element, orig_pos, accelerator); + status = track_elementpass (element, orig_pos, accelerator, nturn); const T& rx = orig_pos.rx; const T& ry = orig_pos.ry; @@ -256,7 +269,7 @@ Status::type track_linepass ( std::vector >& pos, // vector with tracked electron coordinates at start of every element and at the end of last one. unsigned int& element_offset, // index of starting element for tracking Plane::type& lost_plane, // return plane in which particle was lost, if the case. - bool trajectory) { // whether function should return coordinates at all elements + bool trajectory, double nturn = 0.0) { // whether function should return coordinates at all elements std::vector indices; unsigned int nr_elements = accelerator.lattice.size(); @@ -269,7 +282,7 @@ Status::type track_linepass ( } return track_linepass ( - accelerator, orig_pos, pos, element_offset, lost_plane, indices); + accelerator, orig_pos, pos, element_offset, lost_plane, indices, nturn); } @@ -281,7 +294,7 @@ Status::type track_linepass ( unsigned int element_offset, std::vector& lost_plane, std::vector& lost_element, - std::vector& indices) { + std::vector& indices, double nturn = 0.0) { int nr_elements = accelerator.lattice.size(); Status::type status = Status::success; @@ -297,7 +310,7 @@ Status::type track_linepass ( unsigned int le = element_offset; status2 = track_linepass ( - accelerator, orig_pos[i], final_pos, le, lp, indices); + accelerator, orig_pos[i], final_pos, le, lp, indices, nturn); if (status2 != Status::success) status = status2; @@ -340,6 +353,7 @@ Status::type track_ringpass ( Status::type status = Status::success; std::vector > final_pos; + double nturn = 0.0; if (trajectory) pos.reserve(nr_turns+1); @@ -348,7 +362,7 @@ Status::type track_ringpass ( // stores trajectory at beggining of each turn if (trajectory) pos.push_back(orig_pos); - if ((status = track_linepass (accelerator, orig_pos, final_pos, element_offset, lost_plane, false)) != Status::success) { + if ((status = track_linepass (accelerator, orig_pos, final_pos, element_offset, lost_plane, false, nturn)) != Status::success) { // fill last of vector with nans pos.emplace_back( @@ -360,6 +374,7 @@ Status::type track_ringpass ( return status; } final_pos.clear(); + nturn += 1.0; } pos.push_back(orig_pos); diff --git a/python_package/interface.cpp b/python_package/interface.cpp index bcd111f..a74ad12 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -21,7 +21,7 @@ Status::type track_elementpass_wrapper ( const Element& el, double *pos, int n1, int n2, - const Accelerator& accelerator) { + const Accelerator& accelerator, const double nturn) { std::vector > post; @@ -33,7 +33,7 @@ Status::type track_elementpass_wrapper ( pos[4*n2 + i], pos[5*n2 + i]); } Status::type status = track_elementpass( - el, post, accelerator); + el, post, accelerator, nturn); for (unsigned int i=0; i indices; std::vector< unsigned int > lost_plane; std::vector< unsigned int > lost_element; + double nturn = 0.0; }; struct RingPassArgs : public LinePassArgs { @@ -52,7 +53,7 @@ struct String { Status::type track_elementpass_wrapper ( const Element& el, double *pos, int n1, int n2, - const Accelerator& accelerator); + const Accelerator& accelerator, const double nturn = 0.0); Status::type track_linepass_wrapper ( const Accelerator& accelerator, diff --git a/src/accelerator.cpp b/src/accelerator.cpp index 0e83afe..1e5ab9a 100644 --- a/src/accelerator.cpp +++ b/src/accelerator.cpp @@ -27,6 +27,37 @@ double Accelerator::get_length() const { return length; } +inline bool check_time_aware_pm(const int i){ + bool flag = false; + switch (i) + { + case PassMethod::pm_cavity_0_pass: + flag = true; + break; + case PassMethod::pm_cavity_1comp_pass: + flag = true; + break; + case PassMethod::pm_cavity_1frac_pass: + flag = true; + break; + case PassMethod::pm_cavity_2_pass: + flag = true; + break; + default: + break; + } + return flag; +} + +double Accelerator::get_time_aware_frac() const { + double count = 0.0; + for(auto i=0; ienergy != o.energy) return false; diff --git a/src/elements.cpp b/src/elements.cpp index f33eb4c..171b431 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -278,7 +278,7 @@ 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) { - element.pass_method = PassMethod::pm_cavity_pass; + element.pass_method = PassMethod::pm_cavity_0_pass; element.frequency = frequency; element.voltage = voltage; element.phase_lag = phase_lag; diff --git a/src/flat_file.cpp b/src/flat_file.cpp index 52ff4f3..cebb970 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -365,7 +365,7 @@ static Status::type read_flat_file_tracy(const std::string& filename, Accelerato }; break; case FlatFileType::cavity: { - e.pass_method = PassMethod::pm_cavity_pass; + e.pass_method = PassMethod::pm_cavity_0_pass; int hnumber; double energy; fp >> e.voltage >> e.frequency >> hnumber >> energy; e.voltage *= energy; e.frequency *= light_speed / (2*M_PI); diff --git a/src/lattice.cpp b/src/lattice.cpp index 910a841..bfec733 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -78,7 +78,7 @@ void latt_setcavity(std::vector& lattice, const std::string& state) { for(unsigned int i=0; i& lattice = accelerator.lattice; Pos fp = fixed_point; - + const int radsts = accelerator.radiation_on; if (radsts == RadiationState::full){ accelerator.radiation_on = RadiationState::damping; @@ -107,9 +107,9 @@ Status::type track_findm66 (Accelerator& accelerator, // constant term of the final map v0.rx = map.rx.c[0]; v0.px = map.px.c[0]; v0.ry = map.ry.c[0]; v0.py = map.py.c[0]; v0.de = map.de.c[0]; v0.dl = map.dl.c[0]; - + accelerator.radiation_on = radsts; - + return status; } @@ -220,6 +220,89 @@ Status::type track_findorbit6( } +Status::type track_findorbit6_dct( + Accelerator& accelerator, + std::vector >& closed_orbit, + const Pos& fixed_point_guess, const double dct) { + + const std::vector& the_ring = accelerator.lattice; + + double delta = 1e-9; // [m],[rad],[dE/E] + double tolerance = 2.22044604925e-14; + int max_nr_iters = 50; + + const int radsts = accelerator.radiation_on; + if (radsts == RadiationState::full){ + accelerator.radiation_on = RadiationState::damping; + } + // calcs longitudinal fixed point + + // temporary vectors and matrices + std::vector > co(7,0); + for(auto i=0; i<7; ++i) co[i] = fixed_point_guess; + std::vector > co2(7,0); + std::vector > D(7,0); + std::vector > M(6,0); + Pos dco(1.0,1.0,1.0,1.0,1.0,1.0); + Pos theta(0.0,0.0,0.0,0.0,0.0,dct); + matrix6_set_identity_posvec(D, delta); + + int nr_iter = 0; + while ((get_max(dco) > tolerance) and (nr_iter <= max_nr_iters)) { + co = co + D; + Pos Ri = co[6]; + std::vector > co2; + unsigned int element_offset = 0; + Plane::type lost_plane; + Status::type status = Status::success; + + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[0], co2, element_offset, lost_plane, false)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[1], co2, element_offset, lost_plane, false)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[2], co2, element_offset, lost_plane, false)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[3], co2, element_offset, lost_plane, false)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[4], co2, element_offset, lost_plane, false)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[5], co2, element_offset, lost_plane, false)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[6], co2, element_offset, lost_plane, false)); + + if (status != Status::success) { + return Status::findorbit_one_turn_matrix_problem; + } + + Pos Rf = co2[6]; + + M[0] = (co2[0] - Rf) / delta; + M[1] = (co2[1] - Rf) / delta; + M[2] = (co2[2] - Rf) / delta; + M[3] = (co2[3] - Rf) / delta; + M[4] = (co2[4] - Rf) / delta; + M[5] = (co2[5] - Rf) / delta; + + Pos b = Rf - Ri - theta; + std::vector > M_1(6,0); + matrix6_set_identity_posvec(M_1); + M_1 = M_1 - M; + dco = linalg_solve6_posvec(M_1, b); + co[6] = dco + Ri; + co[0] = co[6]; co[1] = co[6]; + co[2] = co[6]; co[3] = co[6]; + co[4] = co[6]; co[5] = co[6]; + nr_iter++; + } + + if (nr_iter > max_nr_iters) { + return Status::findorbit_not_converged; + } + + // propagates fixed point throught the_ring + closed_orbit.clear(); + unsigned int element_offset = 0; + Plane::type lost_plane; + track_linepass(accelerator, co[6], closed_orbit, element_offset, lost_plane, true); + accelerator.radiation_on = radsts; + return Status::success; + +} + Status::type track_findorbit4( Accelerator& accelerator, std::vector >& closed_orbit, From 6129eddbc4af1356db2558a9f510f140eda6ff09 Mon Sep 17 00:00:00 2001 From: vitor Date: Mon, 3 Jun 2024 10:54:17 -0300 Subject: [PATCH 02/20] restore comments in makefile --- Makefile | 42 +++++++++++++++++++++--------------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/Makefile b/Makefile index afbbfa0..d573057 100644 --- a/Makefile +++ b/Makefile @@ -85,27 +85,27 @@ PYTHON_PACKAGE_DIR = python_package $(shell touch $(SRCDIR)/output.cpp) # this is so that last compilation time always goes into executable -# WARNINGS_CFLAGS += -Wall -# WARNINGS_CFLAGS += -Wextra # reasonable and standard -# WARNINGS_CFLAGS += -Wshadow # warn the user if a variable declaration shadows one from a parent context -# WARNINGS_CFLAGS += -Wnon-virtual-dtor # warn the user if a class with virtual functions has a non-virtual destructor. -# # This helps catch hard to track down memory errors -# # WARNINGS_CFLAGS += -Wdouble-promotion # warn if float is implicit promoted to double -# # WARNINGS_CFLAGS += -Wformat=2 # warn on security issues around functions that format output (ie printf) -# # WARNINGS_CFLAGS += -Wold-style-cast # warn for c-style casts -# # WARNINGS_CFLAGS += -Wsign-conversion # warn on sign conversions -# WARNINGS_CFLAGS += -Wcast-align # warn for potential performance problem casts -# WARNINGS_CFLAGS += -Wconversion # warn on type conversions that may lose data -# # WARNINGS_CFLAGS += -Wduplicated-branches # warn if if / else branches have duplicated code // not available in g++ 6.3.0 -# WARNINGS_CFLAGS += -Wduplicated-cond # warn if if / else chain has duplicated conditions -# # WARNINGS_CFLAGS += -Wimplicit-fallthrough # warn on statements that fallthrough without an explicit annotation // not available in g++ 6.3.0 -# WARNINGS_CFLAGS += -Wlogical-op # warn about logical operations being used where bitwise were probably wanted -# WARNINGS_CFLAGS += -Wmisleading-indentation # warn if indentation implies blocks where blocks do not exist -# WARNINGS_CFLAGS += -Wnull-dereference # warn if a null dereference is detected -# WARNINGS_CFLAGS += -Woverloaded-virtual # warn if you overload (not override) a virtual function -# WARNINGS_CFLAGS += -Wpedantic # warn if non-standard C++ is used -# WARNINGS_CFLAGS += -Wunused # warn on anything being unused -# WARNINGS_CFLAGS += -Wuseless-cast # warn if you perform a cast to the same type +WARNINGS_CFLAGS += -Wall +WARNINGS_CFLAGS += -Wextra # reasonable and standard +WARNINGS_CFLAGS += -Wshadow # warn the user if a variable declaration shadows one from a parent context +WARNINGS_CFLAGS += -Wnon-virtual-dtor # warn the user if a class with virtual functions has a non-virtual destructor. + # This helps catch hard to track down memory errors +# WARNINGS_CFLAGS += -Wdouble-promotion # warn if float is implicit promoted to double +# WARNINGS_CFLAGS += -Wformat=2 # warn on security issues around functions that format output (ie printf) +# WARNINGS_CFLAGS += -Wold-style-cast # warn for c-style casts +# WARNINGS_CFLAGS += -Wsign-conversion # warn on sign conversions +WARNINGS_CFLAGS += -Wcast-align # warn for potential performance problem casts +WARNINGS_CFLAGS += -Wconversion # warn on type conversions that may lose data +# WARNINGS_CFLAGS += -Wduplicated-branches # warn if if / else branches have duplicated code // not available in g++ 6.3.0 +WARNINGS_CFLAGS += -Wduplicated-cond # warn if if / else chain has duplicated conditions +# WARNINGS_CFLAGS += -Wimplicit-fallthrough # warn on statements that fallthrough without an explicit annotation // not available in g++ 6.3.0 +WARNINGS_CFLAGS += -Wlogical-op # warn about logical operations being used where bitwise were probably wanted +WARNINGS_CFLAGS += -Wmisleading-indentation # warn if indentation implies blocks where blocks do not exist +WARNINGS_CFLAGS += -Wnull-dereference # warn if a null dereference is detected +WARNINGS_CFLAGS += -Woverloaded-virtual # warn if you overload (not override) a virtual function +WARNINGS_CFLAGS += -Wpedantic # warn if non-standard C++ is used +WARNINGS_CFLAGS += -Wunused # warn on anything being unused +WARNINGS_CFLAGS += -Wuseless-cast # warn if you perform a cast to the same type ifeq ($(MAKECMDGOALS),trackcpp-debug) CFLAGS = $(MACHINE) $(DBG_FLAG) $(DFLAGS) -pthread From b1860924d328affbdbbf0929a36d53db2f1a0a79 Mon Sep 17 00:00:00 2001 From: vitor Date: Mon, 3 Jun 2024 10:54:31 -0300 Subject: [PATCH 03/20] adapt AT cavitypass mod --- include/trackcpp/passmethods.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index f151033..552e33e 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -518,7 +518,8 @@ Status::type pm_cavity_2_pass(Pos &pos, const Element &elem, double frf = elem.frequency; double harmonic_number = accelerator.harmonic_number; double L0 = accelerator.get_length(); - double ddl = nturn * (light_speed*harmonic_number/frf - L0); + double time_aware_frac = accelerator.get_time_aware_frac(); + double ddl = (light_speed*harmonic_number/frf - L0) * time_aware_frac * nturn; if (elem.length == 0) { T &de = pos.de, &dl = pos.dl; de += -nv * sin(TWOPI*frf *(dl-ddl)/ light_speed - elem.phase_lag); From b6dd5754b16341820e0767b6846a17bccd5be499 Mon Sep 17 00:00:00 2001 From: vitor Date: Thu, 6 Jun 2024 12:26:16 -0300 Subject: [PATCH 04/20] adjust in cav pass 1 comp --- include/trackcpp/passmethods.hpp | 85 +++++++++++++++++++++++++++++--- 1 file changed, 79 insertions(+), 6 deletions(-) diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 552e33e..bfbac65 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -426,6 +426,66 @@ Status::type pm_cavity_0_pass(Pos &pos, const Element &elem, return Status::success; } +// template +// Status::type pm_cavity_1comp_pass(Pos &pos, const Element &elem, +// const Accelerator& accelerator) { + +// if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); + +// global_2_local(pos, elem); +// double nv = elem.voltage / accelerator.energy; +// double frf = elem.frequency; +// double harmonic_number = accelerator.harmonic_number; +// double L0 = accelerator.get_length(); +// double s = 0.0; +// double accum_s = 0.0; +// unsigned int last_cav_idx = 0; +// for (unsigned int i = 0; i < accelerator.lattice.size(); i++){ +// auto elem2 = accelerator.lattice[i]; +// if (elem2 == elem){ +// accum_s += s; +// } +// if (elem2.frequency != 0.0){ +// s = 0.0; +// last_cav_idx = i; +// } +// s += elem2.length; +// } +// double ddl = (light_speed*harmonic_number/frf - L0); +// if (elem.length == 0) { +// T &de = pos.de, &dl = pos.dl; +// dl -= ddl*(accum_s/L0); +// de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); +// if (elem == accelerator.lattice[last_cav_idx]){ +// dl -= ddl*(s/L0); +// } +// } else { +// T &rx = pos.rx, &px = pos.px; +// T &ry = pos.ry, &py = pos.py; +// T &de = pos.de, &dl = pos.dl; +// // drift half length +// T pnorm = 1 / (1 + de); +// T norml = (0.5 * elem.length) * pnorm; +// rx += norml * px; +// ry += norml * py; +// dl += 0.5 * norml * pnorm * (px*px + py*py); +// // longitudinal momentum kick +// dl -= ddl*(accum_s/L0); +// de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); +// if (elem == accelerator.lattice[last_cav_idx]){ +// dl -= ddl*(s/L0); +// } +// // drift half length +// pnorm = 1.0 / (1.0 + de); +// norml = (0.5 * elem.length) * pnorm; +// rx += norml * px; +// ry += norml * py; +// dl += 0.5 * norml * pnorm * (px*px + py*py); +// } +// local_2_global(pos, elem); +// return Status::success; +// } + template Status::type pm_cavity_1comp_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { @@ -437,11 +497,21 @@ Status::type pm_cavity_1comp_pass(Pos &pos, const Element &elem, double frf = elem.frequency; double harmonic_number = accelerator.harmonic_number; double L0 = accelerator.get_length(); + unsigned int last_cav_idx = 0; + for (unsigned int i = 0; i < accelerator.lattice.size(); i++){ + auto elem2 = accelerator.lattice[i]; + if (elem2.frequency != 0.0){ + last_cav_idx = i; + } + } double ddl = (light_speed*harmonic_number/frf - L0); if (elem.length == 0) { T &de = pos.de, &dl = pos.dl; + de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); - dl -= ddl; + if (elem == accelerator.lattice[last_cav_idx]){ + dl -= ddl; + } } else { T &rx = pos.rx, &px = pos.px; T &ry = pos.ry, &py = pos.py; @@ -453,8 +523,11 @@ Status::type pm_cavity_1comp_pass(Pos &pos, const Element &elem, ry += norml * py; dl += 0.5 * norml * pnorm * (px*px + py*py); // longitudinal momentum kick + de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); - dl -= ddl; + if (elem == accelerator.lattice[last_cav_idx]){ + dl -= ddl; + } // drift half length pnorm = 1.0 / (1.0 + de); norml = (0.5 * elem.length) * pnorm; @@ -481,8 +554,8 @@ Status::type pm_cavity_1frac_pass(Pos &pos, const Element &elem, double ddl = (light_speed*harmonic_number/frf - L0)*time_aware_frac; if (elem.length == 0) { T &de = pos.de, &dl = pos.dl; - de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); dl -= ddl; + de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); } else { T &rx = pos.rx, &px = pos.px; T &ry = pos.ry, &py = pos.py; @@ -494,8 +567,8 @@ Status::type pm_cavity_1frac_pass(Pos &pos, const Element &elem, ry += norml * py; dl += 0.5 * norml * pnorm * (px*px + py*py); // longitudinal momentum kick - de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); dl -= ddl; + de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); // drift half length pnorm = 1.0 / (1.0 + de); norml = (0.5 * elem.length) * pnorm; @@ -518,8 +591,8 @@ Status::type pm_cavity_2_pass(Pos &pos, const Element &elem, double frf = elem.frequency; double harmonic_number = accelerator.harmonic_number; double L0 = accelerator.get_length(); - double time_aware_frac = accelerator.get_time_aware_frac(); - double ddl = (light_speed*harmonic_number/frf - L0) * time_aware_frac * nturn; + // double time_aware_frac = accelerator.get_time_aware_frac(); + double ddl = (light_speed*harmonic_number/frf - L0)* nturn;// * time_aware_frac ; if (elem.length == 0) { T &de = pos.de, &dl = pos.dl; de += -nv * sin(TWOPI*frf *(dl-ddl)/ light_speed - elem.phase_lag); From 7e9fec6ca3d7a567c782ba308663f3bcb54ba5d2 Mon Sep 17 00:00:00 2001 From: vitor Date: Fri, 7 Jun 2024 14:28:20 -0300 Subject: [PATCH 05/20] test dl cav0pass --- include/trackcpp/tracking.h | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index 3f36dd1..c9df57f 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -157,6 +157,8 @@ Status::type track_linepass ( Status::type status = Status::success; lost_plane = Plane::no_plane; + double timectc = 0; + double tawfrac = accelerator.get_time_aware_frac(); const std::vector& line = accelerator.lattice; int nr_elements = line.size(); @@ -174,12 +176,20 @@ Status::type track_linepass ( // syntactic-sugar for read-only access to element object parameters const Element& element = line[element_offset]; + if (element.frequency != 0.0) { + timectc += tawfrac; + } // stores trajectory at entrance of each element if (indcs[i]) pos.push_back(orig_pos); status = track_elementpass (element, orig_pos, accelerator, nturn); + if ((timectc >= 0.98) && (element.pass_method == PassMethod::pm_cavity_0_pass)) { + orig_pos.dl -= (light_speed*accelerator.harmonic_number/element.frequency - accelerator.get_length()); + timectc = 0.0; + } + const T& rx = orig_pos.rx; const T& ry = orig_pos.ry; From 86ca3e8cab05a171b52942f50c8e1bcf45e8d334 Mon Sep 17 00:00:00 2001 From: vitor Date: Mon, 10 Jun 2024 13:15:16 -0300 Subject: [PATCH 06/20] up --- include/trackcpp/auxiliary.h | 44 +++++--- include/trackcpp/passmethods.h | 8 +- include/trackcpp/passmethods.hpp | 175 ++++++++++++++++++------------- include/trackcpp/tracking.h | 21 ++-- src/accelerator.cpp | 12 ++- 5 files changed, 162 insertions(+), 98 deletions(-) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 597dab3..9fa6d9e 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -37,10 +37,12 @@ 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_cavity_1comp_pass = 11; - static const int pm_cavity_1frac_pass = 12; - static const int pm_cavity_2_pass = 13; - static const int pm_nr_pms = 14; // counter for number of passmethods + static const int pm_cavity_1_pass = 11; + static const int pm_cavity_2_pass = 12; + static const int pm_cavity_3_pass = 13; + static const int pm_cavity_4_pass = 14; + static const int pm_cavity_5_pass = 15; + static const int pm_nr_pms = 16; // counter for number of passmethods PassMethodsClass() { passmethods.push_back("identity_pass"); passmethods.push_back("drift_pass"); @@ -53,9 +55,11 @@ class PassMethodsClass { passmethods.push_back("kicktable_pass"); passmethods.push_back("matrix_pass"); passmethods.push_back("drift_g2l_pass"); - passmethods.push_back("cavity_1comp_pass"); - passmethods.push_back("cavity_1frac_pass"); + passmethods.push_back("cavity_1_pass"); passmethods.push_back("cavity_2_pass"); + passmethods.push_back("cavity_3_pass"); + passmethods.push_back("cavity_4_pass"); + passmethods.push_back("cavity_5_pass"); } int size() const { return passmethods.size(); } std::string operator[](const int i) const { return passmethods[i]; } @@ -66,13 +70,19 @@ class PassMethodsClass { case pm_cavity_0_pass: flag = true; break; - case pm_cavity_1comp_pass: + case pm_cavity_1_pass: flag = true; break; - case pm_cavity_1frac_pass: + case pm_cavity_2_pass: flag = true; break; - case pm_cavity_2_pass: + case pm_cavity_3_pass: + flag = true; + break; + case pm_cavity_4_pass: + flag = true; + break; + case pm_cavity_5_pass: flag = true; break; default: @@ -98,10 +108,12 @@ struct PassMethod { pm_kickmap_pass = 8, pm_matrix_pass = 9, pm_drift_g2l_pass = 10, - pm_cavity_1comp_pass = 11, - pm_cavity_1frac_pass = 12, - pm_cavity_2_pass = 13, - pm_nr_pms = 14, + pm_cavity_1_pass = 11, + pm_cavity_2_pass = 12, + pm_cavity_3_pass = 13, + pm_cavity_4_pass = 14, + pm_cavity_5_pass = 15, + pm_nr_pms = 16, }; }; @@ -119,9 +131,11 @@ const std::vector pm_dict = { "kicktable_pass", "matrix_pass", "drift_g2l_pass", - "cavity_1comp_pass", - "cavity_1frac_pass", + "cavity_1_pass", "cavity_2_pass", + "cavity_3_pass", + "cavity_4_pass", + "cavity_5_pass", }; struct RadiationState { diff --git a/include/trackcpp/passmethods.h b/include/trackcpp/passmethods.h index 1ad9bbe..6aa0a3a 100644 --- a/include/trackcpp/passmethods.h +++ b/include/trackcpp/passmethods.h @@ -61,9 +61,11 @@ template Status::type pm_str_mpole_symplectic4_pass (Pos &pos, c template Status::type pm_bnd_mpole_symplectic4_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_corrector_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_cavity_0_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); -template Status::type pm_cavity_1comp_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); -template Status::type pm_cavity_1frac_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); -template Status::type pm_cavity_2_pass (Pos &pos, const Element &elem, const Accelerator& accelerator, const double nturn = 0.0); +template Status::type pm_cavity_1_pass (Pos &pos, const Element &elem, const Accelerator& accelerator, const double nturn = 0.0); +template Status::type pm_cavity_2_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_cavity_3_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_cavity_4_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_cavity_5_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_thinquad_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_thinsext_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_kickmap_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index bfbac65..056bf7e 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -426,68 +426,47 @@ Status::type pm_cavity_0_pass(Pos &pos, const Element &elem, return Status::success; } -// template -// Status::type pm_cavity_1comp_pass(Pos &pos, const Element &elem, -// const Accelerator& accelerator) { - -// if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); - -// global_2_local(pos, elem); -// double nv = elem.voltage / accelerator.energy; -// double frf = elem.frequency; -// double harmonic_number = accelerator.harmonic_number; -// double L0 = accelerator.get_length(); -// double s = 0.0; -// double accum_s = 0.0; -// unsigned int last_cav_idx = 0; -// for (unsigned int i = 0; i < accelerator.lattice.size(); i++){ -// auto elem2 = accelerator.lattice[i]; -// if (elem2 == elem){ -// accum_s += s; -// } -// if (elem2.frequency != 0.0){ -// s = 0.0; -// last_cav_idx = i; -// } -// s += elem2.length; -// } -// double ddl = (light_speed*harmonic_number/frf - L0); -// if (elem.length == 0) { -// T &de = pos.de, &dl = pos.dl; -// dl -= ddl*(accum_s/L0); -// de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); -// if (elem == accelerator.lattice[last_cav_idx]){ -// dl -= ddl*(s/L0); -// } -// } else { -// T &rx = pos.rx, &px = pos.px; -// T &ry = pos.ry, &py = pos.py; -// T &de = pos.de, &dl = pos.dl; -// // drift half length -// T pnorm = 1 / (1 + de); -// T norml = (0.5 * elem.length) * pnorm; -// rx += norml * px; -// ry += norml * py; -// dl += 0.5 * norml * pnorm * (px*px + py*py); -// // longitudinal momentum kick -// dl -= ddl*(accum_s/L0); -// de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); -// if (elem == accelerator.lattice[last_cav_idx]){ -// dl -= ddl*(s/L0); -// } -// // drift half length -// pnorm = 1.0 / (1.0 + de); -// norml = (0.5 * elem.length) * pnorm; -// rx += norml * px; -// ry += norml * py; -// dl += 0.5 * norml * pnorm * (px*px + py*py); -// } -// local_2_global(pos, elem); -// return Status::success; -// } +template +Status::type pm_cavity_1_pass(Pos &pos, const Element &elem, + const Accelerator& accelerator, const double nturn) { + + if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); + + global_2_local(pos, elem); + double nv = elem.voltage / accelerator.energy; + double frf = elem.frequency; + double harmonic_number = accelerator.harmonic_number; + double L0 = accelerator.get_length(); + // double time_aware_frac = accelerator.get_time_aware_frac(); + double ddl = (light_speed*harmonic_number/frf - L0)* nturn; + if (elem.length == 0) { + T &de = pos.de, &dl = pos.dl; + de += -nv * sin(TWOPI*frf *(dl-ddl)/ light_speed - elem.phase_lag); + } else { + T &rx = pos.rx, &px = pos.px; + T &ry = pos.ry, &py = pos.py; + T &de = pos.de, &dl = pos.dl; + // drift half length + T pnorm = 1 / (1 + de); + T norml = (0.5 * elem.length) * pnorm; + rx += norml * px; + ry += norml * py; + dl += 0.5 * norml * pnorm * (px*px + py*py); + // longitudinal momentum kick + de += -nv * sin(TWOPI*frf *(dl-ddl)/ light_speed - elem.phase_lag); + // drift half length + pnorm = 1.0 / (1.0 + de); + norml = (0.5 * elem.length) * pnorm; + rx += norml * px; + ry += norml * py; + dl += 0.5 * norml * pnorm * (px*px + py*py); + } + local_2_global(pos, elem); + return Status::success; +} template -Status::type pm_cavity_1comp_pass(Pos &pos, const Element &elem, +Status::type pm_cavity_2_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); @@ -540,7 +519,68 @@ Status::type pm_cavity_1comp_pass(Pos &pos, const Element &elem, } template -Status::type pm_cavity_1frac_pass(Pos &pos, const Element &elem, +Status::type pm_cavity_3_pass(Pos &pos, const Element &elem, + const Accelerator& accelerator) { + + if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); + + global_2_local(pos, elem); + double nv = elem.voltage / accelerator.energy; + double frf = elem.frequency; + double harmonic_number = accelerator.harmonic_number; + double L0 = accelerator.get_length(); + double s = 0.0; + double accum_s = 0.0; + unsigned int last_cav_idx = 0; + for (unsigned int i = 0; i < accelerator.lattice.size(); i++){ + auto elem2 = accelerator.lattice[i]; + if (elem2 == elem){ + accum_s += s; + } + if (elem2.frequency != 0.0){ + s = 0.0; + last_cav_idx = i; + } + s += elem2.length; + } + double ddl = (light_speed*harmonic_number/frf - L0); + // std::cout << "accum_s = " << accum_s << ", s = " << s << ", last_cav_index: " << last_cav_idx << std::endl; + if (elem.length == 0) { + T &de = pos.de, &dl = pos.dl; + dl -= ddl*(accum_s/L0); + de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); + if (elem == accelerator.lattice[last_cav_idx]){ + dl -= ddl*(s/L0); + } + } else { + T &rx = pos.rx, &px = pos.px; + T &ry = pos.ry, &py = pos.py; + T &de = pos.de, &dl = pos.dl; + // drift half length + T pnorm = 1 / (1 + de); + T norml = (0.5 * elem.length) * pnorm; + rx += norml * px; + ry += norml * py; + dl += 0.5 * norml * pnorm * (px*px + py*py); + // longitudinal momentum kick + dl -= ddl*(accum_s/L0); + de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); + if (elem == accelerator.lattice[last_cav_idx]){ + dl -= ddl*(s/L0); + } + // drift half length + pnorm = 1.0 / (1.0 + de); + norml = (0.5 * elem.length) * pnorm; + rx += norml * px; + ry += norml * py; + dl += 0.5 * norml * pnorm * (px*px + py*py); + } + local_2_global(pos, elem); + return Status::success; +} + +template +Status::type pm_cavity_4_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); @@ -581,21 +621,16 @@ Status::type pm_cavity_1frac_pass(Pos &pos, const Element &elem, } template -Status::type pm_cavity_2_pass(Pos &pos, const Element &elem, - const Accelerator& accelerator, const double nturn) { +Status::type pm_cavity_5_pass(Pos &pos, const Element &elem, + const Accelerator& accelerator) { if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); global_2_local(pos, elem); double nv = elem.voltage / accelerator.energy; - double frf = elem.frequency; - double harmonic_number = accelerator.harmonic_number; - double L0 = accelerator.get_length(); - // double time_aware_frac = accelerator.get_time_aware_frac(); - double ddl = (light_speed*harmonic_number/frf - L0)* nturn;// * time_aware_frac ; if (elem.length == 0) { T &de = pos.de, &dl = pos.dl; - de += -nv * sin(TWOPI*frf *(dl-ddl)/ light_speed - elem.phase_lag); + de += -nv * sin(TWOPI*elem.frequency * dl/ light_speed - elem.phase_lag); } else { T &rx = pos.rx, &px = pos.px; T &ry = pos.ry, &py = pos.py; @@ -607,7 +642,7 @@ Status::type pm_cavity_2_pass(Pos &pos, const Element &elem, ry += norml * py; dl += 0.5 * norml * pnorm * (px*px + py*py); // longitudinal momentum kick - de += -nv * sin(TWOPI*frf *(dl-ddl)/ light_speed - elem.phase_lag); + de += -nv * sin(TWOPI*elem.frequency*dl/light_speed - elem.phase_lag); // drift half length pnorm = 1.0 / (1.0 + de); norml = (0.5 * elem.length) * pnorm; diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index c9df57f..aa3a7b5 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -97,14 +97,20 @@ 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_cavity_1comp_pass: - if ((status = pm_cavity_1comp_pass(orig_pos, el, accelerator)) != Status::success) return status; - break; - case PassMethod::pm_cavity_1frac_pass: - if ((status = pm_cavity_1frac_pass(orig_pos, el, accelerator)) != Status::success) return status; + case PassMethod::pm_cavity_1_pass: + if ((status = pm_cavity_1_pass(orig_pos, el, accelerator, nturn)) != Status::success) return status; break; case PassMethod::pm_cavity_2_pass: - if ((status = pm_cavity_2_pass(orig_pos, el, accelerator, nturn)) != Status::success) return status; + if ((status = pm_cavity_2_pass(orig_pos, el, accelerator)) != Status::success) return status; + break; + case PassMethod::pm_cavity_3_pass: + if ((status = pm_cavity_3_pass(orig_pos, el, accelerator)) != Status::success) return status; + break; + case PassMethod::pm_cavity_4_pass: + if ((status = pm_cavity_4_pass(orig_pos, el, accelerator)) != Status::success) return status; + break; + case PassMethod::pm_cavity_5_pass: + if ((status = pm_cavity_5_pass(orig_pos, el, accelerator)) != Status::success) return status; break; default: return Status::passmethod_not_defined; @@ -176,6 +182,7 @@ Status::type track_linepass ( // syntactic-sugar for read-only access to element object parameters const Element& element = line[element_offset]; + if (element.frequency != 0.0) { timectc += tawfrac; } @@ -185,7 +192,7 @@ Status::type track_linepass ( status = track_elementpass (element, orig_pos, accelerator, nturn); - if ((timectc >= 0.98) && (element.pass_method == PassMethod::pm_cavity_0_pass)) { + if ((timectc >= 0.98) && (element.pass_method == PassMethod::pm_cavity_5_pass)) { orig_pos.dl -= (light_speed*accelerator.harmonic_number/element.frequency - accelerator.get_length()); timectc = 0.0; } diff --git a/src/accelerator.cpp b/src/accelerator.cpp index 1e5ab9a..79cee26 100644 --- a/src/accelerator.cpp +++ b/src/accelerator.cpp @@ -34,13 +34,19 @@ inline bool check_time_aware_pm(const int i){ case PassMethod::pm_cavity_0_pass: flag = true; break; - case PassMethod::pm_cavity_1comp_pass: + case PassMethod::pm_cavity_1_pass: flag = true; break; - case PassMethod::pm_cavity_1frac_pass: + case PassMethod::pm_cavity_2_pass: flag = true; break; - case PassMethod::pm_cavity_2_pass: + case PassMethod::pm_cavity_3_pass: + flag = true; + break; + case PassMethod::pm_cavity_4_pass: + flag = true; + break; + case PassMethod::pm_cavity_5_pass: flag = true; break; default: From 8be19a5190d846fdfae408f3ecf37819e4ee2317 Mon Sep 17 00:00:00 2001 From: vitor Date: Mon, 10 Jun 2024 17:06:57 -0300 Subject: [PATCH 07/20] saving --- include/trackcpp/accelerator.h | 2 +- include/trackcpp/auxiliary.h | 56 ++----- include/trackcpp/passmethods.h | 7 +- include/trackcpp/passmethods.hpp | 252 ++++--------------------------- include/trackcpp/tracking.h | 192 ++++++++++++++++++----- python_package/interface.cpp | 4 +- python_package/interface.h | 3 +- src/accelerator.cpp | 44 +----- src/flat_file.cpp | 2 +- src/kicktable.cpp | 2 +- src/lattice.cpp | 2 +- src/tracking.cpp | 92 +---------- 12 files changed, 204 insertions(+), 454 deletions(-) diff --git a/include/trackcpp/accelerator.h b/include/trackcpp/accelerator.h index 98c424f..8a4922d 100644 --- a/include/trackcpp/accelerator.h +++ b/include/trackcpp/accelerator.h @@ -39,7 +39,7 @@ class Accelerator { bool operator!=(const Accelerator& o) const { return !(*this == o); }; bool isequal(const Accelerator& a) const { return *this == a; } // necessary for python_package double get_length() const; - double get_time_aware_frac() const; + // double get_time_aware_frac() const; friend std::ostream& operator<< (std::ostream &out, const Accelerator& a); }; diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 9fa6d9e..8058f5e 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -31,64 +31,36 @@ class PassMethodsClass { static const int pm_str_mpole_symplectic4_pass = 2; static const int pm_bnd_mpole_symplectic4_pass = 3; static const int pm_corrector_pass = 4; - static const int pm_cavity_0_pass = 5; + static const int pm_cavity_pass = 5; static const int pm_thinquad_pass = 6; static const int pm_thinsext_pass = 7; 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_cavity_1_pass = 11; - static const int pm_cavity_2_pass = 12; - static const int pm_cavity_3_pass = 13; - static const int pm_cavity_4_pass = 14; - static const int pm_cavity_5_pass = 15; - static const int pm_nr_pms = 16; // counter for number of passmethods + static const int pm_nr_pms = 11; // counter for number of passmethods PassMethodsClass() { passmethods.push_back("identity_pass"); passmethods.push_back("drift_pass"); passmethods.push_back("str_mpole_symplectic4_pass"); passmethods.push_back("bnd_mpole_symplectic4_pass"); passmethods.push_back("corrector_pass"); - passmethods.push_back("cavity_0_pass"); + passmethods.push_back("cavity_pass"); passmethods.push_back("thinquad_pass"); passmethods.push_back("thinsext_pass"); passmethods.push_back("kicktable_pass"); passmethods.push_back("matrix_pass"); passmethods.push_back("drift_g2l_pass"); - passmethods.push_back("cavity_1_pass"); - passmethods.push_back("cavity_2_pass"); - passmethods.push_back("cavity_3_pass"); - passmethods.push_back("cavity_4_pass"); - passmethods.push_back("cavity_5_pass"); } int size() const { return passmethods.size(); } std::string operator[](const int i) const { return passmethods[i]; } bool is_time_aware_pm(const int i) const { - bool flag = false; switch (i) { - case pm_cavity_0_pass: - flag = true; - break; - case pm_cavity_1_pass: - flag = true; - break; - case pm_cavity_2_pass: - flag = true; - break; - case pm_cavity_3_pass: - flag = true; - break; - case pm_cavity_4_pass: - flag = true; - break; - case pm_cavity_5_pass: - flag = true; - break; + case pm_cavity_pass: + return true; default: - break; + return false; } - return flag; }; private: std::vector passmethods; @@ -102,18 +74,13 @@ struct PassMethod { pm_str_mpole_symplectic4_pass = 2, pm_bnd_mpole_symplectic4_pass = 3, pm_corrector_pass = 4, - pm_cavity_0_pass = 5, + pm_cavity_pass = 5, pm_thinquad_pass = 6, pm_thinsext_pass = 7, pm_kickmap_pass = 8, pm_matrix_pass = 9, pm_drift_g2l_pass = 10, - pm_cavity_1_pass = 11, - pm_cavity_2_pass = 12, - pm_cavity_3_pass = 13, - pm_cavity_4_pass = 14, - pm_cavity_5_pass = 15, - pm_nr_pms = 16, + pm_nr_pms = 11, }; }; @@ -125,17 +92,12 @@ const std::vector pm_dict = { "str_mpole_symplectic4_pass", "bnd_mpole_symplectic4_pass", "corrector_pass", - "cavity_0_pass", + "cavity_pass", "thinquad_pass", "thinsext_pass", "kicktable_pass", "matrix_pass", "drift_g2l_pass", - "cavity_1_pass", - "cavity_2_pass", - "cavity_3_pass", - "cavity_4_pass", - "cavity_5_pass", }; struct RadiationState { diff --git a/include/trackcpp/passmethods.h b/include/trackcpp/passmethods.h index 6aa0a3a..2ab02a3 100644 --- a/include/trackcpp/passmethods.h +++ b/include/trackcpp/passmethods.h @@ -60,12 +60,7 @@ template Status::type pm_drift_pass (Pos &pos, c template Status::type pm_str_mpole_symplectic4_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_bnd_mpole_symplectic4_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_corrector_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); -template Status::type pm_cavity_0_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); -template Status::type pm_cavity_1_pass (Pos &pos, const Element &elem, const Accelerator& accelerator, const double nturn = 0.0); -template Status::type pm_cavity_2_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); -template Status::type pm_cavity_3_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); -template Status::type pm_cavity_4_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); -template Status::type pm_cavity_5_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_cavity_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_thinquad_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_thinsext_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_kickmap_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 056bf7e..0578331 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -392,81 +392,9 @@ Status::type pm_corrector_pass(Pos &pos, const Element &elem, return Status::success; } -template -Status::type pm_cavity_0_pass(Pos &pos, const Element &elem, - const Accelerator& accelerator) { - - if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); - - global_2_local(pos, elem); - double nv = elem.voltage / accelerator.energy; - if (elem.length == 0) { - T &de = pos.de, &dl = pos.dl; - de += -nv * sin(TWOPI*elem.frequency * dl/ light_speed - elem.phase_lag); - } else { - T &rx = pos.rx, &px = pos.px; - T &ry = pos.ry, &py = pos.py; - T &de = pos.de, &dl = pos.dl; - // drift half length - T pnorm = 1 / (1 + de); - T norml = (0.5 * elem.length) * pnorm; - rx += norml * px; - ry += norml * py; - dl += 0.5 * norml * pnorm * (px*px + py*py); - // longitudinal momentum kick - de += -nv * sin(TWOPI*elem.frequency*dl/light_speed - elem.phase_lag); - // drift half length - pnorm = 1.0 / (1.0 + de); - norml = (0.5 * elem.length) * pnorm; - rx += norml * px; - ry += norml * py; - dl += 0.5 * norml * pnorm * (px*px + py*py); - } - local_2_global(pos, elem); - return Status::success; -} - -template -Status::type pm_cavity_1_pass(Pos &pos, const Element &elem, - const Accelerator& accelerator, const double nturn) { - - if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); - - global_2_local(pos, elem); - double nv = elem.voltage / accelerator.energy; - double frf = elem.frequency; - double harmonic_number = accelerator.harmonic_number; - double L0 = accelerator.get_length(); - // double time_aware_frac = accelerator.get_time_aware_frac(); - double ddl = (light_speed*harmonic_number/frf - L0)* nturn; - if (elem.length == 0) { - T &de = pos.de, &dl = pos.dl; - de += -nv * sin(TWOPI*frf *(dl-ddl)/ light_speed - elem.phase_lag); - } else { - T &rx = pos.rx, &px = pos.px; - T &ry = pos.ry, &py = pos.py; - T &de = pos.de, &dl = pos.dl; - // drift half length - T pnorm = 1 / (1 + de); - T norml = (0.5 * elem.length) * pnorm; - rx += norml * px; - ry += norml * py; - dl += 0.5 * norml * pnorm * (px*px + py*py); - // longitudinal momentum kick - de += -nv * sin(TWOPI*frf *(dl-ddl)/ light_speed - elem.phase_lag); - // drift half length - pnorm = 1.0 / (1.0 + de); - norml = (0.5 * elem.length) * pnorm; - rx += norml * px; - ry += norml * py; - dl += 0.5 * norml * pnorm * (px*px + py*py); - } - local_2_global(pos, elem); - return Status::success; -} template -Status::type pm_cavity_2_pass(Pos &pos, const Element &elem, +Status::type pm_cavity_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); @@ -474,84 +402,31 @@ Status::type pm_cavity_2_pass(Pos &pos, const Element &elem, global_2_local(pos, elem); double nv = elem.voltage / accelerator.energy; double frf = elem.frequency; - double harmonic_number = accelerator.harmonic_number; - double L0 = accelerator.get_length(); - unsigned int last_cav_idx = 0; - for (unsigned int i = 0; i < accelerator.lattice.size(); i++){ - auto elem2 = accelerator.lattice[i]; - if (elem2.frequency != 0.0){ - last_cav_idx = i; - } - } - double ddl = (light_speed*harmonic_number/frf - L0); - if (elem.length == 0) { - T &de = pos.de, &dl = pos.dl; - - de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); - if (elem == accelerator.lattice[last_cav_idx]){ - dl -= ddl; - } - } else { - T &rx = pos.rx, &px = pos.px; - T &ry = pos.ry, &py = pos.py; - T &de = pos.de, &dl = pos.dl; - // drift half length - T pnorm = 1 / (1 + de); - T norml = (0.5 * elem.length) * pnorm; - rx += norml * px; - ry += norml * py; - dl += 0.5 * norml * pnorm * (px*px + py*py); - // longitudinal momentum kick - - de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); - if (elem == accelerator.lattice[last_cav_idx]){ - dl -= ddl; - } - // drift half length - pnorm = 1.0 / (1.0 + de); - norml = (0.5 * elem.length) * pnorm; - rx += norml * px; - ry += norml * py; - dl += 0.5 * norml * pnorm * (px*px + py*py); - } - local_2_global(pos, elem); - return Status::success; -} - -template -Status::type pm_cavity_3_pass(Pos &pos, const Element &elem, - const Accelerator& accelerator) { - - if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); - - global_2_local(pos, elem); - double nv = elem.voltage / accelerator.energy; - double frf = elem.frequency; - double harmonic_number = accelerator.harmonic_number; - double L0 = accelerator.get_length(); - double s = 0.0; - double accum_s = 0.0; - unsigned int last_cav_idx = 0; - for (unsigned int i = 0; i < accelerator.lattice.size(); i++){ - auto elem2 = accelerator.lattice[i]; - if (elem2 == elem){ - accum_s += s; - } - if (elem2.frequency != 0.0){ - s = 0.0; - last_cav_idx = i; - } - s += elem2.length; - } - double ddl = (light_speed*harmonic_number/frf - L0); + // double harmonic_number = accelerator.harmonic_number; + // double L0 = accelerator.get_length(); + // double s = 0.0; + // double accum_s = 0.0; + // unsigned int last_cav_idx = 0; + // for (unsigned int i = 0; i < accelerator.lattice.size(); i++){ + // auto elem2 = accelerator.lattice[i]; + // if (elem2 == elem){ + // accum_s += s; + // } + // if (elem2.frequency != 0.0){ + // s = 0.0; + // last_cav_idx = i; + // } + // s += elem2.length; + // } + // double ddl = (light_speed*harmonic_number/frf - L0); // std::cout << "accum_s = " << accum_s << ", s = " << s << ", last_cav_index: " << last_cav_idx << std::endl; if (elem.length == 0) { T &de = pos.de, &dl = pos.dl; - dl -= ddl*(accum_s/L0); + // dl -= ddl*(accum_s/L0); de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); - if (elem == accelerator.lattice[last_cav_idx]){ - dl -= ddl*(s/L0); - } + // if (elem == accelerator.lattice[last_cav_idx]){ + // dl -= ddl*(s/L0); + // } } else { T &rx = pos.rx, &px = pos.px; T &ry = pos.ry, &py = pos.py; @@ -563,86 +438,11 @@ Status::type pm_cavity_3_pass(Pos &pos, const Element &elem, ry += norml * py; dl += 0.5 * norml * pnorm * (px*px + py*py); // longitudinal momentum kick - dl -= ddl*(accum_s/L0); + // dl -= ddl*(accum_s/L0); de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); - if (elem == accelerator.lattice[last_cav_idx]){ - dl -= ddl*(s/L0); - } - // drift half length - pnorm = 1.0 / (1.0 + de); - norml = (0.5 * elem.length) * pnorm; - rx += norml * px; - ry += norml * py; - dl += 0.5 * norml * pnorm * (px*px + py*py); - } - local_2_global(pos, elem); - return Status::success; -} - -template -Status::type pm_cavity_4_pass(Pos &pos, const Element &elem, - const Accelerator& accelerator) { - - if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); - - global_2_local(pos, elem); - double nv = elem.voltage / accelerator.energy; - double frf = elem.frequency; - double harmonic_number = accelerator.harmonic_number; - double L0 = accelerator.get_length(); - double time_aware_frac = accelerator.get_time_aware_frac(); - double ddl = (light_speed*harmonic_number/frf - L0)*time_aware_frac; - if (elem.length == 0) { - T &de = pos.de, &dl = pos.dl; - dl -= ddl; - de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); - } else { - T &rx = pos.rx, &px = pos.px; - T &ry = pos.ry, &py = pos.py; - T &de = pos.de, &dl = pos.dl; - // drift half length - T pnorm = 1 / (1 + de); - T norml = (0.5 * elem.length) * pnorm; - rx += norml * px; - ry += norml * py; - dl += 0.5 * norml * pnorm * (px*px + py*py); - // longitudinal momentum kick - dl -= ddl; - de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); - // drift half length - pnorm = 1.0 / (1.0 + de); - norml = (0.5 * elem.length) * pnorm; - rx += norml * px; - ry += norml * py; - dl += 0.5 * norml * pnorm * (px*px + py*py); - } - local_2_global(pos, elem); - return Status::success; -} - -template -Status::type pm_cavity_5_pass(Pos &pos, const Element &elem, - const Accelerator& accelerator) { - - if (not accelerator.cavity_on) return pm_drift_pass(pos, elem, accelerator); - - global_2_local(pos, elem); - double nv = elem.voltage / accelerator.energy; - if (elem.length == 0) { - T &de = pos.de, &dl = pos.dl; - de += -nv * sin(TWOPI*elem.frequency * dl/ light_speed - elem.phase_lag); - } else { - T &rx = pos.rx, &px = pos.px; - T &ry = pos.ry, &py = pos.py; - T &de = pos.de, &dl = pos.dl; - // drift half length - T pnorm = 1 / (1 + de); - T norml = (0.5 * elem.length) * pnorm; - rx += norml * px; - ry += norml * py; - dl += 0.5 * norml * pnorm * (px*px + py*py); - // longitudinal momentum kick - de += -nv * sin(TWOPI*elem.frequency*dl/light_speed - elem.phase_lag); + // if (elem == accelerator.lattice[last_cav_idx]){ + // dl -= ddl*(s/L0); + // } // drift half length pnorm = 1.0 / (1.0 + de); norml = (0.5 * elem.length) * pnorm; diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index aa3a7b5..64ba801 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -44,10 +44,6 @@ Status::type track_findorbit6( Accelerator& accelerator, std::vector >& closed_orbit, const Pos& fixed_point_guess = Pos(0)); -Status::type track_findorbit6_dct( - Accelerator& accelerator, std::vector >& closed_orbit, - const Pos& fixed_point_guess = Pos(0), const double dct = 0.0); - Pos linalg_solve4_posvec( const std::vector >& M, const Pos& b); @@ -59,7 +55,7 @@ template Status::type track_elementpass ( const Element& el, // element through which to track particle Pos &orig_pos, // initial electron coordinates - const Accelerator& accelerator, const double nturn = 0.0) { + const Accelerator& accelerator) { Status::type status = Status::success; @@ -79,8 +75,8 @@ Status::type track_elementpass ( case PassMethod::pm_corrector_pass: if ((status = pm_corrector_pass(orig_pos, el, accelerator)) != Status::success) return status; break; - case PassMethod::pm_cavity_0_pass: - if ((status = pm_cavity_0_pass(orig_pos, el, accelerator)) != Status::success) return status; + case PassMethod::pm_cavity_pass: + if ((status = pm_cavity_pass(orig_pos, el, accelerator)) != Status::success) return status; break; case PassMethod::pm_thinquad_pass: if ((status = pm_thinquad_pass(orig_pos, el, accelerator)) != Status::success) return status; @@ -97,21 +93,6 @@ 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_cavity_1_pass: - if ((status = pm_cavity_1_pass(orig_pos, el, accelerator, nturn)) != Status::success) return status; - break; - case PassMethod::pm_cavity_2_pass: - if ((status = pm_cavity_2_pass(orig_pos, el, accelerator)) != Status::success) return status; - break; - case PassMethod::pm_cavity_3_pass: - if ((status = pm_cavity_3_pass(orig_pos, el, accelerator)) != Status::success) return status; - break; - case PassMethod::pm_cavity_4_pass: - if ((status = pm_cavity_4_pass(orig_pos, el, accelerator)) != Status::success) return status; - break; - case PassMethod::pm_cavity_5_pass: - if ((status = pm_cavity_5_pass(orig_pos, el, accelerator)) != Status::success) return status; - break; default: return Status::passmethod_not_defined; } @@ -124,7 +105,7 @@ template Status::type track_elementpass ( const Element& el, // element through which to track particle std::vector >& orig_pos, // initial electron coordinates - const Accelerator& accelerator, const double nturn) { + const Accelerator& accelerator) { Status::type status = Status::success; @@ -151,7 +132,102 @@ Status::type track_elementpass ( // pos: Pos vector of particles' final coordinates (or trajetory) // element_offset: in case of problems with passmethods, '*element_offset' is the index of the corresponding element // RETURN: status do tracking (see 'auxiliary.h') - +// template +// Status::type track_linepass ( +// const Accelerator& accelerator, +// Pos& orig_pos, // initial electron coordinates +// std::vector >& pos, // vector with tracked electron coordinates at start of every element and at the end of last one. +// unsigned int& element_offset, // index of starting element for tracking +// Plane::type& lost_plane, // return plane in which particle was lost, if the case. +// std::vector& indices, double nturn = 0.0) {// indices to return; + +// Status::type status = Status::success; +// lost_plane = Plane::no_plane; + +// const std::vector& line = accelerator.lattice; +// int nr_elements = line.size(); + +// //pos.clear(); other functions assume pos is not clearedin linepass! +// pos.reserve(pos.size() + indices.size()); + +// // create vector of booleans help determine when to store position +// std::vector indcs; +// indcs.reserve(nr_elements+1); +// for (unsigned int i=0; i<=nr_elements; ++i) indcs[i] = false; +// for (auto&& i: indices) if (i<=nr_elements) indcs[i] = true; + +// for(int i=0; i= element.hmax))) { +// lost_plane = Plane::x; +// status = Status::particle_lost; +// } +// if (((ry <= element.vmin) or (ry >= element.vmax))) { +// if (status != Status::particle_lost) { +// lost_plane = Plane::y; +// status = Status::particle_lost; +// } else { +// lost_plane = Plane::xy; +// } +// } +// } else if (get_norm_amp_in_vchamber(element, rx, ry) > 1) { +// // lost in rhombus, elliptic and all finite p-norm shapes +// lost_plane = Plane::xy; +// status = Status::particle_lost; +// } +// } + +// if (status != Status::success) { +// // fill rest of vector with nans +// for(int ii=i+1; ii<=nr_elements; ++ii) { +// if (indcs[ii]) pos.emplace_back( +// nan(""),nan(""),nan(""),nan(""),nan(""),nan("")); +// } +// return status; +// } +// // moves to next element index +// element_offset = (element_offset + 1) % nr_elements; +// } + +// // stores final particle position at the end of the line +// if (indcs[nr_elements]) pos.push_back(orig_pos); + +// return (status == Status::success) ? status: Status::particle_lost; +// } + +// accelerator_length, ddl, TAW_positions template Status::type track_linepass ( const Accelerator& accelerator, @@ -159,12 +235,14 @@ Status::type track_linepass ( std::vector >& pos, // vector with tracked electron coordinates at start of every element and at the end of last one. unsigned int& element_offset, // index of starting element for tracking Plane::type& lost_plane, // return plane in which particle was lost, if the case. - std::vector& indices, double nturn = 0.0) {// indices to return; + std::vector& indices, // indices to return; + double line_length, + double ddl, + std::vector time_aware_element_indices, + std::vector time_aware_element_positions) { Status::type status = Status::success; lost_plane = Plane::no_plane; - double timectc = 0; - double tawfrac = accelerator.get_time_aware_frac(); const std::vector& line = accelerator.lattice; int nr_elements = line.size(); @@ -178,23 +256,24 @@ Status::type track_linepass ( for (unsigned int i=0; i<=nr_elements; ++i) indcs[i] = false; for (auto&& i: indices) if (i<=nr_elements) indcs[i] = true; + unsigned int TAW_pivot = 0; for(int i=0; i= 0.98) && (element.pass_method == PassMethod::pm_cavity_5_pass)) { - orig_pos.dl -= (light_speed*accelerator.harmonic_number/element.frequency - accelerator.get_length()); - timectc = 0.0; + if (element_offset == time_aware_element_indices.back()) { + orig_pos.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; } const T& rx = orig_pos.rx; @@ -286,7 +365,11 @@ Status::type track_linepass ( std::vector >& pos, // vector with tracked electron coordinates at start of every element and at the end of last one. unsigned int& element_offset, // index of starting element for tracking Plane::type& lost_plane, // return plane in which particle was lost, if the case. - bool trajectory, double nturn = 0.0) { // whether function should return coordinates at all elements + bool trajectory, // whether function should return coordinates at all elements + double line_length, + double ddl, + std::vector time_aware_element_indices, + std::vector time_aware_element_positions) { std::vector indices; unsigned int nr_elements = accelerator.lattice.size(); @@ -299,7 +382,8 @@ Status::type track_linepass ( } return track_linepass ( - accelerator, orig_pos, pos, element_offset, lost_plane, indices, nturn); + accelerator, orig_pos, pos, element_offset, lost_plane, indices, + line_length, ddl, time_aware_element_indices, time_aware_element_positions); } @@ -311,7 +395,11 @@ Status::type track_linepass ( unsigned int element_offset, std::vector& lost_plane, std::vector& lost_element, - std::vector& indices, double nturn = 0.0) { + std::vector& indices, + double line_length, + double ddl, + std::vector time_aware_element_indices, + std::vector time_aware_element_positions) { int nr_elements = accelerator.lattice.size(); Status::type status = Status::success; @@ -327,7 +415,8 @@ Status::type track_linepass ( unsigned int le = element_offset; status2 = track_linepass ( - accelerator, orig_pos[i], final_pos, le, lp, indices, nturn); + accelerator, orig_pos[i], final_pos, le, lp, indices, + line_length, ddl, time_aware_element_indices, time_aware_element_positions); if (status2 != Status::success) status = status2; @@ -370,7 +459,30 @@ Status::type track_ringpass ( Status::type status = Status::success; std::vector > final_pos; - double nturn = 0.0; + + // for longitudinal kick before RF cavities + std::vector TAW_positions = {0.0, }; + std::vector TAW_indices = {}; + unsigned int eoff = 0; eoff += element_offset; + size_t nr_elements = accelerator.lattice.size(); + PassMethodsClass PMClass = PassMethodsClass(); + double s_pos = 0.0; + for (size_t i = 0; i < nr_elements; i+=) { + auto elem = accelerator.lattice[eoff]; + if (PMClass.is_time_aware_pm(elem.pass_method)) { + s_pos += (elem.length*0.5); + TAW_positions.push_back(s_pos); + TAW_indices.push_back(eoff); + s_pos += (elem.length*0.5); + } + else { + s_pos += elem.length; + } + eoff = (eoff + 1) % nr_elements; + } + TAW_positions.push_back(s_pos); + double accelerator_length = s_pos; + double ddl = light_speed*accelerator.harmonic_number/element.frequency - accelerator_length; if (trajectory) pos.reserve(nr_turns+1); @@ -379,7 +491,7 @@ Status::type track_ringpass ( // stores trajectory at beggining of each turn if (trajectory) pos.push_back(orig_pos); - if ((status = track_linepass (accelerator, orig_pos, final_pos, element_offset, lost_plane, false, nturn)) != Status::success) { + if ((status = track_linepass (accelerator, orig_pos, final_pos, element_offset, lost_plane, false, accelerator_length, ddl, TAW_positions)) != Status::success) { // fill last of vector with nans pos.emplace_back( diff --git a/python_package/interface.cpp b/python_package/interface.cpp index a74ad12..2717bc1 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -21,7 +21,7 @@ Status::type track_elementpass_wrapper ( const Element& el, double *pos, int n1, int n2, - const Accelerator& accelerator, const double nturn) { + const Accelerator& accelerator) { std::vector > post; @@ -33,7 +33,7 @@ Status::type track_elementpass_wrapper ( pos[4*n2 + i], pos[5*n2 + i]); } Status::type status = track_elementpass( - el, post, accelerator, nturn); + el, post, accelerator); for (unsigned int i=0; i indices; std::vector< unsigned int > lost_plane; std::vector< unsigned int > lost_element; - double nturn = 0.0; }; struct RingPassArgs : public LinePassArgs { @@ -53,7 +52,7 @@ struct String { Status::type track_elementpass_wrapper ( const Element& el, double *pos, int n1, int n2, - const Accelerator& accelerator, const double nturn = 0.0); + const Accelerator& accelerator); Status::type track_linepass_wrapper ( const Accelerator& accelerator, diff --git a/src/accelerator.cpp b/src/accelerator.cpp index 79cee26..85e1ebb 100644 --- a/src/accelerator.cpp +++ b/src/accelerator.cpp @@ -27,42 +27,14 @@ double Accelerator::get_length() const { return length; } -inline bool check_time_aware_pm(const int i){ - bool flag = false; - switch (i) - { - case PassMethod::pm_cavity_0_pass: - flag = true; - break; - case PassMethod::pm_cavity_1_pass: - flag = true; - break; - case PassMethod::pm_cavity_2_pass: - flag = true; - break; - case PassMethod::pm_cavity_3_pass: - flag = true; - break; - case PassMethod::pm_cavity_4_pass: - flag = true; - break; - case PassMethod::pm_cavity_5_pass: - flag = true; - break; - default: - break; - } - return flag; -} - -double Accelerator::get_time_aware_frac() const { - double count = 0.0; - for(auto i=0; i> e.voltage >> e.frequency >> hnumber >> energy; e.voltage *= energy; e.frequency *= light_speed / (2*M_PI); diff --git a/src/kicktable.cpp b/src/kicktable.cpp index aace606..0034cba 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -104,7 +104,7 @@ Status::type Kicktable::load_from_file(const std::string& filename_) { } int add_kicktable(const std::string& filename) { - + // looks through vector of kicktables... for(unsigned int i=0; i& lattice, const std::string& state) { for(unsigned int i=0; i cav_idx = latt_findcells_frequency(the_ring, 0, true); - double frf = the_ring[cav_idx[0]].frequency; - double fixedpoint = light_speed*((1.0*accelerator.harmonic_number)/frf - T0); // temporary vectors and matrices std::vector > co(7,0); @@ -160,8 +154,6 @@ Status::type track_findorbit6( std::vector > D(7,0); std::vector > M(6,0); Pos dco(1.0,1.0,1.0,1.0,1.0,1.0); - Pos theta(0.0,0.0,0.0,0.0,0.0,0.0); - theta.dl = fixedpoint; matrix6_set_identity_posvec(D, delta); int nr_iter = 0; @@ -194,7 +186,7 @@ Status::type track_findorbit6( M[4] = (co2[4] - Rf) / delta; M[5] = (co2[5] - Rf) / delta; - Pos b = Rf - Ri - theta; + Pos b = Rf - Ri; std::vector > M_1(6,0); matrix6_set_identity_posvec(M_1); M_1 = M_1 - M; @@ -220,88 +212,6 @@ Status::type track_findorbit6( } -Status::type track_findorbit6_dct( - Accelerator& accelerator, - std::vector >& closed_orbit, - const Pos& fixed_point_guess, const double dct) { - - const std::vector& the_ring = accelerator.lattice; - - double delta = 1e-9; // [m],[rad],[dE/E] - double tolerance = 2.22044604925e-14; - int max_nr_iters = 50; - - const int radsts = accelerator.radiation_on; - if (radsts == RadiationState::full){ - accelerator.radiation_on = RadiationState::damping; - } - // calcs longitudinal fixed point - - // temporary vectors and matrices - std::vector > co(7,0); - for(auto i=0; i<7; ++i) co[i] = fixed_point_guess; - std::vector > co2(7,0); - std::vector > D(7,0); - std::vector > M(6,0); - Pos dco(1.0,1.0,1.0,1.0,1.0,1.0); - Pos theta(0.0,0.0,0.0,0.0,0.0,dct); - matrix6_set_identity_posvec(D, delta); - - int nr_iter = 0; - while ((get_max(dco) > tolerance) and (nr_iter <= max_nr_iters)) { - co = co + D; - Pos Ri = co[6]; - std::vector > co2; - unsigned int element_offset = 0; - Plane::type lost_plane; - Status::type status = Status::success; - - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[0], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[1], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[2], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[3], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[4], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[5], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[6], co2, element_offset, lost_plane, false)); - - if (status != Status::success) { - return Status::findorbit_one_turn_matrix_problem; - } - - Pos Rf = co2[6]; - - M[0] = (co2[0] - Rf) / delta; - M[1] = (co2[1] - Rf) / delta; - M[2] = (co2[2] - Rf) / delta; - M[3] = (co2[3] - Rf) / delta; - M[4] = (co2[4] - Rf) / delta; - M[5] = (co2[5] - Rf) / delta; - - Pos b = Rf - Ri - theta; - std::vector > M_1(6,0); - matrix6_set_identity_posvec(M_1); - M_1 = M_1 - M; - dco = linalg_solve6_posvec(M_1, b); - co[6] = dco + Ri; - co[0] = co[6]; co[1] = co[6]; - co[2] = co[6]; co[3] = co[6]; - co[4] = co[6]; co[5] = co[6]; - nr_iter++; - } - - if (nr_iter > max_nr_iters) { - return Status::findorbit_not_converged; - } - - // propagates fixed point throught the_ring - closed_orbit.clear(); - unsigned int element_offset = 0; - Plane::type lost_plane; - track_linepass(accelerator, co[6], closed_orbit, element_offset, lost_plane, true); - accelerator.radiation_on = radsts; - return Status::success; - -} Status::type track_findorbit4( Accelerator& accelerator, From a976b4bd95a8f490178de467fbf41db875f98a04 Mon Sep 17 00:00:00 2001 From: vitor Date: Mon, 10 Jun 2024 17:26:25 -0300 Subject: [PATCH 08/20] saving2 --- include/trackcpp/tracking.h | 37 +++++++++++++++++------------------- python_package/interface.cpp | 4 +++- python_package/interface.h | 3 +++ src/elements.cpp | 2 +- 4 files changed, 24 insertions(+), 22 deletions(-) diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index 64ba801..20f5873 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -110,7 +110,7 @@ Status::type track_elementpass ( Status::type status = Status::success; for(auto&& pos: orig_pos) { - Status::type status2 = track_elementpass(el, pos, accelerator, nturn); + Status::type status2 = track_elementpass(el, pos, accelerator); if (status2 != Status::success) status = status2; } return status; @@ -236,13 +236,13 @@ Status::type track_linepass ( unsigned int& element_offset, // index of starting element for tracking Plane::type& lost_plane, // return plane in which particle was lost, if the case. std::vector& indices, // indices to return; - double line_length, - double ddl, - std::vector time_aware_element_indices, - std::vector time_aware_element_positions) { + double line_length = 0.0, + std::vector time_aware_element_indices = {0, }, + std::vector time_aware_element_positions = {0.0, }) { Status::type status = Status::success; lost_plane = Plane::no_plane; + double ddl = 0.0; const std::vector& line = accelerator.lattice; int nr_elements = line.size(); @@ -266,11 +266,12 @@ Status::type track_linepass ( if (indcs[i]) pos.push_back(orig_pos); if (element_offset == time_aware_element_indices[TAW_pivot]) { + ddl = light_speed*accelerator.harmonic_number/element.frequency - line_length; orig_pos.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; TAW_pivot++; } - status = track_elementpass (element, orig_pos, accelerator, nturn); + status = track_elementpass (element, orig_pos, accelerator); if (element_offset == time_aware_element_indices.back()) { orig_pos.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; @@ -366,10 +367,9 @@ Status::type track_linepass ( unsigned int& element_offset, // index of starting element for tracking Plane::type& lost_plane, // return plane in which particle was lost, if the case. bool trajectory, // whether function should return coordinates at all elements - double line_length, - double ddl, - std::vector time_aware_element_indices, - std::vector time_aware_element_positions) { + double line_length = 0.0, + std::vector time_aware_element_indices = {0,}, + std::vector time_aware_element_positions = {0,}) { std::vector indices; unsigned int nr_elements = accelerator.lattice.size(); @@ -383,7 +383,7 @@ Status::type track_linepass ( return track_linepass ( accelerator, orig_pos, pos, element_offset, lost_plane, indices, - line_length, ddl, time_aware_element_indices, time_aware_element_positions); + line_length, time_aware_element_indices, time_aware_element_positions); } @@ -396,10 +396,9 @@ Status::type track_linepass ( std::vector& lost_plane, std::vector& lost_element, std::vector& indices, - double line_length, - double ddl, - std::vector time_aware_element_indices, - std::vector time_aware_element_positions) { + double line_length = 0.0, + std::vector time_aware_element_indices = {0, }, + std::vector time_aware_element_positions = {0.0, }) { int nr_elements = accelerator.lattice.size(); Status::type status = Status::success; @@ -416,7 +415,7 @@ Status::type track_linepass ( status2 = track_linepass ( accelerator, orig_pos[i], final_pos, le, lp, indices, - line_length, ddl, time_aware_element_indices, time_aware_element_positions); + line_length, time_aware_element_indices, time_aware_element_positions); if (status2 != Status::success) status = status2; @@ -467,7 +466,7 @@ Status::type track_ringpass ( size_t nr_elements = accelerator.lattice.size(); PassMethodsClass PMClass = PassMethodsClass(); double s_pos = 0.0; - for (size_t i = 0; i < nr_elements; i+=) { + for (size_t i = 0; i < nr_elements; i++) { auto elem = accelerator.lattice[eoff]; if (PMClass.is_time_aware_pm(elem.pass_method)) { s_pos += (elem.length*0.5); @@ -482,7 +481,6 @@ Status::type track_ringpass ( } TAW_positions.push_back(s_pos); double accelerator_length = s_pos; - double ddl = light_speed*accelerator.harmonic_number/element.frequency - accelerator_length; if (trajectory) pos.reserve(nr_turns+1); @@ -491,7 +489,7 @@ Status::type track_ringpass ( // stores trajectory at beggining of each turn if (trajectory) pos.push_back(orig_pos); - if ((status = track_linepass (accelerator, orig_pos, final_pos, element_offset, lost_plane, false, accelerator_length, ddl, TAW_positions)) != Status::success) { + if ((status = track_linepass (accelerator, orig_pos, final_pos, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)) != Status::success) { // fill last of vector with nans pos.emplace_back( @@ -503,7 +501,6 @@ Status::type track_ringpass ( return status; } final_pos.clear(); - nturn += 1.0; } pos.push_back(orig_pos); diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 2717bc1..6c0d7ec 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -66,7 +66,9 @@ Status::type track_linepass_wrapper( args.lost_plane, args.lost_element, args.indices, - args.nturn); + args.line_length, + args.time_aware_element_indices, + args.time_aware_element_positions); for (unsigned int i=0; i indices; std::vector< unsigned int > lost_plane; std::vector< unsigned int > lost_element; + double line_length; + std::vector time_aware_element_indices; + std::vector time_aware_element_positions; }; struct RingPassArgs : public LinePassArgs { diff --git a/src/elements.cpp b/src/elements.cpp index 21540c5..5d96fb1 100644 --- a/src/elements.cpp +++ b/src/elements.cpp @@ -296,7 +296,7 @@ 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) { - element.pass_method = PassMethod::pm_cavity_0_pass; + element.pass_method = PassMethod::pm_cavity_pass; element.frequency = frequency; element.voltage = voltage; element.phase_lag = phase_lag; From d6fda162e80f0e8b95357ae9c714bd3cfa48866b Mon Sep 17 00:00:00 2001 From: vitor Date: Tue, 11 Jun 2024 11:44:55 -0300 Subject: [PATCH 09/20] working --- include/trackcpp/accelerator.h | 2 +- include/trackcpp/tracking.h | 138 ++++----------------------------- python_package/interface.cpp | 8 +- python_package/interface.h | 5 +- src/accelerator.cpp | 26 +++++++ src/optics.cpp | 13 +++- src/tracking.cpp | 65 ++++++++++++---- 7 files changed, 108 insertions(+), 149 deletions(-) diff --git a/include/trackcpp/accelerator.h b/include/trackcpp/accelerator.h index 8a4922d..9ad28d3 100644 --- a/include/trackcpp/accelerator.h +++ b/include/trackcpp/accelerator.h @@ -39,8 +39,8 @@ class Accelerator { bool operator!=(const Accelerator& o) const { return !(*this == o); }; bool isequal(const Accelerator& a) const { return *this == a; } // necessary for python_package double get_length() const; - // double get_time_aware_frac() const; friend std::ostream& operator<< (std::ostream &out, const Accelerator& a); + double get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions, unsigned int element_offset = 0) const; }; #endif diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index 20f5873..0485046 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -29,12 +29,16 @@ Status::type track_findm66( Accelerator& accelerator, const Pos& fixed_point, - std::vector& tm, Matrix& m66, Pos& v0); + std::vector& tm, Matrix& m66, Pos& v0, double line_length, + std::vector time_aware_element_indices, + std::vector time_aware_element_positions); Status::type track_findm66( Accelerator& accelerator, const Pos& fixed_point, std::vector& tm, Matrix& m66, Pos& v0, - std::vector& indices); + std::vector& indices, double line_length, + std::vector time_aware_element_indices, + std::vector time_aware_element_positions); Status::type track_findorbit4( Accelerator& accelerator, std::vector >& closed_orbit, @@ -132,102 +136,6 @@ Status::type track_elementpass ( // pos: Pos vector of particles' final coordinates (or trajetory) // element_offset: in case of problems with passmethods, '*element_offset' is the index of the corresponding element // RETURN: status do tracking (see 'auxiliary.h') -// template -// Status::type track_linepass ( -// const Accelerator& accelerator, -// Pos& orig_pos, // initial electron coordinates -// std::vector >& pos, // vector with tracked electron coordinates at start of every element and at the end of last one. -// unsigned int& element_offset, // index of starting element for tracking -// Plane::type& lost_plane, // return plane in which particle was lost, if the case. -// std::vector& indices, double nturn = 0.0) {// indices to return; - -// Status::type status = Status::success; -// lost_plane = Plane::no_plane; - -// const std::vector& line = accelerator.lattice; -// int nr_elements = line.size(); - -// //pos.clear(); other functions assume pos is not clearedin linepass! -// pos.reserve(pos.size() + indices.size()); - -// // create vector of booleans help determine when to store position -// std::vector indcs; -// indcs.reserve(nr_elements+1); -// for (unsigned int i=0; i<=nr_elements; ++i) indcs[i] = false; -// for (auto&& i: indices) if (i<=nr_elements) indcs[i] = true; - -// for(int i=0; i= element.hmax))) { -// lost_plane = Plane::x; -// status = Status::particle_lost; -// } -// if (((ry <= element.vmin) or (ry >= element.vmax))) { -// if (status != Status::particle_lost) { -// lost_plane = Plane::y; -// status = Status::particle_lost; -// } else { -// lost_plane = Plane::xy; -// } -// } -// } else if (get_norm_amp_in_vchamber(element, rx, ry) > 1) { -// // lost in rhombus, elliptic and all finite p-norm shapes -// lost_plane = Plane::xy; -// status = Status::particle_lost; -// } -// } - -// if (status != Status::success) { -// // fill rest of vector with nans -// for(int ii=i+1; ii<=nr_elements; ++ii) { -// if (indcs[ii]) pos.emplace_back( -// nan(""),nan(""),nan(""),nan(""),nan(""),nan("")); -// } -// return status; -// } -// // moves to next element index -// element_offset = (element_offset + 1) % nr_elements; -// } - -// // stores final particle position at the end of the line -// if (indcs[nr_elements]) pos.push_back(orig_pos); - -// return (status == Status::success) ? status: Status::particle_lost; -// } - -// accelerator_length, ddl, TAW_positions template Status::type track_linepass ( const Accelerator& accelerator, @@ -236,7 +144,7 @@ Status::type track_linepass ( unsigned int& element_offset, // index of starting element for tracking Plane::type& lost_plane, // return plane in which particle was lost, if the case. std::vector& indices, // indices to return; - double line_length = 0.0, + double line_length = 1.0, std::vector time_aware_element_indices = {0, }, std::vector time_aware_element_positions = {0.0, }) { @@ -367,9 +275,9 @@ Status::type track_linepass ( unsigned int& element_offset, // index of starting element for tracking Plane::type& lost_plane, // return plane in which particle was lost, if the case. bool trajectory, // whether function should return coordinates at all elements - double line_length = 0.0, + double line_length = 1.0, std::vector time_aware_element_indices = {0,}, - std::vector time_aware_element_positions = {0,}) { + std::vector time_aware_element_positions = {0.0,}) { std::vector indices; unsigned int nr_elements = accelerator.lattice.size(); @@ -396,7 +304,7 @@ Status::type track_linepass ( std::vector& lost_plane, std::vector& lost_element, std::vector& indices, - double line_length = 0.0, + double line_length = 1.0, std::vector time_aware_element_indices = {0, }, std::vector time_aware_element_positions = {0.0, }) { @@ -443,8 +351,6 @@ Status::type track_linepass ( // turn_idx: in case of problems with passmethods, '*turn_idx' is the index of the corresponding turn // element_offset: in case of problems with passmethods, '*element_offset' is the index of the corresponding element // RETURN: status do tracking (see 'auxiliary.h') - - template Status::type track_ringpass ( const Accelerator& accelerator, @@ -460,27 +366,9 @@ Status::type track_ringpass ( std::vector > final_pos; // for longitudinal kick before RF cavities - std::vector TAW_positions = {0.0, }; - std::vector TAW_indices = {}; - unsigned int eoff = 0; eoff += element_offset; - size_t nr_elements = accelerator.lattice.size(); - PassMethodsClass PMClass = PassMethodsClass(); - double s_pos = 0.0; - for (size_t i = 0; i < nr_elements; i++) { - auto elem = accelerator.lattice[eoff]; - if (PMClass.is_time_aware_pm(elem.pass_method)) { - s_pos += (elem.length*0.5); - TAW_positions.push_back(s_pos); - TAW_indices.push_back(eoff); - s_pos += (elem.length*0.5); - } - else { - s_pos += elem.length; - } - eoff = (eoff + 1) % nr_elements; - } - TAW_positions.push_back(s_pos); - double accelerator_length = s_pos; + std::vector TAW_positions; + std::vector TAW_indices; + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, element_offset); if (trajectory) pos.reserve(nr_turns+1); diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 6c0d7ec..9a0e3c2 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -223,13 +223,17 @@ 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, + double line_length, + std::vector TAW_indices, + std::vector TAW_positions) { std::vector vec_tm; Matrix vec_m66; Status::type status = track_findm66( - accelerator, fixed_point, vec_tm, vec_m66, v0, indices); + accelerator, fixed_point, vec_tm, vec_m66, v0, indices, + line_length, TAW_indices, TAW_positions); for (unsigned int i=0; i& 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); + Pos& v0, std::vector& indices, + double line_length, + std::vector TAW_indices, + std::vector TAW_positions); Status::type track_diffusionmatrix_wrapper( diff --git a/src/accelerator.cpp b/src/accelerator.cpp index 85e1ebb..12083a8 100644 --- a/src/accelerator.cpp +++ b/src/accelerator.cpp @@ -61,3 +61,29 @@ std::ostream& operator<< (std::ostream &out, const Accelerator& a) { out << std::endl << "lattice_version: " << a.lattice_version; return out; } + +double Accelerator::get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions, unsigned int element_offset) const { + // for longitudinal kick before RF cavities + TAW_positions.clear(); + TAW_indices.clear(); + unsigned int eoff = 0; eoff += element_offset; + size_t nr_elements = this->lattice.size(); + PassMethodsClass PMClass = PassMethodsClass(); + double s_pos = 0.0; + TAW_positions.push_back(s_pos); + for (size_t i = 0; i < nr_elements; i++) { + auto elem = this->lattice[eoff]; + if (PMClass.is_time_aware_pm(elem.pass_method)) { + s_pos += (elem.length*0.5); + TAW_positions.push_back(s_pos); + TAW_indices.push_back(i); + s_pos += (elem.length*0.5); + } + else { + s_pos += elem.length; + } + eoff = (eoff + 1) % nr_elements; + } + TAW_positions.push_back(s_pos); + return s_pos; +} diff --git a/src/optics.cpp b/src/optics.cpp index 3b3cf56..2dab1e7 100644 --- a/src/optics.cpp +++ b/src/optics.cpp @@ -81,7 +81,14 @@ Status::type calc_twiss(Accelerator& accelerator, std::vector> closed_orbit; Plane::type lost_plane; unsigned int element_offset = 0; - Status::type status = track_linepass(accelerator, fp, closed_orbit, element_offset, lost_plane, true); + + // for longitudinal kick before RF cavities + std::vector TAW_positions; + std::vector TAW_indices; + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); + + + Status::type status = track_linepass(accelerator, fp, closed_orbit, element_offset, lost_plane, true, accelerator_length, TAW_indices, TAW_positions); if (status != Status::success) return status; @@ -96,7 +103,7 @@ Status::type calc_twiss(Accelerator& accelerator, std::vector atm; Pos v0; - status = track_findm66(accelerator, closed_orbit[0], atm, m66, v0); + status = track_findm66(accelerator, closed_orbit[0], atm, m66, v0, accelerator_length, TAW_indices, TAW_positions); if (status != Status::success) return status; #ifdef TIMEIT @@ -148,7 +155,7 @@ Status::type calc_twiss(Accelerator& accelerator, fpp.px += twiss0.etax[1] * dpp; fpp.ry += twiss0.etay[0] * dpp; fpp.py += twiss0.etay[1] * dpp; - Status::type status = track_linepass(accelerator, fpp, codp, element_offset, lost_plane, true); + Status::type status = track_linepass(accelerator, fpp, codp, element_offset, lost_plane, true, accelerator_length, TAW_indices, TAW_positions); if (status != Status::success) return status; } diff --git a/src/tracking.cpp b/src/tracking.cpp index 0e20b8d..7bfb96e 100644 --- a/src/tracking.cpp +++ b/src/tracking.cpp @@ -41,7 +41,10 @@ Status::type track_findm66 (Accelerator& accelerator, std::vector& tm, Matrix& m66, Pos& v0, - std::vector& indices) { + std::vector& indices, + double line_length = 1.0, + std::vector time_aware_element_indices = {0,}, + std::vector time_aware_element_positions = {0.0,}) { Status::type status = Status::success; const std::vector& lattice = accelerator.lattice; @@ -61,12 +64,16 @@ Status::type track_findm66 (Accelerator& accelerator, for (unsigned int i=0; i<=nr_elements; ++i) indcs[i] = false; for (auto&& i: indices) if (i<=nr_elements) indcs[i] = true; + Pos > map; map.rx = Tpsa<6,1>(fp.rx, 0); map.px = Tpsa<6,1>(fp.px, 1); map.ry = Tpsa<6,1>(fp.ry, 2); map.py = Tpsa<6,1>(fp.py, 3); map.de = Tpsa<6,1>(fp.de, 4); map.dl = Tpsa<6,1>(fp.dl, 5); tm.clear(); tm.reserve(indices.size()); + + unsigned int TAW_pivot = 0; + double ddl = 0; for(unsigned int i=0; i& fixed_point, std::vector& tm, Matrix& m66, - Pos& v0) { + Pos& v0, + double line_length = 1.0, + std::vector time_aware_element_indices = {0,}, + std::vector time_aware_element_positions = {0.0,}) { std::vector indices; unsigned int nr_elements = accelerator.lattice.size(); @@ -127,7 +146,7 @@ Status::type track_findm66 (Accelerator& accelerator, indices.reserve(nr_elements + 1); for (unsigned int i=0; i<=nr_elements; ++i) indices.push_back(i); - return track_findm66 (accelerator, fixed_point, tm, m66, v0, indices); + return track_findm66 (accelerator, fixed_point, tm, m66, v0, indices, line_length, time_aware_element_indices, time_aware_element_positions); } @@ -147,6 +166,12 @@ Status::type track_findorbit6( accelerator.radiation_on = RadiationState::damping; } + // for longitudinal kick before RF cavities + std::vector TAW_positions; + std::vector TAW_indices; + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); + + // temporary vectors and matrices std::vector > co(7,0); for(auto i=0; i<7; ++i) co[i] = fixed_point_guess; @@ -165,13 +190,13 @@ Status::type track_findorbit6( Plane::type lost_plane; Status::type status = Status::success; - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[0], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[1], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[2], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[3], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[4], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[5], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[6], co2, element_offset, lost_plane, false)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[0], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[1], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[2], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[3], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[4], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[5], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[6], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); if (status != Status::success) { return Status::findorbit_one_turn_matrix_problem; @@ -206,7 +231,7 @@ Status::type track_findorbit6( closed_orbit.clear(); unsigned int element_offset = 0; Plane::type lost_plane; - track_linepass(accelerator, co[6], closed_orbit, element_offset, lost_plane, true); + track_linepass(accelerator, co[6], closed_orbit, element_offset, lost_plane, true, accelerator_length, TAW_indices, TAW_positions); accelerator.radiation_on = radsts; return Status::success; @@ -228,6 +253,12 @@ Status::type track_findorbit4( if (radsts == RadiationState::full){ accelerator.radiation_on = RadiationState::damping; } + + // for longitudinal kick before RF cavities + std::vector TAW_positions; + std::vector TAW_indices; + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); + // temporary vectors and matrices // std::vector > co(7,0); // no initial guess std::vector > co(7,fixed_point_guess); @@ -246,11 +277,11 @@ Status::type track_findorbit4( unsigned int element_offset = 0; Plane::type lost_plane; Status::type status = Status::success; - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[0], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[1], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[2], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[3], co2, element_offset, lost_plane, false)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[6], co2, element_offset, lost_plane, false)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[0], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[1], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[2], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[3], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[6], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); if (status != Status::success) { return Status::findorbit_one_turn_matrix_problem; } @@ -279,7 +310,7 @@ Status::type track_findorbit4( closed_orbit.clear(); unsigned int element_offset = 0; Plane::type lost_plane; - track_linepass(accelerator, co[6], closed_orbit, element_offset, lost_plane, true); + track_linepass(accelerator, co[6], closed_orbit, element_offset, lost_plane, true, accelerator_length, TAW_indices, TAW_positions); accelerator.radiation_on = radsts; return Status::success; From a8d2b72be0b1c1d655ae6a54ef675253a1165c38 Mon Sep 17 00:00:00 2001 From: vitor Date: Tue, 11 Jun 2024 12:25:45 -0300 Subject: [PATCH 10/20] restoring spaces and deleting comments --- include/trackcpp/auxiliary.h | 4 ++-- include/trackcpp/passmethods.h | 2 +- include/trackcpp/passmethods.hpp | 28 ---------------------------- include/trackcpp/tracking.h | 1 + src/accelerator.cpp | 9 --------- src/kicktable.cpp | 1 - src/optics.cpp | 2 -- src/tracking.cpp | 9 ++------- 8 files changed, 6 insertions(+), 50 deletions(-) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 8058f5e..6e368c4 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -31,7 +31,7 @@ class PassMethodsClass { static const int pm_str_mpole_symplectic4_pass = 2; static const int pm_bnd_mpole_symplectic4_pass = 3; static const int pm_corrector_pass = 4; - static const int pm_cavity_pass = 5; + static const int pm_cavity_pass = 5; static const int pm_thinquad_pass = 6; static const int pm_thinsext_pass = 7; static const int pm_kickmap_pass = 8; @@ -74,7 +74,7 @@ struct PassMethod { pm_str_mpole_symplectic4_pass = 2, pm_bnd_mpole_symplectic4_pass = 3, pm_corrector_pass = 4, - pm_cavity_pass = 5, + pm_cavity_pass = 5, pm_thinquad_pass = 6, pm_thinsext_pass = 7, pm_kickmap_pass = 8, diff --git a/include/trackcpp/passmethods.h b/include/trackcpp/passmethods.h index 2ab02a3..51a38aa 100644 --- a/include/trackcpp/passmethods.h +++ b/include/trackcpp/passmethods.h @@ -60,7 +60,7 @@ template Status::type pm_drift_pass (Pos &pos, c template Status::type pm_str_mpole_symplectic4_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_bnd_mpole_symplectic4_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_corrector_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); -template Status::type pm_cavity_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); +template Status::type pm_cavity_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_thinquad_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_thinsext_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); template Status::type pm_kickmap_pass (Pos &pos, const Element &elem, const Accelerator& accelerator); diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 0578331..64f7a97 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -392,7 +392,6 @@ Status::type pm_corrector_pass(Pos &pos, const Element &elem, return Status::success; } - template Status::type pm_cavity_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { @@ -402,31 +401,9 @@ Status::type pm_cavity_pass(Pos &pos, const Element &elem, global_2_local(pos, elem); double nv = elem.voltage / accelerator.energy; double frf = elem.frequency; - // double harmonic_number = accelerator.harmonic_number; - // double L0 = accelerator.get_length(); - // double s = 0.0; - // double accum_s = 0.0; - // unsigned int last_cav_idx = 0; - // for (unsigned int i = 0; i < accelerator.lattice.size(); i++){ - // auto elem2 = accelerator.lattice[i]; - // if (elem2 == elem){ - // accum_s += s; - // } - // if (elem2.frequency != 0.0){ - // s = 0.0; - // last_cav_idx = i; - // } - // s += elem2.length; - // } - // double ddl = (light_speed*harmonic_number/frf - L0); - // std::cout << "accum_s = " << accum_s << ", s = " << s << ", last_cav_index: " << last_cav_idx << std::endl; if (elem.length == 0) { T &de = pos.de, &dl = pos.dl; - // dl -= ddl*(accum_s/L0); de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); - // if (elem == accelerator.lattice[last_cav_idx]){ - // dl -= ddl*(s/L0); - // } } else { T &rx = pos.rx, &px = pos.px; T &ry = pos.ry, &py = pos.py; @@ -438,11 +415,7 @@ Status::type pm_cavity_pass(Pos &pos, const Element &elem, ry += norml * py; dl += 0.5 * norml * pnorm * (px*px + py*py); // longitudinal momentum kick - // dl -= ddl*(accum_s/L0); de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); - // if (elem == accelerator.lattice[last_cav_idx]){ - // dl -= ddl*(s/L0); - // } // drift half length pnorm = 1.0 / (1.0 + de); norml = (0.5 * elem.length) * pnorm; @@ -468,7 +441,6 @@ Status::type pm_thinsext_pass(Pos &pos, const Element &elem, return Status::passmethod_not_implemented; } - template Status::type pm_kickmap_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index 0485046..f481910 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -136,6 +136,7 @@ Status::type track_elementpass ( // pos: Pos vector of particles' final coordinates (or trajetory) // element_offset: in case of problems with passmethods, '*element_offset' is the index of the corresponding element // RETURN: status do tracking (see 'auxiliary.h') + template Status::type track_linepass ( const Accelerator& accelerator, diff --git a/src/accelerator.cpp b/src/accelerator.cpp index 12083a8..676dfcc 100644 --- a/src/accelerator.cpp +++ b/src/accelerator.cpp @@ -27,15 +27,6 @@ double Accelerator::get_length() const { return length; } -// double Accelerator::get_time_aware_frac() const { -// double count = 0.0; -// for(auto i=0; ienergy != o.energy) return false; diff --git a/src/kicktable.cpp b/src/kicktable.cpp index 0034cba..4495bfc 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -104,7 +104,6 @@ Status::type Kicktable::load_from_file(const std::string& filename_) { } int add_kicktable(const std::string& filename) { - // looks through vector of kicktables... for(unsigned int i=0; i TAW_indices; double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); - Status::type status = track_linepass(accelerator, fp, closed_orbit, element_offset, lost_plane, true, accelerator_length, TAW_indices, TAW_positions); if (status != Status::success) return status; - #ifdef TIMEIT end = std::chrono::steady_clock::now(); diff = end - start; std::cout << "calc_twiss::close_orbit: " << std::chrono::duration (diff).count() << " ms" << std::endl; diff --git a/src/tracking.cpp b/src/tracking.cpp index 7bfb96e..32ddc3f 100644 --- a/src/tracking.cpp +++ b/src/tracking.cpp @@ -43,12 +43,11 @@ Status::type track_findm66 (Accelerator& accelerator, Pos& v0, std::vector& indices, double line_length = 1.0, - std::vector time_aware_element_indices = {0,}, - std::vector time_aware_element_positions = {0.0,}) { + std::vector time_aware_element_indices = {0,}, + std::vector time_aware_element_positions = {0.0,}) { Status::type status = Status::success; const std::vector& lattice = accelerator.lattice; - Pos fp = fixed_point; const int radsts = accelerator.radiation_on; @@ -64,14 +63,12 @@ Status::type track_findm66 (Accelerator& accelerator, for (unsigned int i=0; i<=nr_elements; ++i) indcs[i] = false; for (auto&& i: indices) if (i<=nr_elements) indcs[i] = true; - Pos > map; map.rx = Tpsa<6,1>(fp.rx, 0); map.px = Tpsa<6,1>(fp.px, 1); map.ry = Tpsa<6,1>(fp.ry, 2); map.py = Tpsa<6,1>(fp.py, 3); map.de = Tpsa<6,1>(fp.de, 4); map.dl = Tpsa<6,1>(fp.dl, 5); tm.clear(); tm.reserve(indices.size()); - unsigned int TAW_pivot = 0; double ddl = 0; for(unsigned int i=0; i Date: Tue, 11 Jun 2024 14:34:09 -0300 Subject: [PATCH 11/20] saving corrections --- include/trackcpp/accelerator.h | 2 +- include/trackcpp/tracking.h | 46 +++++++++++++++++++++------------- src/accelerator.cpp | 6 ++--- src/commands.cpp | 4 +-- src/tests.cpp | 6 ++--- src/tracking.cpp | 12 ++++----- 6 files changed, 42 insertions(+), 34 deletions(-) diff --git a/include/trackcpp/accelerator.h b/include/trackcpp/accelerator.h index 9ad28d3..5ce1a2a 100644 --- a/include/trackcpp/accelerator.h +++ b/include/trackcpp/accelerator.h @@ -40,7 +40,7 @@ class Accelerator { bool isequal(const Accelerator& a) const { return *this == a; } // necessary for python_package double get_length() const; friend std::ostream& operator<< (std::ostream &out, const Accelerator& a); - double get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions, unsigned int element_offset = 0) const; + double get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions) const; }; #endif diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index f481910..38ada43 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -29,16 +29,18 @@ Status::type track_findm66( Accelerator& accelerator, const Pos& fixed_point, - std::vector& tm, Matrix& m66, Pos& v0, double line_length, - std::vector time_aware_element_indices, - std::vector time_aware_element_positions); + std::vector& tm, Matrix& m66, Pos& v0, + const double line_length, + const std::vector& time_aware_element_indices, + const std::vector& time_aware_element_positions); Status::type track_findm66( Accelerator& accelerator, const Pos& fixed_point, std::vector& tm, Matrix& m66, Pos& v0, - std::vector& indices, double line_length, - std::vector time_aware_element_indices, - std::vector time_aware_element_positions); + std::vector& indices, + const double line_length, + const std::vector& time_aware_element_indices, + const std::vector& time_aware_element_positions); Status::type track_findorbit4( Accelerator& accelerator, std::vector >& closed_orbit, @@ -145,9 +147,9 @@ Status::type track_linepass ( unsigned int& element_offset, // index of starting element for tracking Plane::type& lost_plane, // return plane in which particle was lost, if the case. std::vector& indices, // indices to return; - double line_length = 1.0, - std::vector time_aware_element_indices = {0, }, - std::vector time_aware_element_positions = {0.0, }) { + const double line_length, + const std::vector& time_aware_element_indices, + const std::vector& time_aware_element_positions) { Status::type status = Status::success; lost_plane = Plane::no_plane; @@ -165,6 +167,12 @@ Status::type track_linepass ( for (unsigned int i=0; i<=nr_elements; ++i) indcs[i] = false; for (auto&& i: indices) if (i<=nr_elements) indcs[i] = true; + // std::cout << "len acc = " << line_length << std::endl; + // for (size_t i; i& args) { std::vector> pos_list; Plane::type lost_plane; unsigned int offset_element = start_element; - track_linepass(accelerator, pos, pos_list, offset_element, lost_plane, true); + track_linepass(accelerator, pos, pos_list, offset_element, lost_plane, true, 0, {0,}, {0.0, 0.0, }); std::cout << get_timestamp() << " saving track_linepass data to file" << std::endl; status = print_tracking_linepass(accelerator, pos_list, start_element, "track_linepass_out.txt"); diff --git a/src/tests.cpp b/src/tests.cpp index 0696a5b..d30a451 100644 --- a/src/tests.cpp +++ b/src/tests.cpp @@ -36,7 +36,7 @@ int test_linepass(const Accelerator& accelerator) { std::vector > new_pos; unsigned int element_offset = 0; Plane::type lost_plane; - Status::type status = track_linepass(accelerator, pos, new_pos, element_offset, lost_plane, true); + Status::type status = track_linepass(accelerator, pos, new_pos, element_offset, lost_plane, true, 0, {0,}, {0.0, 0.0, }); std::cout << "status: " << string_error_messages[status] << std::endl; @@ -65,7 +65,7 @@ int test_linepass_tpsa(const Accelerator& accelerator, const std::vector > > new_tpsa; unsigned int element_offset = 0; Plane::type lost_plane; - track_linepass(accelerator, tpsa, new_tpsa, element_offset, lost_plane, false); + track_linepass(accelerator, tpsa, new_tpsa, element_offset, lost_plane, false, 0, {0,}, {0.0, 0.0, }); for(unsigned int i=0; i >& c = new_particles[i]; //std::cout << c.rx << std::endl; @@ -403,7 +403,7 @@ int test_linepass2() { pos, // vector with electron coordinates from tracking at every element. element_offset, // index of starting element for tracking lost_plane, // return plane in which particle was lost, if the case. - trajectory); + trajectory, 0, {0,}, {0.0, 0.0, }); for(unsigned int i=0; i<10; ++i) { std::cout << std::endl; diff --git a/src/tracking.cpp b/src/tracking.cpp index 32ddc3f..446be00 100644 --- a/src/tracking.cpp +++ b/src/tracking.cpp @@ -42,9 +42,9 @@ Status::type track_findm66 (Accelerator& accelerator, Matrix& m66, Pos& v0, std::vector& indices, - double line_length = 1.0, - std::vector time_aware_element_indices = {0,}, - std::vector time_aware_element_positions = {0.0,}) { + const double line_length, + const std::vector& time_aware_element_indices, + const std::vector& time_aware_element_positions) { Status::type status = Status::success; const std::vector& lattice = accelerator.lattice; @@ -131,9 +131,9 @@ Status::type track_findm66 (Accelerator& accelerator, std::vector& tm, Matrix& m66, Pos& v0, - double line_length = 1.0, - std::vector time_aware_element_indices = {0,}, - std::vector time_aware_element_positions = {0.0,}) { + const double line_length, + const std::vector& time_aware_element_indices, + const std::vector& time_aware_element_positions) { std::vector indices; unsigned int nr_elements = accelerator.lattice.size(); From e174bf287edfcc415241ec6bddb3553086eec9bb Mon Sep 17 00:00:00 2001 From: vitor Date: Sat, 31 Aug 2024 18:34:10 -0300 Subject: [PATCH 12/20] merge master (manually) --- include/trackcpp/tracking.h | 415 ++++++++++++++++++++++------------- python_package/interface.cpp | 80 +++++-- python_package/interface.h | 52 +++-- python_package/trackcpp.i | 2 + src/commands.cpp | 2 +- src/diffusion_matrix.cpp | 2 +- src/optics.cpp | 5 +- src/tests.cpp | 37 +++- src/tracking.cpp | 58 +++-- tests/test-kickmap.cpp | 16 +- 10 files changed, 440 insertions(+), 229 deletions(-) diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index 38ada43..40ce76c 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -59,9 +59,10 @@ Pos linalg_solve6_posvec( template Status::type track_elementpass ( - const Element& el, // element through which to track particle - Pos &orig_pos, // initial electron coordinates - const Accelerator& accelerator) { + const Accelerator& accelerator, + const Element& el, // element through which to track particle + Pos &orig_pos // initial electron coordinates +) { Status::type status = Status::success; @@ -109,12 +110,12 @@ Status::type track_elementpass ( template Status::type track_elementpass ( - const Element& el, // element through which to track particle - std::vector >& orig_pos, // initial electron coordinates - const Accelerator& accelerator) { + const Accelerator& accelerator, + const Element& el, // element through which to track particle + std::vector >& orig_pos // initial electron coordinates +) { Status::type status = Status::success; - for(auto&& pos: orig_pos) { Status::type status2 = track_elementpass(el, pos, accelerator); if (status2 != Status::success) status = status2; @@ -129,34 +130,40 @@ Status::type track_elementpass ( // tracks particles along a beam transport line // // inputs: -// line: Element vector representing the beam transport line -// orig_pos: Pos vector representing initial positions of particles -// element_offset: equivalent to shifting the lattice so that '*element_offset' is the index for the first element -// trajectory: flag indicating that trajectory is to be recorded at entrance of all elements -// (otherwise only the coordinates at the exit of last element is recorded) +// accelerator: Element vector representing the beam transport line. +// orig_pos: Pos vector representing initial positions of particles. +// indices: vector of integers defining at which elements to return. +// element_offset: equivalent to shifting the lattice so that +// '*element_offset' is the index for the first element. +// // outputs: -// pos: Pos vector of particles' final coordinates (or trajetory) -// element_offset: in case of problems with passmethods, '*element_offset' is the index of the corresponding element -// RETURN: status do tracking (see 'auxiliary.h') - +// element_offset: in case of problems with passmethods, such as particle +// loss '*element_offset' is the index of the corresponding element. +// pos: Pos vector of particles coordinates at the start of every element +// selected by indices and at the end of the last element. +// lost_plane: which plane the particle was lost. +// lost_pos: return coordinates of lost particle. +// +// RETURN: +// status do tracking (see 'auxiliary.h') +// template Status::type track_linepass ( const Accelerator& accelerator, - Pos& orig_pos, // initial electron coordinates - std::vector >& pos, // vector with tracked electron coordinates at start of every element and at the end of last one. - unsigned int& element_offset, // index of starting element for tracking - Plane::type& lost_plane, // return plane in which particle was lost, if the case. - std::vector& indices, // indices to return; + Pos& orig_pos, + const std::vector& indices, + unsigned int& element_offset, + std::vector >& pos, + Plane::type& lost_plane, const double line_length, const std::vector& time_aware_element_indices, - const std::vector& time_aware_element_positions) { + const std::vector& time_aware_element_positions +) { Status::type status = Status::success; - lost_plane = Plane::no_plane; - double ddl = 0.0; - const std::vector& line = accelerator.lattice; int nr_elements = line.size(); + double ddl = 0.0; //pos.clear(); other functions assume pos is not clearedin linepass! pos.reserve(pos.size() + indices.size()); @@ -189,55 +196,15 @@ Status::type track_linepass ( TAW_pivot++; } - status = track_elementpass (element, orig_pos, accelerator); + status = track_elementpass(accelerator, element, orig_pos); if (i == time_aware_element_indices.back()) { orig_pos.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; // std::cout << "taw pivot = " << TAW_pivot << std::endl; } - const T& rx = orig_pos.rx; - const T& ry = orig_pos.ry; - - // checks if particle is lost - if (not isfinite(rx)) { - lost_plane = Plane::x; - status = Status::particle_lost; - } - if (not isfinite(ry)) { - if (status != Status::particle_lost) { - lost_plane = Plane::y; - status = Status::particle_lost; - } else { - lost_plane = Plane::xy; - } - } - if ((status != Status::particle_lost) and accelerator.vchamber_on) { - if (element.vchamber < 0) { - // invalid p-norm shape (negative p) - // safely signals lost particle - lost_plane = Plane::xy; - status = Status::particle_lost; - } else if (element.vchamber == VChamberShape::rectangle) { - // rectangular vacuum chamber - if (((rx <= element.hmin) or (rx >= element.hmax))) { - lost_plane = Plane::x; - status = Status::particle_lost; - } - if (((ry <= element.vmin) or (ry >= element.vmax))) { - if (status != Status::particle_lost) { - lost_plane = Plane::y; - status = Status::particle_lost; - } else { - lost_plane = Plane::xy; - } - } - } else if (get_norm_amp_in_vchamber(element, rx, ry) > 1) { - // lost in rhombus, elliptic and all finite p-norm shapes - lost_plane = Plane::xy; - status = Status::particle_lost; - } - } + lost_plane = check_particle_loss(accelerator, element, orig_pos); + if (lost_plane != Plane::no_plane) status = Status::particle_lost; if (status != Status::success) { // fill rest of vector with nans @@ -257,42 +224,42 @@ Status::type track_linepass ( return (status == Status::success) ? status: Status::particle_lost; } - -template -T get_norm_amp_in_vchamber(const Element& elem, const T& rx, const T& ry) { - double lx = (elem.hmax - elem.hmin) / 2; - double ly = (elem.vmax - elem.vmin) / 2; - double xc = (elem.hmax + elem.hmin) / 2; - double yc = (elem.vmax + elem.vmin) / 2; - T xn = abs((rx - xc)/lx); - T yn = abs((ry - yc)/ly); - T amplitude; - if (elem.vchamber == VChamberShape::rhombus) { - amplitude = xn + yn; - } else if (elem.vchamber == VChamberShape::ellipse) { - amplitude = xn*xn + yn*yn; - } else { - amplitude = pow(xn, elem.vchamber) + pow(yn, elem.vchamber); - } - return amplitude; -} - - +// linepass +// -------- +// tracks particles along a beam transport line +// +// inputs: +// accelerator: Element vector representing the beam transport line. +// orig_pos: Pos vector representing initial positions of particles. +// trajectory: whether to return coordinates at all elements. +// element_offset: equivalent to shifting the lattice so that +// '*element_offset' is the index for the first element. +// +// outputs: +// element_offset: in case of problems with passmethods, such as particle +// loss '*element_offset' is the index of the corresponding element. +// pos: Pos vector of particles coordinates at the start of every element, +// if trajectory == true, and at the end of the last element. +// lost_plane: which plane the particle was lost. +// lost_pos: return coordinates of lost particle. +// +// RETURN: +// status do tracking (see 'auxiliary.h') +// template Status::type track_linepass ( const Accelerator& accelerator, - Pos& orig_pos, // initial electron coordinates - std::vector >& pos, // vector with tracked electron coordinates at start of every element and at the end of last one. - unsigned int& element_offset, // index of starting element for tracking - Plane::type& lost_plane, // return plane in which particle was lost, if the case. - bool trajectory, // whether function should return coordinates at all elements + Pos& orig_pos, + const bool trajectory, + unsigned int& element_offset, + std::vector >& pos, + Plane::type& lost_plane, const double line_length, const std::vector& time_aware_element_indices, - const std::vector& time_aware_element_positions) { - + const std::vector& time_aware_element_positions +) { std::vector indices; unsigned int nr_elements = accelerator.lattice.size(); - if (trajectory){ indices.reserve(nr_elements + 1); for (unsigned int i=0; i<=nr_elements; ++i) indices.push_back(i); @@ -301,23 +268,56 @@ Status::type track_linepass ( } return track_linepass ( - accelerator, orig_pos, pos, element_offset, lost_plane, indices, - line_length, time_aware_element_indices, time_aware_element_positions); + accelerator, + orig_pos, + indices, + element_offset, + pos, + lost_plane, + line_length, + time_aware_element_indices, + time_aware_element_positions + ); } +// linepass +// -------- +// tracks particles along a beam transport line +// +// inputs: +// accelerator: Element vector representing the beam transport line. +// orig_pos: Pos vector representing initial positions of particles. +// indices: vector of integers defining at which elements to return. +// element_offset: equivalent to shifting the lattice so that +// '*element_offset' is the index for the first element. +// +// outputs: +// element_offset: in case of problems with passmethods, such as particle +// loss '*element_offset' is the index of the corresponding element. +// pos: Pos vector of particles coordinates at the start of every element +// selected by indices and at the end of the last element. +// lost_plane: which plane each particle was lost. +// lost_flag: whether or not each particle was lost. +// lost_element: which element each particle was lost. +// +// RETURN: +// status do tracking (see 'auxiliary.h') +// template Status::type track_linepass ( const Accelerator& accelerator, - std::vector > &orig_pos, - std::vector > &pos, - unsigned int element_offset, + std::vector> &orig_pos, + const std::vector& indices, + const unsigned int element_offset, + std::vector> &pos, std::vector& lost_plane, - std::vector& lost_element, - std::vector& indices, + std::vector& lost_flag, + std::vector& lost_element, const double line_length, const std::vector& time_aware_element_indices, - const std::vector& time_aware_element_positions) { + const std::vector& time_aware_element_positions +) { int nr_elements = accelerator.lattice.size(); Status::type status = Status::success; @@ -326,6 +326,7 @@ Status::type track_linepass ( Plane::type lp; pos.reserve(indices.size() * orig_pos.size()); + lost_flag.reserve(orig_pos.size()); lost_plane.reserve(orig_pos.size()); lost_element.reserve(orig_pos.size()); @@ -333,13 +334,17 @@ Status::type track_linepass ( unsigned int le = element_offset; status2 = track_linepass ( - accelerator, orig_pos[i], final_pos, le, lp, indices, - line_length, time_aware_element_indices, time_aware_element_positions); - - if (status2 != Status::success) status = status2; - + accelerator, orig_pos[i], indices, le, final_pos, lp,line_length, time_aware_element_indices, time_aware_element_positions); + + if (status2 != Status::success){ + status = status2; + lost_element.push_back(le); + lost_flag.push_back(true); + } else { + lost_element.push_back(-1); + lost_flag.push_back(false); + } lost_plane.push_back(lp); - lost_element.push_back(le); for (auto&& p: final_pos) pos.push_back(p); final_pos.clear(); } @@ -347,31 +352,63 @@ Status::type track_linepass ( } +// linepass +// -------- +// tracks particles along a beam transport line +// +// inputs: +// accelerator: Element vector representing the beam transport line. +// orig_pos: Pos vector representing initial positions of particles. +// indices: vector of integers defining at which elements to return. +// element_offset: equivalent to shifting the lattice so that +// '*element_offset' is the index for the first element. +// +// outputs: +// element_offset: in case of problems with passmethods, such as particle +// loss '*element_offset' is the index of the corresponding element. +// pos: Pos vector of particles coordinates at the start of every element, +// if trajectory == true, and at the end of the last element. +// lost_plane: which plane each particle was lost. +// lost_flag: whether or not each particle was lost. +// lost_element: which element each particle was lost. +// +// RETURN: +// status do tracking (see 'auxiliary.h') +// + // ringpass // -------- // tracks particles around a ring // // inputs: -// ring: Element vector representing the ring -// orig_pos: Pos vector representing initial positions of particles -// element_offset: equivalent to shifting the lattice so that '*element_offset' is the index for the first element -// nr_turns: number of turns for tracking +// accelerator: Element vector representing the beam transport line. +// orig_pos: Pos vector representing initial positions of particles. +// nr_turns: number of turns for tracking +// turn_by_turn: whether to return coordinates at all elements. +// element_offset: equivalent to shifting the lattice so that +// '*element_offset' is the index for the first element // // outputs: -// pos: Pos vector of particles' coordinates of the turn_by_turn data (at end of the ring at each turn) -// turn_idx: in case of problems with passmethods, '*turn_idx' is the index of the corresponding turn -// element_offset: in case of problems with passmethods, '*element_offset' is the index of the corresponding element -// RETURN: status do tracking (see 'auxiliary.h') +// element_offset: in case of problems with passmethods, such as particle +// loss '*element_offset' is the index of the corresponding element. +// pos: Pos vector of particles coordinates of the turn_by_turn data +// (at end of the ring at each turn) +// lost_plane: which plane each particle was lost. +// lost_turn: which turn each particle was lost. +// +// RETURN: +// status do tracking (see 'auxiliary.h') template Status::type track_ringpass ( - const Accelerator& accelerator, - Pos &orig_pos, - std::vector >& pos, + const Accelerator& accelerator, + Pos& orig_pos, const unsigned int nr_turns, - unsigned int& lost_turn, + const bool turn_by_turn, unsigned int& element_offset, + std::vector >& pos, Plane::type& lost_plane, - bool trajectory) { + unsigned int& lost_turn +) { Status::type status = Status::success; std::vector > final_pos; @@ -381,21 +418,29 @@ Status::type track_ringpass ( std::vector TAW_indices; double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); - if (trajectory) pos.reserve(nr_turns+1); + if (turn_by_turn) pos.reserve(nr_turns+1); for(lost_turn=0; lost_turn Status::type track_ringpass ( - const Accelerator& accelerator, - std::vector > &orig_pos, - std::vector > &pos, - const unsigned int nr_turns, - std::vector &lost_turn, - unsigned int element_offset, - std::vector& lost_plane, - std::vector& lost_element, - bool trajectory) { + const Accelerator& accelerator, + std::vector > &orig_pos, + const unsigned int nr_turns, + bool turn_by_turn, + unsigned int element_offset, + std::vector > &pos, + std::vector& lost_plane, + std::vector& lost_turn, + std::vector& lost_flag, + std::vector& lost_element +) { Status::type status = Status::success; Status::type status2 = Status::success; @@ -425,28 +472,102 @@ Status::type track_ringpass ( unsigned int lt; Plane::type lp; - if (trajectory) pos.reserve((nr_turns + 1) * orig_pos.size()); + if (turn_by_turn) pos.reserve((nr_turns + 1) * orig_pos.size()); else pos.reserve(orig_pos.size()); lost_turn.reserve(orig_pos.size()); lost_plane.reserve(orig_pos.size()); lost_element.reserve(orig_pos.size()); + lost_flag.reserve(orig_pos.size()); for(unsigned int i=0; i +Plane::type check_particle_loss( + const Accelerator& accelerator, + const Element& ele, + const Pos& orig_pos +){ + + Plane::type lost_plane = Plane::no_plane; + const T& rx = orig_pos.rx; + const T& ry = orig_pos.ry; + + if (not isfinite(rx)) lost_plane = Plane::x; + if (not isfinite(ry)) { + if (lost_plane == Plane::no_plane) return Plane::y; + else return Plane::xy; + } + if (lost_plane != Plane::no_plane) return lost_plane; + + if (not accelerator.vchamber_on) return Plane::no_plane; + + // invalid p-norm shape (negative p). Safely signals lost particle + if (ele.vchamber < 0) return Plane::xy; + // rectangular vacuum chamber + if (ele.vchamber == VChamberShape::rectangle) { + if (((rx <= ele.hmin) or (rx >= ele.hmax))) lost_plane = Plane::x; + if (((ry <= ele.vmin) or (ry >= ele.vmax))) { + if (lost_plane == Plane::no_plane) return Plane::y; + else return Plane::xy; + } + if (lost_plane != Plane::no_plane) return lost_plane; + } + // lost in rhombus, elliptic and all finite p-norm shapes + else if (get_norm_amp_in_vchamber(ele, rx, ry) > 1) return Plane::xy; + return lost_plane; +} + + +template +T get_norm_amp_in_vchamber(const Element& elem, const T& rx, const T& ry) { + double lx = (elem.hmax - elem.hmin) / 2; + double ly = (elem.vmax - elem.vmin) / 2; + double xc = (elem.hmax + elem.hmin) / 2; + double yc = (elem.vmax + elem.vmin) / 2; + T xn = abs((rx - xc)/lx); + T yn = abs((ry - yc)/ly); + T amplitude; + if (elem.vchamber == VChamberShape::rhombus) { + amplitude = xn + yn; + } else if (elem.vchamber == VChamberShape::ellipse) { + amplitude = xn*xn + yn*yn; + } else { + amplitude = pow(xn, elem.vchamber) + pow(yn, elem.vchamber); + } + return amplitude; +} + #endif diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 9a0e3c2..51906df 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -19,9 +19,10 @@ Status::type track_elementpass_wrapper ( - const Element& el, - double *pos, int n1, int n2, - const Accelerator& accelerator) { + const Accelerator& accelerator, + const Element& el, + double *pos, int n1, int n2 +) { std::vector > post; @@ -32,8 +33,7 @@ Status::type track_elementpass_wrapper ( pos[2*n2 + i], pos[3*n2 + i], pos[4*n2 + i], pos[5*n2 + i]); } - Status::type status = track_elementpass( - el, post, accelerator); + Status::type status = track_elementpass(accelerator, el, post); for (unsigned int i=0; i indices; - std::vector< unsigned int > lost_plane; - std::vector< unsigned int > lost_element; + std::vector indices; + std::vector lost_plane; + std::vector lost_element; + std::vector lost_flag; double line_length; std::vector time_aware_element_indices; std::vector time_aware_element_positions; @@ -39,11 +40,12 @@ struct LinePassArgs { struct RingPassArgs : public LinePassArgs { unsigned int nr_turns; - std::vector< unsigned int > lost_turn; + unsigned int turn_by_turn; unsigned int element_offset; - std::vector< unsigned int > lost_plane; - std::vector< unsigned int > lost_element; - unsigned int trajectory; + std::vector lost_plane; + std::vector lost_turn; + std::vector lost_element; + std::vector lost_flag; }; struct String { @@ -53,28 +55,32 @@ struct String { }; Status::type track_elementpass_wrapper ( - const Element& el, - double *pos, int n1, int n2, - const Accelerator& accelerator); + const Accelerator& accelerator, + const Element& el, + double *pos, int n1, int n2 +); Status::type track_linepass_wrapper ( - const Accelerator& accelerator, - double *orig_pos, int ni1, int ni2, - double *pos, int n1, int n2, - LinePassArgs& args); + const Accelerator& accelerator, + double *orig_pos, int ni1, int ni2, + double *pos, int n1, int n2, + LinePassArgs& args +); Status::type track_ringpass_wrapper ( - const Accelerator& accelerator, - double *orig_pos, int ni1, int ni2, - double *pos, int n1, int n2, - RingPassArgs& args); + const Accelerator& accelerator, + double *orig_pos, int ni1, int ni2, + double *pos, int n1, int n2, + RingPassArgs& args +); 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 +); Element marker_wrapper(const std::string& fam_name_); Element bpm_wrapper(const std::string& fam_name_); diff --git a/python_package/trackcpp.i b/python_package/trackcpp.i index 240627a..766a189 100644 --- a/python_package/trackcpp.i +++ b/python_package/trackcpp.i @@ -42,6 +42,8 @@ namespace std { %template(CppStringVector) vector; %template(CppDoubleVector) vector; %template(CppUnsigIntVector) vector; + %template(CppIntVector) vector; + %template(CppBoolVector) vector; %template(CppElementVector) vector; %template(CppDoublePosVector) vector< Pos >; %template(CppDoubleMatrix) vector< vector >; diff --git a/src/commands.cpp b/src/commands.cpp index c2d6944..2f553a7 100644 --- a/src/commands.cpp +++ b/src/commands.cpp @@ -1084,7 +1084,7 @@ int cmd_track_linepass(const std::vector& args) { std::vector> pos_list; Plane::type lost_plane; unsigned int offset_element = start_element; - track_linepass(accelerator, pos, pos_list, offset_element, lost_plane, true, 0, {0,}, {0.0, 0.0, }); + track_linepass(accelerator, pos, true, offset_element, pos_list, lost_plane, 0, {0,}, {0.0, 0.0, }); std::cout << get_timestamp() << " saving track_linepass data to file" << std::endl; status = print_tracking_linepass(accelerator, pos_list, start_element, "track_linepass_out.txt"); diff --git a/src/diffusion_matrix.cpp b/src/diffusion_matrix.cpp index c994bf4..a8ed280 100644 --- a/src/diffusion_matrix.cpp +++ b/src/diffusion_matrix.cpp @@ -256,7 +256,7 @@ Status::type track_diffusionmatrix (const Accelerator& accelerator, bdiff.sandwichme_with_matrix(tm[i]); bmat.push_back(bdiff); // track through element - track_elementpass (ele, fp, accelerator); + track_elementpass(accelerator, ele, fp); } return status; } diff --git a/src/optics.cpp b/src/optics.cpp index eb16a28..de4e31b 100644 --- a/src/optics.cpp +++ b/src/optics.cpp @@ -87,7 +87,7 @@ Status::type calc_twiss(Accelerator& accelerator, std::vector TAW_indices; double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); - Status::type status = track_linepass(accelerator, fp, closed_orbit, element_offset, lost_plane, true, accelerator_length, TAW_indices, TAW_positions); + Status::type status = track_linepass(accelerator, fp, true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions); if (status != Status::success) return status; #ifdef TIMEIT @@ -153,7 +153,8 @@ Status::type calc_twiss(Accelerator& accelerator, fpp.px += twiss0.etax[1] * dpp; fpp.ry += twiss0.etay[0] * dpp; fpp.py += twiss0.etay[1] * dpp; - Status::type status = track_linepass(accelerator, fpp, codp, element_offset, lost_plane, true, accelerator_length, TAW_indices, TAW_positions); + Status::type status = track_linepass(accelerator, fpp, true, element_offset, codp, lost_plane, accelerator_length, TAW_indices, TAW_positions + ); if (status != Status::success) return status; } diff --git a/src/tests.cpp b/src/tests.cpp index d30a451..10736c5 100644 --- a/src/tests.cpp +++ b/src/tests.cpp @@ -65,7 +65,7 @@ int test_linepass_tpsa(const Accelerator& accelerator, const std::vector > > new_tpsa; unsigned int element_offset = 0; Plane::type lost_plane; - track_linepass(accelerator, tpsa, new_tpsa, element_offset, lost_plane, false, 0, {0,}, {0.0, 0.0, }); + track_linepass(accelerator, tpsa, false, element_offset, new_tpsa, lost_plane, 0, {0,}, {0.0, 0.0, }); for(unsigned int i=0; i >& c = new_particles[i]; //std::cout << c.rx << std::endl; @@ -84,13 +84,23 @@ int test_ringpass(const Accelerator& accelerator) { pos.ry = 0.00010; std::vector > new_pos; - unsigned int element_offset = 0, lost_turn = 0; + unsigned int element_offset = 0 + unsigned int lost_turn = 0; Plane::type lost_plane; clock_t begin, end; double time_spent; begin = clock(); - Status::type status = track_ringpass(accelerator, pos, new_pos, 5000, lost_turn, element_offset, lost_plane, true); + Status::type status = track_ringpass( + accelerator, + pos, + 5000, + true, + element_offset, + new_pos, + lost_plane, + lost_turn + ); end = clock(); time_spent = (double)(end - begin) / CLOCKS_PER_SEC; @@ -339,7 +349,7 @@ int test_simple_drift() { accelerator.lattice.push_back(ds); Pos pos(0.001,0.002,0.003,0.004,0.005,0.006); - track_elementpass (ds, pos, accelerator); + track_elementpass (accelerator, ds, pos); fprintf(stdout, "test_simple_drift\n"); fprintf(stdout, "rx: %+.16f\n", pos.rx); @@ -367,7 +377,7 @@ int test_simple_quadrupole() { accelerator.lattice.push_back(ds); Pos pos(0.001,0.002,0.003,0.004,0.005,0.006); - track_elementpass (ds, pos, accelerator); + track_elementpass (accelerator, ds, pos); fprintf(stdout, "test_simple_quadrupole\n"); fprintf(stdout, "rx: %+.16f\n", pos.rx); @@ -398,12 +408,17 @@ int test_linepass2() { orig_pos.rx = 0.0001; orig_pos.px = 0.0001; - Status::type status = track_linepass (accelerator, - orig_pos, // initial electron coordinates - pos, // vector with electron coordinates from tracking at every element. - element_offset, // index of starting element for tracking - lost_plane, // return plane in which particle was lost, if the case. - trajectory, 0, {0,}, {0.0, 0.0, }); + Status::type status = track_linepass ( + accelerator, + orig_pos, + trajectory, + element_offset, + pos, + lost_plane, + 0, + {0,}, + {0.0, 0.0, } + ); for(unsigned int i=0; i<10; ++i) { std::cout << std::endl; diff --git a/src/tracking.cpp b/src/tracking.cpp index 446be00..a908145 100644 --- a/src/tracking.cpp +++ b/src/tracking.cpp @@ -94,7 +94,7 @@ Status::type track_findm66 (Accelerator& accelerator, map.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; TAW_pivot++; } - if ((status = track_elementpass (lattice[i], map, accelerator)) != Status::success) return status; + if ((status = track_elementpass (accelerator, lattice[i], map)) != Status::success) return status; if (i == time_aware_element_indices.back()) { map.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; @@ -185,13 +185,27 @@ Status::type track_findorbit6( Plane::type lost_plane; Status::type status = Status::success; - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[0], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[1], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[2], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[3], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[4], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[5], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[6], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[0], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[1], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[2], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[3], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[4], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[5], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[6], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); if (status != Status::success) { return Status::findorbit_one_turn_matrix_problem; @@ -226,7 +240,9 @@ Status::type track_findorbit6( closed_orbit.clear(); unsigned int element_offset = 0; Plane::type lost_plane; - track_linepass(accelerator, co[6], closed_orbit, element_offset, lost_plane, true, accelerator_length, TAW_indices, TAW_positions); + track_linepass( + accelerator, co[6], true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions + ); accelerator.radiation_on = radsts; return Status::success; @@ -272,11 +288,21 @@ Status::type track_findorbit4( unsigned int element_offset = 0; Plane::type lost_plane; Status::type status = Status::success; - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[0], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[1], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[2], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[3], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); - status = (Status::type) ((int) status | (int) track_linepass(accelerator, co[6], co2, element_offset, lost_plane, false, accelerator_length, TAW_indices, TAW_positions)); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[0], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[1], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[2], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[3], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); + status = (Status::type) ((int) status | (int) track_linepass( + accelerator, co[6], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions + )); if (status != Status::success) { return Status::findorbit_one_turn_matrix_problem; } @@ -305,7 +331,9 @@ Status::type track_findorbit4( closed_orbit.clear(); unsigned int element_offset = 0; Plane::type lost_plane; - track_linepass(accelerator, co[6], closed_orbit, element_offset, lost_plane, true, accelerator_length, TAW_indices, TAW_positions); + track_linepass( + accelerator, co[6], true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions + ); accelerator.radiation_on = radsts; return Status::success; diff --git a/tests/test-kickmap.cpp b/tests/test-kickmap.cpp index 0808742..7492d42 100644 --- a/tests/test-kickmap.cpp +++ b/tests/test-kickmap.cpp @@ -52,7 +52,7 @@ int main() { fixed_point.rx = 0.001; std::cout << fixed_point << std::endl; - track_elementpass(accelerator.lattice[2], fixed_point, accelerator); + track_elementpass(accelerator, accelerator.lattice[2], fixed_point); std::cout << fixed_point << std::endl; return 0; @@ -65,12 +65,14 @@ int main() { Plane::type lost_plane; unsigned int element_offset = 0; - status = track_linepass(accelerator, fp, closed_orbit, element_offset, lost_plane, true); + status = track_linepass( + accelerator, fp, true, element_offset, closed_orbit, lost_plane, 0, {0,}, {0.0, 0.0, } + ); if (status != Status::success) return status; std::vector atm; Pos v0; - status = track_findm66(accelerator, closed_orbit[0], atm, m66, v0); + status = track_findm66(accelerator, closed_orbit[0], atm, m66, v0, 0, {0,}, {0.0, 0.0, }); if (status != Status::success) return status; std::cout << m66[0][0] + m66[1][1] << std::endl; std::cout << m66[2][2] + m66[3][3] << std::endl; @@ -116,7 +118,9 @@ int main() { fpp.px += twiss0.etax[1] * dpp; fpp.ry += twiss0.etay[0] * dpp; fpp.py += twiss0.etay[1] * dpp; - Status::type status = track_linepass(accelerator, fpp, codp, element_offset, lost_plane, true); + Status::type status = track_linepass( + accelerator, fpp, true, element_offset, codp, lost_plane, 0, {0,}, {0.0, 0.0, } + ); std::cout << "h2" << std::endl; if (status != Status::success) return status; } @@ -132,7 +136,9 @@ int main() { // std::vector> closed_orbit, traj; // Plane::type lost_plane; // unsigned int element_offset = 0; - // Status::type status1 = track_linepass(sirius, fp, traj, element_offset, lost_plane, true); + // Status::type status1 = track_linepass( + // sirius, fp, true, element_offset, traj, lost_plane, 0, {0,}, {0.0, 0.0, } + // ); // for(auto i=0; i <= closed_orbit) // std::cout << status1 << std::endl; From 3011a916fb26f5f424b168cec2133c046eed0630 Mon Sep 17 00:00:00 2001 From: vitor Date: Sat, 31 Aug 2024 19:07:30 -0300 Subject: [PATCH 13/20] FIX: duplicated lines --- python_package/interface.cpp | 44 +----------------------------------- python_package/interface.h | 16 ------------- src/tests.cpp | 13 ++--------- 3 files changed, 3 insertions(+), 70 deletions(-) diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 7dadcbe..8e8a3ff 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -23,10 +23,6 @@ Status::type track_elementpass_wrapper ( const Element& el, double *pos, int n1, int n2 ) { - const Accelerator& accelerator, - const Element& el, - double *pos, int n1, int n2 -) { std::vector > post; @@ -37,7 +33,7 @@ Status::type track_elementpass_wrapper ( pos[2*n2 + i], pos[3*n2 + i], pos[4*n2 + i], pos[5*n2 + i]); } - Status::type status = track_elementpass(accelerator, el, post); + Status::type status = track_elementpass(accelerator, el, post); for (unsigned int i=0; i& fixed_point, - Matrix& m66, - double *twiss, int n1, int n2, - Twiss twiss0 -); Element marker_wrapper(const std::string& fam_name_); Element bpm_wrapper(const std::string& fam_name_); diff --git a/src/tests.cpp b/src/tests.cpp index b5e9053..a4ff9f3 100644 --- a/src/tests.cpp +++ b/src/tests.cpp @@ -36,7 +36,8 @@ int test_linepass(const Accelerator& accelerator) { std::vector > new_pos; unsigned int element_offset = 0; Plane::type lost_plane; - Status::type status = track_linepass(accelerator, pos, new_pos, element_offset, lost_plane, true, 0, {0,}, {0.0, 0.0, }); + // Status::type status = track_linepass(accelerator, pos, new_pos, element_offset, lost_plane, true, 0, {0,}, {0.0, 0.0, }); + Status::type status = track_linepass(accelerator, pos, true, element_offset, new_pos, lost_plane, 0, {0,}, {0.0, 0.0, }); std::cout << "status: " << string_error_messages[status] << std::endl; @@ -101,16 +102,6 @@ int test_ringpass(const Accelerator& accelerator) { lost_plane, lost_turn ); - Status::type status = track_ringpass( - accelerator, - pos, - 5000, - true, - element_offset, - new_pos, - lost_plane, - lost_turn - ); end = clock(); time_spent = (double)(end - begin) / CLOCKS_PER_SEC; From b41341a824c4acdb13c8de1f22ff84896e05632f Mon Sep 17 00:00:00 2001 From: vitor Date: Mon, 2 Sep 2024 10:26:55 -0300 Subject: [PATCH 14/20] Cleaning (diff) --- include/trackcpp/passmethods.hpp | 5 +- include/trackcpp/tracking.h | 100 ++++++++++++++----------------- python_package/interface.cpp | 4 +- src/commands.cpp | 1 - src/kicktable.cpp | 1 + src/tests.cpp | 5 +- src/tracking.cpp | 4 +- 7 files changed, 54 insertions(+), 66 deletions(-) diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 64f7a97..cd2a9df 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -400,10 +400,9 @@ Status::type pm_cavity_pass(Pos &pos, const Element &elem, global_2_local(pos, elem); double nv = elem.voltage / accelerator.energy; - double frf = elem.frequency; if (elem.length == 0) { T &de = pos.de, &dl = pos.dl; - de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); + de += -nv * sin(TWOPI*elem.frequency *dl/ light_speed - elem.phase_lag); } else { T &rx = pos.rx, &px = pos.px; T &ry = pos.ry, &py = pos.py; @@ -415,7 +414,7 @@ Status::type pm_cavity_pass(Pos &pos, const Element &elem, ry += norml * py; dl += 0.5 * norml * pnorm * (px*px + py*py); // longitudinal momentum kick - de += -nv * sin(TWOPI*frf *dl/ light_speed - elem.phase_lag); + de += -nv * sin(TWOPI*elem.frequency*dl/light_speed - elem.phase_lag); // drift half length pnorm = 1.0 / (1.0 + de); norml = (0.5 * elem.length) * pnorm; diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index b26f665..a87b651 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -59,9 +59,9 @@ Pos linalg_solve6_posvec( template Status::type track_elementpass ( - const Accelerator& accelerator, - const Element& el, // element through which to track particle - Pos &orig_pos // initial electron coordinates + const Accelerator& accelerator, + const Element& el, // element through which to track particle + Pos &orig_pos // initial electron coordinates ) { Status::type status = Status::success; @@ -108,9 +108,9 @@ Status::type track_elementpass ( template Status::type track_elementpass ( - const Accelerator& accelerator, - const Element& el, // element through which to track particle - std::vector >& orig_pos // initial electron coordinates + const Accelerator& accelerator, + const Element& el, // element through which to track particle + std::vector >& orig_pos // initial electron coordinates ) { Status::type status = Status::success; @@ -145,15 +145,15 @@ Status::type track_elementpass ( // template Status::type track_linepass ( - const Accelerator& accelerator, - Pos& orig_pos, - const std::vector& indices, - unsigned int& element_offset, - std::vector >& pos, - Plane::type& lost_plane, - const double line_length, - const std::vector& time_aware_element_indices, - const std::vector& time_aware_element_positions + const Accelerator& accelerator, + Pos& orig_pos, + const std::vector& indices, + unsigned int& element_offset, + std::vector >& pos, + Plane::type& lost_plane, + const double line_length, + const std::vector& time_aware_element_indices, + const std::vector& time_aware_element_positions ) { Status::type status = Status::success; @@ -170,12 +170,6 @@ Status::type track_linepass ( for (unsigned int i=0; i<=nr_elements; ++i) indcs[i] = false; for (auto&& i: indices) if (i<=nr_elements) indcs[i] = true; - // std::cout << "len acc = " << line_length << std::endl; - // for (size_t i; i& args) { } std::cout << get_timestamp() << " input file with flat lattice read." << std::endl; - // builds accelerator accelerator.energy = ring_energy; accelerator.harmonic_number = harmonic_number; diff --git a/src/kicktable.cpp b/src/kicktable.cpp index 4495bfc..0034cba 100644 --- a/src/kicktable.cpp +++ b/src/kicktable.cpp @@ -104,6 +104,7 @@ Status::type Kicktable::load_from_file(const std::string& filename_) { } int add_kicktable(const std::string& filename) { + // looks through vector of kicktables... for(unsigned int i=0; i > new_pos; unsigned int element_offset = 0; Plane::type lost_plane; - // Status::type status = track_linepass(accelerator, pos, new_pos, element_offset, lost_plane, true, 0, {0,}, {0.0, 0.0, }); Status::type status = track_linepass(accelerator, pos, true, element_offset, new_pos, lost_plane, 0, {0,}, {0.0, 0.0, }); std::cout << "status: " << string_error_messages[status] << std::endl; @@ -350,7 +349,7 @@ int test_simple_drift() { accelerator.lattice.push_back(ds); Pos pos(0.001,0.002,0.003,0.004,0.005,0.006); - track_elementpass (accelerator, ds, pos); + track_elementpass(accelerator, ds, pos); fprintf(stdout, "test_simple_drift\n"); fprintf(stdout, "rx: %+.16f\n", pos.rx); @@ -378,7 +377,7 @@ int test_simple_quadrupole() { accelerator.lattice.push_back(ds); Pos pos(0.001,0.002,0.003,0.004,0.005,0.006); - track_elementpass (accelerator, ds, pos); + track_elementpass(accelerator, ds, pos); fprintf(stdout, "test_simple_quadrupole\n"); fprintf(stdout, "rx: %+.16f\n", pos.rx); diff --git a/src/tracking.cpp b/src/tracking.cpp index 150aa9f..4ccfb26 100644 --- a/src/tracking.cpp +++ b/src/tracking.cpp @@ -50,7 +50,6 @@ Status::type track_findm66 (Accelerator& accelerator, const std::vector& lattice = accelerator.lattice; Pos fp = fixed_point; - const int radsts = accelerator.radiation_on; if (radsts == RadiationState::full){ accelerator.radiation_on = RadiationState::damping; @@ -121,7 +120,9 @@ Status::type track_findm66 (Accelerator& accelerator, // constant term of the final map v0.rx = map.rx.c[0]; v0.px = map.px.c[0]; v0.ry = map.ry.c[0]; v0.py = map.py.c[0]; v0.de = map.de.c[0]; v0.dl = map.dl.c[0]; + accelerator.radiation_on = radsts; + return status; } @@ -167,7 +168,6 @@ Status::type track_findorbit6( std::vector TAW_indices; double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); - // temporary vectors and matrices std::vector > co(7,0); for(auto i=0; i<7; ++i) co[i] = fixed_point_guess; From 3b5052c248419d7be5584008bd87d130ccd44a01 Mon Sep 17 00:00:00 2001 From: vitor Date: Mon, 2 Sep 2024 10:38:06 -0300 Subject: [PATCH 15/20] rm diff --- include/trackcpp/passmethods.hpp | 3 ++- include/trackcpp/tracking.h | 3 +-- src/optics.cpp | 7 +++++-- src/tests.cpp | 8 ++++++-- src/tracking.cpp | 5 ++--- tests/test-kickmap.cpp | 1 - 6 files changed, 16 insertions(+), 11 deletions(-) diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index cd2a9df..9b315ea 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -402,7 +402,7 @@ Status::type pm_cavity_pass(Pos &pos, const Element &elem, double nv = elem.voltage / accelerator.energy; if (elem.length == 0) { T &de = pos.de, &dl = pos.dl; - de += -nv * sin(TWOPI*elem.frequency *dl/ light_speed - elem.phase_lag); + de += -nv * sin(TWOPI*elem.frequency * dl/ light_speed - elem.phase_lag); } else { T &rx = pos.rx, &px = pos.px; T &ry = pos.ry, &py = pos.py; @@ -440,6 +440,7 @@ Status::type pm_thinsext_pass(Pos &pos, const Element &elem, return Status::passmethod_not_implemented; } + template Status::type pm_kickmap_pass(Pos &pos, const Element &elem, const Accelerator& accelerator) { diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index a87b651..f9b6e93 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -111,8 +111,7 @@ Status::type track_elementpass ( const Accelerator& accelerator, const Element& el, // element through which to track particle std::vector >& orig_pos // initial electron coordinates -) { - +){ Status::type status = Status::success; for(auto&& pos: orig_pos) { Status::type status2 = track_elementpass(accelerator, el, pos); diff --git a/src/optics.cpp b/src/optics.cpp index de4e31b..2316eb1 100644 --- a/src/optics.cpp +++ b/src/optics.cpp @@ -87,7 +87,9 @@ Status::type calc_twiss(Accelerator& accelerator, std::vector TAW_indices; double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); - Status::type status = track_linepass(accelerator, fp, true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions); + Status::type status = track_linepass( + accelerator, fp, true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions + ); if (status != Status::success) return status; #ifdef TIMEIT @@ -153,7 +155,8 @@ Status::type calc_twiss(Accelerator& accelerator, fpp.px += twiss0.etax[1] * dpp; fpp.ry += twiss0.etay[0] * dpp; fpp.py += twiss0.etay[1] * dpp; - Status::type status = track_linepass(accelerator, fpp, true, element_offset, codp, lost_plane, accelerator_length, TAW_indices, TAW_positions + Status::type status = track_linepass( + accelerator, fpp, true, element_offset, codp, lost_plane, accelerator_length, TAW_indices, TAW_positions ); if (status != Status::success) return status; } diff --git a/src/tests.cpp b/src/tests.cpp index cd1c51b..009e233 100644 --- a/src/tests.cpp +++ b/src/tests.cpp @@ -36,7 +36,9 @@ int test_linepass(const Accelerator& accelerator) { std::vector > new_pos; unsigned int element_offset = 0; Plane::type lost_plane; - Status::type status = track_linepass(accelerator, pos, true, element_offset, new_pos, lost_plane, 0, {0,}, {0.0, 0.0, }); + Status::type status = track_linepass( + accelerator, pos, true, element_offset, new_pos, lost_plane, 0, {0,}, {0.0, 0.0, } + ); std::cout << "status: " << string_error_messages[status] << std::endl; @@ -65,7 +67,9 @@ int test_linepass_tpsa(const Accelerator& accelerator, const std::vector > > new_tpsa; unsigned int element_offset = 0; Plane::type lost_plane; - track_linepass(accelerator, tpsa, false, element_offset, new_tpsa, lost_plane, 0, {0,}, {0.0, 0.0, }); + track_linepass( + accelerator, tpsa, false, element_offset, new_tpsa, lost_plane, 0, {0,}, {0.0, 0.0, } + ); for(unsigned int i=0; i >& c = new_particles[i]; //std::cout << c.rx << std::endl; diff --git a/src/tracking.cpp b/src/tracking.cpp index 4ccfb26..32bd76e 100644 --- a/src/tracking.cpp +++ b/src/tracking.cpp @@ -188,7 +188,7 @@ Status::type track_findorbit6( status = (Status::type) ((int) status | (int) track_linepass( accelerator, co[0], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions - )); + )); status = (Status::type) ((int) status | (int) track_linepass( accelerator, co[1], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions )); @@ -249,7 +249,6 @@ Status::type track_findorbit6( } - Status::type track_findorbit4( Accelerator& accelerator, std::vector >& closed_orbit, @@ -291,7 +290,7 @@ Status::type track_findorbit4( Status::type status = Status::success; status = (Status::type) ((int) status | (int) track_linepass( accelerator, co[0], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions - )); + )); status = (Status::type) ((int) status | (int) track_linepass( accelerator, co[1], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions )); diff --git a/tests/test-kickmap.cpp b/tests/test-kickmap.cpp index e6f4ade..7492d42 100644 --- a/tests/test-kickmap.cpp +++ b/tests/test-kickmap.cpp @@ -53,7 +53,6 @@ int main() { fixed_point.rx = 0.001; std::cout << fixed_point << std::endl; track_elementpass(accelerator, accelerator.lattice[2], fixed_point); - track_elementpass(accelerator, accelerator.lattice[2], fixed_point); std::cout << fixed_point << std::endl; return 0; From ee035b75eb917e03073e8fe8449347f2440b09cc Mon Sep 17 00:00:00 2001 From: Vitor Date: Thu, 3 Apr 2025 11:05:49 -0300 Subject: [PATCH 16/20] error fix , lost previous local commit --- include/trackcpp/tracking.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index f9b6e93..feae000 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -112,7 +112,7 @@ Status::type track_elementpass ( const Element& el, // element through which to track particle std::vector >& orig_pos // initial electron coordinates ){ - Status::type status = Status::success; + Status::type status = Status::success; for(auto&& pos: orig_pos) { Status::type status2 = track_elementpass(accelerator, el, pos); if (status2 != Status::success) status = status2; @@ -178,7 +178,7 @@ Status::type track_linepass ( // stores trajectory at entrance of each element if (indcs[i]) pos.push_back(orig_pos); - if (i == time_aware_element_indices[TAW_pivot]) { + if (element_offset == time_aware_element_indices[TAW_pivot]) { ddl = light_speed*accelerator.harmonic_number/element.frequency - line_length; orig_pos.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; TAW_pivot++; @@ -186,7 +186,7 @@ Status::type track_linepass ( status = track_elementpass(accelerator, element, orig_pos); - if (i == time_aware_element_indices.back()) { + if (element_offset == time_aware_element_indices.back()) { orig_pos.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; } From d8d84d9313b2b1636f61eb0fcc693c4bbb1916f6 Mon Sep 17 00:00:00 2001 From: Vitor Date: Thu, 3 Apr 2025 14:58:15 -0300 Subject: [PATCH 17/20] avoid empty vector to be acessed --- include/trackcpp/tracking.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index feae000..63b9712 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -186,7 +186,7 @@ Status::type track_linepass ( status = track_elementpass(accelerator, element, orig_pos); - if (element_offset == time_aware_element_indices.back()) { + if (time_aware_element_indices.size() > 0 && element_offset == time_aware_element_indices.back()) { orig_pos.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; } From ea7902d02b94892b696dfd6e0e28be4fa2c36d6d Mon Sep 17 00:00:00 2001 From: Vitor Date: Wed, 16 Apr 2025 18:08:05 -0300 Subject: [PATCH 18/20] element_offset fix --- include/trackcpp/accelerator.h | 2 +- include/trackcpp/tracking.h | 31 +++++---- python_package/interface.cpp | 16 ++--- python_package/interface.h | 8 +-- src/accelerator.cpp | 121 ++++++++++++++++++++++++++++++--- src/commands.cpp | 8 ++- src/optics.cpp | 2 +- src/tests.cpp | 25 +++++-- src/tracking.cpp | 47 ++++++------- tests/test-kickmap.cpp | 13 ++-- 10 files changed, 193 insertions(+), 80 deletions(-) diff --git a/include/trackcpp/accelerator.h b/include/trackcpp/accelerator.h index 5ce1a2a..9ad28d3 100644 --- a/include/trackcpp/accelerator.h +++ b/include/trackcpp/accelerator.h @@ -40,7 +40,7 @@ class Accelerator { bool isequal(const Accelerator& a) const { return *this == a; } // necessary for python_package double get_length() const; friend std::ostream& operator<< (std::ostream &out, const Accelerator& a); - double get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions) const; + double get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions, unsigned int element_offset = 0) const; }; #endif diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index 63b9712..b57cf04 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -29,18 +29,12 @@ Status::type track_findm66( Accelerator& accelerator, const Pos& fixed_point, - std::vector& tm, Matrix& m66, Pos& v0, - const double line_length, - const std::vector& time_aware_element_indices, - const std::vector& time_aware_element_positions); + std::vector& tm, Matrix& m66, Pos& v0); Status::type track_findm66( Accelerator& accelerator, const Pos& fixed_point, std::vector& tm, Matrix& m66, Pos& v0, - std::vector& indices, - const double line_length, - const std::vector& time_aware_element_indices, - const std::vector& time_aware_element_positions); + std::vector& indices); Status::type track_findorbit4( Accelerator& accelerator, std::vector >& closed_orbit, @@ -180,15 +174,24 @@ Status::type track_linepass ( if (element_offset == time_aware_element_indices[TAW_pivot]) { ddl = light_speed*accelerator.harmonic_number/element.frequency - line_length; - orig_pos.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; - TAW_pivot++; + + // std::cout << "elem idx = " << element_offset << " taw_idx[pivot] = " << time_aware_element_indices[TAW_pivot] << std::endl; + // std::cout << "c = " << light_speed << " h = " << accelerator.harmonic_number << " elem freq = " << element.frequency << " line length = " << line_length << std::endl; + // std::cout << "ddl = " << ddl << std::endl; + // std::cout << "dl before = " << orig_pos.dl << std::endl; + + orig_pos.dl -= ddl * time_aware_element_positions[TAW_pivot] / line_length; + + // std::cout << "dl after = " << orig_pos.dl << std::endl; + + if (TAW_pivot < time_aware_element_indices.size() - 1) {TAW_pivot++;} } status = track_elementpass(accelerator, element, orig_pos); - if (time_aware_element_indices.size() > 0 && element_offset == time_aware_element_indices.back()) { - orig_pos.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; - } + // if (time_aware_element_indices.size() > 0 && i == time_aware_element_indices.back()) { + // orig_pos.dl -= ddl * (time_aware_element_positions[TAW_pivot+1]-time_aware_element_positions[TAW_pivot]) / line_length; + // } lost_plane = check_particle_loss(accelerator, element, orig_pos); if (lost_plane != Plane::no_plane) status = Status::particle_lost; @@ -402,7 +405,7 @@ Status::type track_ringpass ( // for longitudinal kick before RF cavities std::vector TAW_positions; std::vector TAW_indices; - double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, element_offset); if (turn_by_turn) pos.reserve(nr_turns+1); diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 8657666..66ab689 100644 --- a/python_package/interface.cpp +++ b/python_package/interface.cpp @@ -52,6 +52,9 @@ Status::type track_linepass_wrapper( std::vector > post; std::vector > orig_post; + std::vector TAW_positions; + std::vector TAW_indices; + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, args.element_offset); orig_post.reserve(ni2); for (unsigned int i=0; i& v0, - std::vector& indices, - double line_length, - std::vector TAW_indices, - std::vector TAW_positions) { + std::vector& indices) { std::vector vec_tm; Matrix vec_m66; Status::type status = track_findm66( - accelerator, fixed_point, vec_tm, vec_m66, v0, indices, line_length, TAW_indices, TAW_positions); + accelerator, fixed_point, vec_tm, vec_m66, v0, indices); for (unsigned int i=0; i lost_plane; std::vector lost_element; std::vector lost_flag; - double line_length; - std::vector time_aware_element_indices; - std::vector time_aware_element_positions; }; struct RingPassArgs : public LinePassArgs { @@ -108,10 +105,7 @@ Status::type track_findm66_wrapper( 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, - double line_length, - std::vector TAW_indices, - std::vector TAW_positions); + Pos& v0, std::vector& indices); Status::type track_diffusionmatrix_wrapper( diff --git a/src/accelerator.cpp b/src/accelerator.cpp index 6abb9e1..2bbbafd 100644 --- a/src/accelerator.cpp +++ b/src/accelerator.cpp @@ -53,26 +53,129 @@ std::ostream& operator<< (std::ostream &out, const Accelerator& a) { return out; } -double Accelerator::get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions) const { +// double Accelerator::get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions, unsigned int element_offset) const { +// // for longitudinal kick before RF cavities +// TAW_positions.clear(); +// TAW_indices.clear(); +// size_t nr_elements = this->lattice.size(); +// PassMethodsClass PMClass = PassMethodsClass(); +// int elem_offset = element_offset; +// for (size_t i = 0; i < nr_elements; i++) { +// auto elem = this->lattice[elem_offset]; +// if (PMClass.is_time_aware_pm(elem.pass_method)) {TAW_indices.push_back(i);} +// elem_offset = (elem_offset + 1) % nr_elements; +// } +// double s_pos = 0.0; +// TAW_positions.push_back(s_pos); +// for (size_t i = 0; i < nr_elements; i++) { +// auto elem = this->lattice[i]; +// if (PMClass.is_time_aware_pm(elem.pass_method)) { +// s_pos += (elem.length*0.5); +// TAW_positions.push_back(s_pos); +// s_pos += (elem.length*0.5); +// } +// else { +// s_pos += elem.length; +// } +// } +// TAW_positions.push_back(s_pos); +// return s_pos; +// } + +// double Accelerator::get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions, unsigned int element_offset) const { +// // for longitudinal kick before RF cavities +// TAW_positions.clear(); +// TAW_indices.clear(); +// size_t nr_elements = this->lattice.size(); +// PassMethodsClass PMClass = PassMethodsClass(); + +// double s_pos = 0.0; +// for (size_t i = element_offset; i < nr_elements; i++) { +// auto elem = this->lattice[i]; +// s_pos += elem.length; +// if (PMClass.is_time_aware_pm(elem.pass_method)) { +// TAW_indices.push_back(i); +// TAW_positions.push_back(s_pos); +// s_pos = 0.0; +// } +// } +// TAW_positions.push_back(s_pos); + +// s_pos = 0.0; +// std::vector temp_pos = {}; +// for (size_t i = 0; i < element_offset; i++) { +// auto elem = this->lattice[i]; +// s_pos += elem.length; +// if (PMClass.is_time_aware_pm(elem.pass_method)) { +// TAW_indices.push_back(i); +// temp_pos.push_back(s_pos); +// s_pos = 0.0; +// } +// } +// temp_pos.push_back(s_pos); + +// TAW_positions[0] += temp_pos.back(); +// temp_pos.pop_back(); + +// if (temp_pos.size() >= 1) { +// TAW_positions.back() += temp_pos[0]; +// temp_pos.erase(temp_pos.begin()); +// } +// if (TAW_positions.size() > TAW_indices.size()) { +// TAW_positions[0] += TAW_positions.back(); +// TAW_positions.pop_back(); +// } +// for (double s: temp_pos) { +// TAW_positions.push_back(s); +// } +// s_pos = 0.0; +// for (double s: TAW_positions) { +// s_pos += s; +// } +// temp_pos.clear(); + +// return s_pos; + +// } + + +double Accelerator::get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions, unsigned int element_offset) const { // for longitudinal kick before RF cavities TAW_positions.clear(); TAW_indices.clear(); size_t nr_elements = this->lattice.size(); PassMethodsClass PMClass = PassMethodsClass(); + + std::vector temp_positions; + std::vector temp_indices; double s_pos = 0.0; - TAW_positions.push_back(s_pos); + double acclen = 0.0; for (size_t i = 0; i < nr_elements; i++) { auto elem = this->lattice[i]; + acclen += elem.length; if (PMClass.is_time_aware_pm(elem.pass_method)) { - s_pos += (elem.length*0.5); - TAW_positions.push_back(s_pos); - TAW_indices.push_back(i); - s_pos += (elem.length*0.5); + temp_indices.push_back(i); + temp_positions.push_back(s_pos + elem.length/2); + s_pos = 0.0 + elem.length/2; } else { - s_pos += elem.length; + s_pos += elem.length; } + + } + temp_positions[0] += s_pos; + + size_t split = 0; + for (; split < temp_indices.size(); split++) { + if (temp_indices[split] >= element_offset) {break;} } - TAW_positions.push_back(s_pos); - return s_pos; + + TAW_indices.insert(TAW_indices.end(), temp_indices.begin() + split, temp_indices.end()); + TAW_indices.insert(TAW_indices.end(), temp_indices.begin(), temp_indices.begin() + split); + + TAW_positions.insert(TAW_positions.end(), temp_positions.begin() + split, temp_positions.end()); + TAW_positions.insert(TAW_positions.end(), temp_positions.begin(), temp_positions.begin() + split); + + return acclen; + } diff --git a/src/commands.cpp b/src/commands.cpp index 2f553a7..3f8e2b8 100644 --- a/src/commands.cpp +++ b/src/commands.cpp @@ -1084,7 +1084,13 @@ int cmd_track_linepass(const std::vector& args) { std::vector> pos_list; Plane::type lost_plane; unsigned int offset_element = start_element; - track_linepass(accelerator, pos, true, offset_element, pos_list, lost_plane, 0, {0,}, {0.0, 0.0, }); + + // for longitudinal kick before RF cavities + std::vector TAW_positions; + std::vector TAW_indices; + double acc_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, start_element); + + track_linepass(accelerator, pos, true, offset_element, pos_list, lost_plane, acc_length, TAW_indices, TAW_positions); std::cout << get_timestamp() << " saving track_linepass data to file" << std::endl; status = print_tracking_linepass(accelerator, pos_list, start_element, "track_linepass_out.txt"); diff --git a/src/optics.cpp b/src/optics.cpp index 2316eb1..76cb15a 100644 --- a/src/optics.cpp +++ b/src/optics.cpp @@ -103,7 +103,7 @@ Status::type calc_twiss(Accelerator& accelerator, std::vector atm; Pos v0; - status = track_findm66(accelerator, closed_orbit[0], atm, m66, v0, accelerator_length, TAW_indices, TAW_positions); + status = track_findm66(accelerator, closed_orbit[0], atm, m66, v0); if (status != Status::success) return status; #ifdef TIMEIT diff --git a/src/tests.cpp b/src/tests.cpp index 009e233..7cdc568 100644 --- a/src/tests.cpp +++ b/src/tests.cpp @@ -36,9 +36,12 @@ int test_linepass(const Accelerator& accelerator) { std::vector > new_pos; unsigned int element_offset = 0; Plane::type lost_plane; + // for longitudinal kick before RF cavities + std::vector TAW_positions; + std::vector TAW_indices; + double acc_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 37 Status::type status = track_linepass( - accelerator, pos, true, element_offset, new_pos, lost_plane, 0, {0,}, {0.0, 0.0, } - ); + accelerator, pos, true, element_offset, new_pos, lost_plane, acc_length, TAW_indices, TAW_positions); std::cout << "status: " << string_error_messages[status] << std::endl; @@ -67,9 +70,12 @@ int test_linepass_tpsa(const Accelerator& accelerator, const std::vector > > new_tpsa; unsigned int element_offset = 0; Plane::type lost_plane; + // for longitudinal kick before RF cavities + std::vector TAW_positions; + std::vector TAW_indices; + double acc_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 71 track_linepass( - accelerator, tpsa, false, element_offset, new_tpsa, lost_plane, 0, {0,}, {0.0, 0.0, } - ); + accelerator, tpsa, false, element_offset, new_tpsa, lost_plane, acc_length, TAW_indices, TAW_positions); for(unsigned int i=0; i >& c = new_particles[i]; //std::cout << c.rx << std::endl; @@ -409,6 +415,11 @@ int test_linepass2() { Plane::type lost_plane; bool trajectory = true; + // for longitudinal kick before RF cavities + std::vector TAW_positions; + std::vector TAW_indices; + double acc_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 414 + orig_pos.rx = 0.0001; orig_pos.px = 0.0001; @@ -419,9 +430,9 @@ int test_linepass2() { element_offset, pos, lost_plane, - 0, - {0,}, - {0.0, 0.0, } + acc_length, + TAW_indices, + TAW_positions ); for(unsigned int i=0; i<10; ++i) { diff --git a/src/tracking.cpp b/src/tracking.cpp index 32bd76e..9c4d732 100644 --- a/src/tracking.cpp +++ b/src/tracking.cpp @@ -41,10 +41,7 @@ Status::type track_findm66 (Accelerator& accelerator, std::vector& tm, Matrix& m66, Pos& v0, - std::vector& indices, - const double line_length, - const std::vector& time_aware_element_indices, - const std::vector& time_aware_element_positions) { + std::vector& indices) { Status::type status = Status::success; const std::vector& lattice = accelerator.lattice; @@ -70,6 +67,9 @@ Status::type track_findm66 (Accelerator& accelerator, tm.clear(); tm.reserve(indices.size()); unsigned int TAW_pivot = 0; + std::vector TAW_positions; + std::vector TAW_indices; + double line_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); double ddl = 0; for(unsigned int i=0; i& fixed_point, std::vector& tm, Matrix& m66, - Pos& v0, - const double line_length, - const std::vector& time_aware_element_indices, - const std::vector& time_aware_element_positions) { + Pos& v0) { std::vector indices; unsigned int nr_elements = accelerator.lattice.size(); @@ -143,7 +136,7 @@ Status::type track_findm66 (Accelerator& accelerator, indices.reserve(nr_elements + 1); for (unsigned int i=0; i<=nr_elements; ++i) indices.push_back(i); - return track_findm66 (accelerator, fixed_point, tm, m66, v0, indices, line_length, time_aware_element_indices, time_aware_element_positions); + return track_findm66 (accelerator, fixed_point, tm, m66, v0, indices); } @@ -163,11 +156,6 @@ Status::type track_findorbit6( accelerator.radiation_on = RadiationState::damping; } - // for longitudinal kick before RF cavities - std::vector TAW_positions; - std::vector TAW_indices; - double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); - // temporary vectors and matrices std::vector > co(7,0); for(auto i=0; i<7; ++i) co[i] = fixed_point_guess; @@ -177,6 +165,11 @@ Status::type track_findorbit6( Pos dco(1.0,1.0,1.0,1.0,1.0,1.0); matrix6_set_identity_posvec(D, delta); + // for longitudinal kick before RF cavities + std::vector TAW_positions; + std::vector TAW_indices; + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 182 + int nr_iter = 0; while ((get_max(dco) > tolerance) and (nr_iter <= max_nr_iters)) { co = co + D; @@ -265,11 +258,6 @@ Status::type track_findorbit4( accelerator.radiation_on = RadiationState::damping; } - // for longitudinal kick before RF cavities - std::vector TAW_positions; - std::vector TAW_indices; - double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions); - // temporary vectors and matrices // std::vector > co(7,0); // no initial guess std::vector > co(7,fixed_point_guess); @@ -280,6 +268,11 @@ Status::type track_findorbit4( Pos theta(0.0,0.0,0.0,0.0,0.0,0.0); matrix6_set_identity_posvec(D, delta); + // for longitudinal kick before RF cavities + std::vector TAW_positions; + std::vector TAW_indices; + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 285 + int nr_iter = 0; while ((get_max(dco) > tolerance) and (nr_iter <= max_nr_iters)) { co = co + D; diff --git a/tests/test-kickmap.cpp b/tests/test-kickmap.cpp index 7492d42..78681f0 100644 --- a/tests/test-kickmap.cpp +++ b/tests/test-kickmap.cpp @@ -65,14 +65,18 @@ int main() { Plane::type lost_plane; unsigned int element_offset = 0; + // for longitudinal kick before RF cavities + std::vector TAW_positions; + std::vector TAW_indices; + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 66 + status = track_linepass( - accelerator, fp, true, element_offset, closed_orbit, lost_plane, 0, {0,}, {0.0, 0.0, } - ); + accelerator, fp, true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions); if (status != Status::success) return status; std::vector atm; Pos v0; - status = track_findm66(accelerator, closed_orbit[0], atm, m66, v0, 0, {0,}, {0.0, 0.0, }); + status = track_findm66(accelerator, closed_orbit[0], atm, m66, v0); if (status != Status::success) return status; std::cout << m66[0][0] + m66[1][1] << std::endl; std::cout << m66[2][2] + m66[3][3] << std::endl; @@ -119,8 +123,7 @@ int main() { fpp.ry += twiss0.etay[0] * dpp; fpp.py += twiss0.etay[1] * dpp; Status::type status = track_linepass( - accelerator, fpp, true, element_offset, codp, lost_plane, 0, {0,}, {0.0, 0.0, } - ); + accelerator, fpp, true, element_offset, codp, lost_plane, accelerator_length, TAW_indices, TAW_positions); std::cout << "h2" << std::endl; if (status != Status::success) return status; } From a6120bc22b6e47cabe3d48a9bd9b3a1fc76a0d60 Mon Sep 17 00:00:00 2001 From: Vitor Date: Thu, 17 Apr 2025 11:12:48 -0300 Subject: [PATCH 19/20] del comment --- src/accelerator.cpp | 85 --------------------------------------------- 1 file changed, 85 deletions(-) diff --git a/src/accelerator.cpp b/src/accelerator.cpp index 2bbbafd..19c5eca 100644 --- a/src/accelerator.cpp +++ b/src/accelerator.cpp @@ -53,91 +53,6 @@ std::ostream& operator<< (std::ostream &out, const Accelerator& a) { return out; } -// double Accelerator::get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions, unsigned int element_offset) const { -// // for longitudinal kick before RF cavities -// TAW_positions.clear(); -// TAW_indices.clear(); -// size_t nr_elements = this->lattice.size(); -// PassMethodsClass PMClass = PassMethodsClass(); -// int elem_offset = element_offset; -// for (size_t i = 0; i < nr_elements; i++) { -// auto elem = this->lattice[elem_offset]; -// if (PMClass.is_time_aware_pm(elem.pass_method)) {TAW_indices.push_back(i);} -// elem_offset = (elem_offset + 1) % nr_elements; -// } -// double s_pos = 0.0; -// TAW_positions.push_back(s_pos); -// for (size_t i = 0; i < nr_elements; i++) { -// auto elem = this->lattice[i]; -// if (PMClass.is_time_aware_pm(elem.pass_method)) { -// s_pos += (elem.length*0.5); -// TAW_positions.push_back(s_pos); -// s_pos += (elem.length*0.5); -// } -// else { -// s_pos += elem.length; -// } -// } -// TAW_positions.push_back(s_pos); -// return s_pos; -// } - -// double Accelerator::get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions, unsigned int element_offset) const { -// // for longitudinal kick before RF cavities -// TAW_positions.clear(); -// TAW_indices.clear(); -// size_t nr_elements = this->lattice.size(); -// PassMethodsClass PMClass = PassMethodsClass(); - -// double s_pos = 0.0; -// for (size_t i = element_offset; i < nr_elements; i++) { -// auto elem = this->lattice[i]; -// s_pos += elem.length; -// if (PMClass.is_time_aware_pm(elem.pass_method)) { -// TAW_indices.push_back(i); -// TAW_positions.push_back(s_pos); -// s_pos = 0.0; -// } -// } -// TAW_positions.push_back(s_pos); - -// s_pos = 0.0; -// std::vector temp_pos = {}; -// for (size_t i = 0; i < element_offset; i++) { -// auto elem = this->lattice[i]; -// s_pos += elem.length; -// if (PMClass.is_time_aware_pm(elem.pass_method)) { -// TAW_indices.push_back(i); -// temp_pos.push_back(s_pos); -// s_pos = 0.0; -// } -// } -// temp_pos.push_back(s_pos); - -// TAW_positions[0] += temp_pos.back(); -// temp_pos.pop_back(); - -// if (temp_pos.size() >= 1) { -// TAW_positions.back() += temp_pos[0]; -// temp_pos.erase(temp_pos.begin()); -// } -// if (TAW_positions.size() > TAW_indices.size()) { -// TAW_positions[0] += TAW_positions.back(); -// TAW_positions.pop_back(); -// } -// for (double s: temp_pos) { -// TAW_positions.push_back(s); -// } -// s_pos = 0.0; -// for (double s: TAW_positions) { -// s_pos += s; -// } -// temp_pos.clear(); - -// return s_pos; - -// } - double Accelerator::get_time_aware_elements_info(std::vector& TAW_indices, std::vector& TAW_positions, unsigned int element_offset) const { // for longitudinal kick before RF cavities From d079806518405d52e4df7b7eb9cb9c76f2ec6d75 Mon Sep 17 00:00:00 2001 From: Vitor Date: Fri, 12 Sep 2025 00:30:18 -0300 Subject: [PATCH 20/20] last saves --- include/trackcpp/auxiliary.h | 16 +++++------ include/trackcpp/tracking.h | 52 +++++++++++++++++++++--------------- src/accelerator.cpp | 30 +++++++++------------ src/auxiliary.cpp | 2 ++ src/flat_file.cpp | 4 +-- src/tests.cpp | 6 ++--- src/tracking.cpp | 18 ++++++++----- tests/test-kickmap.cpp | 2 +- 8 files changed, 70 insertions(+), 60 deletions(-) diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 6e368c4..99378e3 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -38,6 +38,7 @@ class PassMethodsClass { 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 std::vector TAW_pms; PassMethodsClass() { passmethods.push_back("identity_pass"); passmethods.push_back("drift_pass"); @@ -47,20 +48,19 @@ class PassMethodsClass { passmethods.push_back("cavity_pass"); passmethods.push_back("thinquad_pass"); passmethods.push_back("thinsext_pass"); - passmethods.push_back("kicktable_pass"); + passmethods.push_back("kickmap_pass"); passmethods.push_back("matrix_pass"); passmethods.push_back("drift_g2l_pass"); } int size() const { return passmethods.size(); } + std::string operator[](const int i) const { return passmethods[i]; } + bool is_time_aware_pm(const int i) const { - switch (i) - { - case pm_cavity_pass: - return true; - default: - return false; + for (int pm: TAW_pms) { + if (i == pm) {return true;} } + return false; }; private: std::vector passmethods; @@ -95,7 +95,7 @@ const std::vector pm_dict = { "cavity_pass", "thinquad_pass", "thinsext_pass", - "kicktable_pass", + "kickmap_pass", "matrix_pass", "drift_g2l_pass", }; diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index b57cf04..3f6896d 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -52,7 +52,7 @@ Pos linalg_solve6_posvec( template -Status::type track_elementpass ( +Status::type track_elementpass( const Accelerator& accelerator, const Element& el, // element through which to track particle Pos &orig_pos // initial electron coordinates @@ -101,7 +101,7 @@ Status::type track_elementpass ( } template -Status::type track_elementpass ( +Status::type track_elementpass( const Accelerator& accelerator, const Element& el, // element through which to track particle std::vector >& orig_pos // initial electron coordinates @@ -137,7 +137,7 @@ Status::type track_elementpass ( // status do tracking (see 'auxiliary.h') // template -Status::type track_linepass ( +Status::type track_linepass( const Accelerator& accelerator, Pos& orig_pos, const std::vector& indices, @@ -152,7 +152,10 @@ Status::type track_linepass ( Status::type status = Status::success; const std::vector& line = accelerator.lattice; int nr_elements = line.size(); + + // For longitudinal RF kick double ddl = 0.0; + unsigned int TAW_pivot = 0; //pos.clear(); other functions assume pos is not clearedin linepass! pos.reserve(pos.size() + indices.size()); @@ -163,7 +166,6 @@ Status::type track_linepass ( for (unsigned int i=0; i<=nr_elements; ++i) indcs[i] = false; for (auto&& i: indices) if (i<=nr_elements) indcs[i] = true; - unsigned int TAW_pivot = 0; for(int i=0; i distr_gauss(0., 1.); static std::uniform_real_distribution distr_uniform(-sqrt(3.0), sqrt(3.0)); int choosen_distribution = Distributions::normal; +const std::vector PassMethodsClass::TAW_pms = { PassMethodsClass::pm_cavity_pass }; + void set_random_distribution(unsigned value){ if (value == Distributions::normal){ choosen_distribution = value; diff --git a/src/flat_file.cpp b/src/flat_file.cpp index a6ecc4b..9cd58d3 100644 --- a/src/flat_file.cpp +++ b/src/flat_file.cpp @@ -98,10 +98,10 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator) fp << std::setw(pw) << "fam_name" << e.fam_name << '\n'; fp << std::setw(pw) << "length" << e.length << '\n'; fp << std::setw(pw) << "pass_method" << pm_dict[e.pass_method] << '\n'; - if (pm_dict[e.pass_method].compare("kicktable_pass") == 0) { + if (pm_dict[e.pass_method].compare("kickmap_pass") == 0) { unsigned int idx = e.kicktable_idx; const Kicktable& ktable = kicktable_list[idx]; - fp << std::setw(pw) << "kicktable_fname" << ktable.filename << '\n'; + fp << std::setw(pw) << "kickmap_fname" << ktable.filename << '\n'; } if (e.nr_steps != 1) { fp.unsetf(std::ios_base::showpos); diff --git a/src/tests.cpp b/src/tests.cpp index 7cdc568..e213a35 100644 --- a/src/tests.cpp +++ b/src/tests.cpp @@ -39,7 +39,7 @@ int test_linepass(const Accelerator& accelerator) { // for longitudinal kick before RF cavities std::vector TAW_positions; std::vector TAW_indices; - double acc_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 37 + double acc_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, element_offset); Status::type status = track_linepass( accelerator, pos, true, element_offset, new_pos, lost_plane, acc_length, TAW_indices, TAW_positions); std::cout << "status: " << string_error_messages[status] << std::endl; @@ -73,7 +73,7 @@ int test_linepass_tpsa(const Accelerator& accelerator, const std::vector TAW_positions; std::vector TAW_indices; - double acc_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 71 + double acc_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, element_offset); track_linepass( accelerator, tpsa, false, element_offset, new_tpsa, lost_plane, acc_length, TAW_indices, TAW_positions); for(unsigned int i=0; i TAW_positions; std::vector TAW_indices; - double acc_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 414 + double acc_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, element_offset); orig_pos.rx = 0.0001; orig_pos.px = 0.0001; diff --git a/src/tracking.cpp b/src/tracking.cpp index 9c4d732..ff97917 100644 --- a/src/tracking.cpp +++ b/src/tracking.cpp @@ -66,11 +66,14 @@ Status::type track_findm66 (Accelerator& accelerator, map.de = Tpsa<6,1>(fp.de, 4); map.dl = Tpsa<6,1>(fp.dl, 5); tm.clear(); tm.reserve(indices.size()); + + // For longitudinal RF kick unsigned int TAW_pivot = 0; std::vector TAW_positions; std::vector TAW_indices; double line_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); double ddl = 0; + for(unsigned int i=0; i& fixed_point_guess) { const std::vector& the_ring = accelerator.lattice; + unsigned int element_offset = 0; double delta = 1e-9; // [m],[rad],[dE/E] double tolerance = 2.22044604925e-14; @@ -168,14 +175,13 @@ Status::type track_findorbit6( // for longitudinal kick before RF cavities std::vector TAW_positions; std::vector TAW_indices; - double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 182 + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, element_offset); int nr_iter = 0; while ((get_max(dco) > tolerance) and (nr_iter <= max_nr_iters)) { co = co + D; Pos Ri = co[6]; std::vector > co2; - unsigned int element_offset = 0; Plane::type lost_plane; Status::type status = Status::success; @@ -232,7 +238,6 @@ Status::type track_findorbit6( // propagates fixed point throught the_ring closed_orbit.clear(); - unsigned int element_offset = 0; Plane::type lost_plane; track_linepass( accelerator, co[6], true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions @@ -248,6 +253,7 @@ Status::type track_findorbit4( const Pos& fixed_point_guess) { const std::vector& the_ring = accelerator.lattice; + unsigned int element_offset = 0; double delta = 1e-9; // [m],[rad],[dE/E] double tolerance = 2.22044604925e-14; @@ -271,14 +277,13 @@ Status::type track_findorbit4( // for longitudinal kick before RF cavities std::vector TAW_positions; std::vector TAW_indices; - double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 285 + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, element_offset); int nr_iter = 0; while ((get_max(dco) > tolerance) and (nr_iter <= max_nr_iters)) { co = co + D; Pos Ri = co[6]; std::vector > co2; - unsigned int element_offset = 0; Plane::type lost_plane; Status::type status = Status::success; status = (Status::type) ((int) status | (int) track_linepass( @@ -322,7 +327,6 @@ Status::type track_findorbit4( // propagates fixed point throught the_ring closed_orbit.clear(); - unsigned int element_offset = 0; Plane::type lost_plane; track_linepass( accelerator, co[6], true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions diff --git a/tests/test-kickmap.cpp b/tests/test-kickmap.cpp index 78681f0..ccb8110 100644 --- a/tests/test-kickmap.cpp +++ b/tests/test-kickmap.cpp @@ -68,7 +68,7 @@ int main() { // for longitudinal kick before RF cavities std::vector TAW_positions; std::vector TAW_indices; - double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, 0); // 0 -> element_offset = 0 #line 66 + double accelerator_length = accelerator.get_time_aware_elements_info(TAW_indices, TAW_positions, element_offset); status = track_linepass( accelerator, fp, true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions);