Cav studies

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<unsigned int>& TAW_indices, std::vector<double>& TAW_positions, unsigned int element_offset = 0) const;
 };
 
 #endif

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<int> 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<std::string> passmethods;
 };
@@ -86,7 +95,7 @@ const std::vector<std::string> pm_dict = {
     "cavity_pass",
     "thinquad_pass",
     "thinsext_pass",
-    "kicktable_pass",
+    "kickmap_pass",
     "matrix_pass",
     "drift_g2l_pass",
 };

include/trackcpp/tracking.h

@@ -52,7 +52,7 @@ Pos<double>  linalg_solve6_posvec(
 
 
 template <typename T>
-Status::type track_elementpass (
+Status::type track_elementpass(
     const Accelerator& accelerator,
     const Element& el, // element through which to track particle
     Pos<T> &orig_pos // initial electron coordinates
@@ -101,7 +101,7 @@ Status::type track_elementpass (
 }
 
 template <typename T>
-Status::type track_elementpass (
+Status::type track_elementpass(
     const Accelerator& accelerator,
     const Element& el, // element through which to track particle
     std::vector<Pos<T> >& orig_pos // initial electron coordinates
@@ -137,19 +137,26 @@ Status::type track_elementpass (
 //      status do tracking (see 'auxiliary.h')
 //
 template <typename T>
-Status::type track_linepass (
+Status::type track_linepass(
     const Accelerator& accelerator,
     Pos<T>& orig_pos,
     const std::vector<unsigned int>& indices,
     unsigned int& element_offset,
     std::vector<Pos<T> >& pos,
-    Plane::type& lost_plane
+    Plane::type& lost_plane,
+    const double line_length,
+    const std::vector<unsigned int>& time_aware_element_indices,
+    const std::vector<double>& time_aware_element_positions
 ) {
 
     Status::type status = Status::success;
     const std::vector<Element>& 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 <typename T>
-Status::type track_linepass (
+Status::type track_linepass(
     const Accelerator& accelerator,
     Pos<T>& orig_pos,
     const bool trajectory,
     unsigned int& element_offset,
     std::vector<Pos<T> >& pos,
-    Plane::type& lost_plane
+    Plane::type& lost_plane,
+    const double line_length,
+    const std::vector<unsigned int>& time_aware_element_indices,
+    const std::vector<double>& time_aware_element_positions
 ) {
     std::vector<unsigned int> 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,15 +293,18 @@ Status::type track_linepass (
 //      status do tracking (see 'auxiliary.h')
 //
 template <typename T>
-Status::type track_linepass (
+Status::type track_linepass(
     const Accelerator& accelerator,
     std::vector<Pos<T>> &orig_pos,
     const std::vector<unsigned int >& indices,
     const unsigned int element_offset,
     std::vector<Pos<T>> &pos,
     std::vector<unsigned int >& lost_plane,
     std::vector<bool>& lost_flag,
-    std::vector<int>& lost_element
+    std::vector<int>& lost_element,
+    const double line_length,
+    const std::vector<unsigned int>& time_aware_element_indices,
+    const std::vector<double>& 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 <typename T>
-Status::type track_ringpass (
+Status::type track_ringpass(
     const Accelerator& accelerator,
     Pos<T>& orig_pos,
     const unsigned int nr_turns,
@@ -369,6 +410,11 @@ Status::type track_ringpass (
     Status::type status  = Status::success;
     std::vector<Pos<T> > final_pos;
 
+    // for longitudinal kick before RF cavities
+    std::vector<double> TAW_positions;
+    std::vector<unsigned int> 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<nr_turns; ++lost_turn) {
@@ -382,7 +428,10 @@ Status::type track_ringpass (
             false,
             element_offset,
             final_pos,
-            lost_plane
+            lost_plane,
+            accelerator_length,
+            TAW_indices,
+            TAW_positions
         )) != Status::success) {
 
             // fill last of vector with nans
@@ -401,7 +450,7 @@ Status::type track_ringpass (
 
 
 template <typename T>
-Status::type track_ringpass (
+Status::type track_ringpass(
     const Accelerator& accelerator,
     std::vector<Pos<T> > &orig_pos,
     const unsigned int nr_turns,

python_package/interface.cpp

@@ -52,6 +52,9 @@ Status::type track_linepass_wrapper(
 
     std::vector<Pos<double> > post;
     std::vector<Pos<double> > orig_post;
+    std::vector<double> TAW_positions;
+    std::vector<unsigned int> 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<ni2; ++i){
@@ -69,8 +72,11 @@ Status::type track_linepass_wrapper(
         post,
         args.lost_plane,
         args.lost_flag,
-        args.lost_element
-    );
+        args.lost_element,
+        accelerator_length,
+        TAW_indices,
+        TAW_positions
+);
     for (unsigned int i=0; i<post.size(); ++i){
         pos[0*n2 + i] = post[i].rx; pos[1*n2 + i] = post[i].px;
         pos[2*n2 + i] = post[i].ry; pos[3*n2 + i] = post[i].py;

src/accelerator.cpp

@@ -52,3 +52,41 @@ 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<unsigned int>& TAW_indices, std::vector<double>& 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<double> temp_positions;
+  std::vector<unsigned int> 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;
+
+}

src/auxiliary.cpp

@@ -22,6 +22,8 @@ static std::normal_distribution<double>  distr_gauss(0., 1.);
 static std::uniform_real_distribution<double> distr_uniform(-sqrt(3.0), sqrt(3.0));
 int choosen_distribution = Distributions::normal;
 
+const std::vector<int> PassMethodsClass::TAW_pms = { PassMethodsClass::pm_cavity_pass };
+
 void set_random_distribution(unsigned value){
   if (value == Distributions::normal){
     choosen_distribution = value;

src/commands.cpp

@@ -1084,7 +1084,13 @@ int cmd_track_linepass(const std::vector<std::string>& args) {
   std::vector<Pos<double>> 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<double> TAW_positions;
+  std::vector<unsigned int> 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");

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);

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<kicktable_list.size(); ++i) {
     if (kicktable_list[i].filename == filename) {

src/optics.cpp

@@ -81,12 +81,17 @@ Status::type calc_twiss(Accelerator& accelerator,
   std::vector<Pos<double>> closed_orbit;
   Plane::type lost_plane;
   unsigned int element_offset = 0;
+
+  // for longitudinal kick before RF cavities
+  std::vector<double> TAW_positions;
+  std::vector<unsigned int> 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 <double, std::milli> (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;
   }