diff --git a/include/trackcpp/accelerator.h b/include/trackcpp/accelerator.h index 0653507..9116a61 100644 --- a/include/trackcpp/accelerator.h +++ b/include/trackcpp/accelerator.h @@ -34,12 +34,15 @@ class Accelerator { int harmonic_number = 0; std::vector lattice; std::string lattice_version = ""; + mutable std::vector time_aware_indices; + mutable std::vector time_aware_dl_kicks; bool operator==(const Accelerator& o) const; 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; friend std::ostream& operator<< (std::ostream &out, const Accelerator& a); + void update_time_aware_info(void) const; }; #endif diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 3cdfdc5..35e2956 100644 --- a/include/trackcpp/auxiliary.h +++ b/include/trackcpp/auxiliary.h @@ -23,6 +23,7 @@ #include #include #include +#include class PassMethodsClass { public: @@ -38,6 +39,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 time_aware_passmethods; PassMethodsClass() { passmethods.push_back("identity_pass"); passmethods.push_back("drift_pass"); @@ -53,6 +55,9 @@ class PassMethodsClass { } 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 { + return std::find(time_aware_passmethods.begin(), time_aware_passmethods.end(), i) != time_aware_passmethods.end(); + }; private: std::vector passmethods; }; diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp index 9b315ea..0931064 100644 --- a/include/trackcpp/passmethods.hpp +++ b/include/trackcpp/passmethods.hpp @@ -32,6 +32,7 @@ #include "tpsa.h" #include "linalg.h" #include +#include template inline T SQR(const T& X) { return X*X; } template inline T POW3(const T& X) { return X*X*X; } @@ -499,4 +500,17 @@ Status::type pm_matrix_pass(Pos &pos, const Element &elem, return Status::success; } +template +inline void adjust_path_length( + const Accelerator& accelerator, + unsigned int& element_index, + Pos& pos +) { + auto it = std::find(accelerator.time_aware_indices.begin(), accelerator.time_aware_indices.end(), element_index); + if (it != accelerator.time_aware_indices.end()) { + auto idx = std::distance(accelerator.time_aware_indices.begin(), it); + pos.dl -= accelerator.time_aware_dl_kicks[idx]; + } +} + #endif diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h index a05c7a4..0542992 100644 --- a/include/trackcpp/tracking.h +++ b/include/trackcpp/tracking.h @@ -164,6 +164,11 @@ Status::type track_linepass ( // syntactic-sugar for read-only access to element object parameters const Element& element = line[element_offset]; + // adjust dl to keep the arrival-time in sync with wall clock + // note: for performance reasons, the adjustment of the path length + // is done only at the beginning of "time aware" elements (e.g. RF cavities) + adjust_path_length(accelerator, element_offset, orig_pos); + // stores trajectory at entrance of each element if (indcs[i]) pos.push_back(orig_pos); @@ -289,7 +294,12 @@ Status::type track_linepass ( unsigned int le = element_offset; status2 = track_linepass( - accelerator, orig_pos[i], indices, le, final_pos, lp + accelerator, + orig_pos[i], + indices, + le, + final_pos, + lp ); if (status2 != Status::success){ @@ -369,6 +379,9 @@ Status::type track_ringpass ( Status::type status = Status::success; std::vector > final_pos; + // adjust dl to keep the arrival-time in sync with wall clock + accelerator.update_time_aware_info(); + if (turn_by_turn) pos.reserve(nr_turns+1); for(lost_turn=0; lost_turntime_aware_indices.clear(); + this->time_aware_dl_kicks.clear(); + std::vector time_aware_displacements = {}; + + size_t nr_elements = this->lattice.size(); + PassMethodsClass PMClass = PassMethodsClass(); + + double s_pos = 0.0; + + for (size_t i = 0; i < nr_elements; i++) { + const Element& element = this->lattice[i]; + + if (PMClass.is_time_aware_pm(element.pass_method)) { + this->time_aware_indices.push_back(i); + s_pos += 0.5 * element.length; + time_aware_displacements.push_back(s_pos); + s_pos = 0.5 * element.length; + } + else { + s_pos += element.length; + } + } + + if (this->time_aware_indices.size() > 0) { + time_aware_displacements[0] += s_pos; + } + + //? NOTE : The diference between "line_length" and the sum of "time_aware_displacements" (cum_length) is on the order of + //? 1e-15 ~ 1e-16. The propagation of this tiny error affects the tracking. The following lines avoid losing precision. + double cum_length = std::accumulate(time_aware_displacements.begin(), time_aware_displacements.end(), 0.0); + double line_length = this->get_length(); + time_aware_displacements.back() += line_length - cum_length; + + double ddl = 0.0; + for (size_t i=0; itime_aware_indices.size(); i++) { + ddl = light_speed*this->harmonic_number/this->lattice[this->time_aware_indices[i]].frequency - line_length; + this->time_aware_dl_kicks.push_back(ddl * time_aware_displacements[i] / line_length); + } + +} diff --git a/src/auxiliary.cpp b/src/auxiliary.cpp index 78ddeff..ea3ec74 100644 --- a/src/auxiliary.cpp +++ b/src/auxiliary.cpp @@ -22,6 +22,10 @@ static std::normal_distribution 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::time_aware_passmethods = { + PassMethodsClass::pm_cavity_pass, +}; + void set_random_distribution(unsigned value){ if (value == Distributions::normal){ choosen_distribution = value; diff --git a/src/commands.cpp b/src/commands.cpp index d09f431..6df24f2 100644 --- a/src/commands.cpp +++ b/src/commands.cpp @@ -1084,7 +1084,10 @@ 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); + // adjust dl to keep the arrival-time in sync with wall clock + accelerator.update_time_aware_info(); + track_linepass(accelerator, pos, true, offset_element, pos_list, lost_plane + ); 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 416f850..f38c83f 100644 --- a/src/optics.cpp +++ b/src/optics.cpp @@ -81,6 +81,10 @@ Status::type calc_twiss(Accelerator& accelerator, std::vector> closed_orbit; Plane::type lost_plane; unsigned int element_offset = 0; + + // adjust dl to keep the arrival-time in sync with wall clock + accelerator.update_time_aware_info(); + Status::type status = track_linepass( accelerator, fp, true, element_offset, closed_orbit, lost_plane ); diff --git a/src/tests.cpp b/src/tests.cpp index 8ed7ead..a084365 100644 --- a/src/tests.cpp +++ b/src/tests.cpp @@ -36,6 +36,10 @@ int test_linepass(const Accelerator& accelerator) { std::vector > new_pos; unsigned int element_offset = 0; Plane::type lost_plane; + + // adjust dl to keep the arrival-time in sync with wall clock + accelerator.update_time_aware_info(); + Status::type status = track_linepass( accelerator, pos, true, element_offset, new_pos, lost_plane ); @@ -67,6 +71,10 @@ int test_linepass_tpsa(const Accelerator& accelerator, const std::vector > > new_tpsa; unsigned int element_offset = 0; Plane::type lost_plane; + + // adjust dl to keep the arrival-time in sync with wall clock + accelerator.update_time_aware_info(); + track_linepass( accelerator, tpsa, false, element_offset, new_tpsa, lost_plane ); @@ -409,6 +417,9 @@ int test_linepass2() { Plane::type lost_plane; bool trajectory = true; + // adjust dl to keep the arrival-time in sync with wall clock + accelerator.update_time_aware_info(); + orig_pos.rx = 0.0001; orig_pos.px = 0.0001; diff --git a/src/tracking.cpp b/src/tracking.cpp index f751390..957ee12 100644 --- a/src/tracking.cpp +++ b/src/tracking.cpp @@ -66,8 +66,14 @@ Status::type track_findm66 (Accelerator& accelerator, 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); + // adjust dl to keep the arrival-time in sync with wall clock + accelerator.update_time_aware_info(); + tm.clear(); tm.reserve(indices.size()); for(unsigned int i=0; i& fixed_point_guess) { 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; @@ -146,12 +151,6 @@ Status::type track_findorbit6( if (radsts == RadiationState::full){ accelerator.radiation_on = RadiationState::damping; } - // calcs longitudinal fixed point - double L0 = latt_findspos(the_ring, 1+the_ring.size()); - double T0 = L0 / light_speed; - std::vector 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,10 +159,11 @@ 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); + // adjust dl to keep the arrival-time in sync with wall clock + accelerator.update_time_aware_info(); + int nr_iter = 0; while ((get_max(dco) > tolerance) and (nr_iter <= max_nr_iters)) { co = co + D; @@ -208,7 +208,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; @@ -242,7 +242,6 @@ Status::type track_findorbit4( const Pos& fixed_point_guess) { 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; @@ -261,6 +260,9 @@ Status::type track_findorbit4( Pos theta(0.0,0.0,0.0,0.0,0.0,0.0); matrix6_set_identity_posvec(D, delta); + // adjust dl to keep the arrival-time in sync with wall clock + accelerator.update_time_aware_info(); + 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 efbb3ff..8f8eb1e 100644 --- a/tests/test-kickmap.cpp +++ b/tests/test-kickmap.cpp @@ -65,6 +65,9 @@ int main() { Plane::type lost_plane; unsigned int element_offset = 0; + // adjust dl to keep the arrival-time in sync with wall clock + accelerator.update_time_aware_info(); + status = track_linepass( accelerator, fp, true, element_offset, closed_orbit, lost_plane );