diff --git a/include/trackcpp/accelerator.h b/include/trackcpp/accelerator.h index 0653507..9ad28d3 100644 --- a/include/trackcpp/accelerator.h +++ b/include/trackcpp/accelerator.h @@ -40,6 +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; }; #endif diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h index 3cdfdc5..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,12 +48,20 @@ 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 { + for (int pm: TAW_pms) { + if (i == pm) {return true;} + } + return false; + }; private: std::vector passmethods; }; @@ -86,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 a05c7a4..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,19 +137,26 @@ 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, unsigned int& element_offset, std::vector >& pos, - Plane::type& lost_plane + 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; 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()); @@ -167,7 +174,24 @@ Status::type track_linepass ( // stores trajectory at entrance of each element if (indcs[i]) pos.push_back(orig_pos); + // For longitudinal RF kick + 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] / line_length; + if (TAW_pivot < time_aware_element_indices.size() - 1) {TAW_pivot++;} + } + + // tracking itself status = track_elementpass(accelerator, element, orig_pos); + + //! Superseded by lines 178-181 + // 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; + // } + + // checking particle loss lost_plane = check_particle_loss(accelerator, element, orig_pos); if (lost_plane != Plane::no_plane) status = Status::particle_lost; @@ -211,13 +235,16 @@ Status::type track_linepass ( // status do tracking (see 'auxiliary.h') // template -Status::type track_linepass ( +Status::type track_linepass( const Accelerator& accelerator, Pos& orig_pos, const bool trajectory, unsigned int& element_offset, std::vector >& pos, - Plane::type& lost_plane + Plane::type& lost_plane, + 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(); @@ -228,13 +255,16 @@ Status::type track_linepass ( indices.push_back(nr_elements); } - return track_linepass ( + return track_linepass( accelerator, orig_pos, indices, element_offset, pos, - lost_plane + lost_plane, + line_length, + time_aware_element_indices, + time_aware_element_positions ); } @@ -263,7 +293,7 @@ Status::type track_linepass ( // status do tracking (see 'auxiliary.h') // template -Status::type track_linepass ( +Status::type track_linepass( const Accelerator& accelerator, std::vector> &orig_pos, const std::vector& indices, @@ -271,7 +301,10 @@ Status::type track_linepass ( std::vector> &pos, std::vector& lost_plane, std::vector& lost_flag, - std::vector& lost_element + std::vector& lost_element, + const double line_length, + const std::vector& time_aware_element_indices, + const std::vector& time_aware_element_positions ) { int nr_elements = accelerator.lattice.size(); @@ -289,7 +322,15 @@ 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, + line_length, + time_aware_element_indices, + time_aware_element_positions ); if (status2 != Status::success){ @@ -355,7 +396,7 @@ Status::type track_linepass ( // RETURN: // status do tracking (see 'auxiliary.h') template -Status::type track_ringpass ( +Status::type track_ringpass( const Accelerator& accelerator, Pos& orig_pos, const unsigned int nr_turns, @@ -369,6 +410,11 @@ Status::type track_ringpass ( Status::type status = Status::success; std::vector > final_pos; + // 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, element_offset); + if (turn_by_turn) pos.reserve(nr_turns+1); for(lost_turn=0; lost_turn -Status::type track_ringpass ( +Status::type track_ringpass( const Accelerator& accelerator, std::vector > &orig_pos, const unsigned int nr_turns, diff --git a/python_package/interface.cpp b/python_package/interface.cpp index 383c7fe..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& 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; + double acclen = 0.0; + for (size_t i = 0; i < nr_elements; i++) { + const Element& element = this->lattice[element_offset]; + acclen += element.length; + if (PMClass.is_time_aware_pm(element.pass_method)) { + temp_indices.push_back(i); + temp_positions.push_back(s_pos + element.length/2); + s_pos = 0.0 + element.length/2; + } + else { + s_pos += element.length; + } + element_offset = (element_offset + 1) % nr_elements; + } + + // In case no element is "time aware" + if (temp_indices.size() < 1) { + TAW_indices.push_back(-1); + TAW_positions.push_back(0.0); + } else { + temp_positions[0] += s_pos; + } + + return acclen; + +} diff --git a/src/auxiliary.cpp b/src/auxiliary.cpp index 78ddeff..2b08e34 100644 --- a/src/auxiliary.cpp +++ b/src/auxiliary.cpp @@ -22,6 +22,8 @@ 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::TAW_pms = { 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..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); + + // 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/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/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> closed_orbit; 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); + Status::type status = track_linepass( - accelerator, fp, true, element_offset, closed_orbit, lost_plane + accelerator, fp, true, element_offset, closed_orbit, lost_plane, 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; @@ -151,7 +156,7 @@ Status::type calc_twiss(Accelerator& accelerator, 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, 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 8ed7ead..e213a35 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, element_offset); Status::type status = track_linepass( - accelerator, pos, true, element_offset, new_pos, lost_plane - ); + 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, element_offset); track_linepass( - accelerator, tpsa, false, element_offset, new_tpsa, lost_plane - ); + 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, element_offset); + orig_pos.rx = 0.0001; orig_pos.px = 0.0001; @@ -418,7 +429,10 @@ int test_linepass2() { trajectory, element_offset, pos, - lost_plane + lost_plane, + 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 f751390..ff97917 100644 --- a/src/tracking.cpp +++ b/src/tracking.cpp @@ -45,7 +45,6 @@ Status::type track_findm66 (Accelerator& accelerator, Status::type status = Status::success; const std::vector& lattice = accelerator.lattice; - Pos fp = fixed_point; const int radsts = accelerator.radiation_on; @@ -67,6 +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; @@ -146,12 +162,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,39 +170,41 @@ 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); + // 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, 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( - accelerator, co[0], false, element_offset, co2, lost_plane + 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, 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, 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, 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, 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, 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, co[6], false, element_offset, co2, lost_plane, accelerator_length, TAW_indices, TAW_positions )); if (status != Status::success) { @@ -208,7 +220,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; @@ -226,10 +238,9 @@ 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, co[6], true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions ); accelerator.radiation_on = radsts; return Status::success; @@ -242,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; @@ -251,6 +263,7 @@ Status::type track_findorbit4( if (radsts == RadiationState::full){ accelerator.radiation_on = RadiationState::damping; } + // temporary vectors and matrices // std::vector > co(7,0); // no initial guess std::vector > co(7,fixed_point_guess); @@ -261,28 +274,32 @@ 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, 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( - accelerator, co[0], false, element_offset, co2, lost_plane + 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, 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, 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, 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, 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; @@ -310,10 +327,9 @@ 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, 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 efbb3ff..ccb8110 100644 --- a/tests/test-kickmap.cpp +++ b/tests/test-kickmap.cpp @@ -65,9 +65,13 @@ 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, element_offset); + status = track_linepass( - accelerator, fp, true, element_offset, closed_orbit, lost_plane - ); + accelerator, fp, true, element_offset, closed_orbit, lost_plane, accelerator_length, TAW_indices, TAW_positions); if (status != Status::success) return status; std::vector atm; @@ -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 - ); + 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; } @@ -137,7 +140,7 @@ int main() { // Plane::type lost_plane; // unsigned int element_offset = 0; // Status::type status1 = track_linepass( - // sirius, fp, true, element_offset, traj, lost_plane + // 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;