Longitudinal Phase-Space Rework

include/trackcpp/accelerator.h

@@ -34,12 +34,15 @@ class Accelerator {
   int                     harmonic_number = 0;
   std::vector<Element>    lattice;
   std::string             lattice_version = "";
+  mutable std::vector<unsigned int> time_aware_indices;
+  mutable std::vector<double> 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

include/trackcpp/auxiliary.h

@@ -23,6 +23,7 @@
 #include <iostream>
 #include <cmath>
 #include <random>
+#include <algorithm>
 
 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<int> 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<std::string> passmethods;
 };

include/trackcpp/passmethods.hpp

@@ -32,6 +32,7 @@
 #include "tpsa.h"
 #include "linalg.h"
 #include <cmath>
+#include <algorithm>
 
 template <typename T> inline T SQR(const T& X) { return X*X; }
 template <typename T> inline T POW3(const T& X) { return X*X*X; }
@@ -499,4 +500,17 @@ Status::type pm_matrix_pass(Pos<T> &pos, const Element &elem,
   return Status::success;
 }
 
+template <typename T>
+inline void adjust_path_length(
+  const Accelerator& accelerator,
+  unsigned int& element_index,
+  Pos<T>& 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

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<Pos<T> > 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_turn<nr_turns; ++lost_turn) {

python_package/interface.cpp

@@ -61,6 +61,10 @@ Status::type track_linepass_wrapper(
             orig_pos[4*ni2 + i], orig_pos[5*ni2 + i]
         );
     }
+
+    // adjust dl to keep the arrival-time in sync with wall clock
+    accelerator.update_time_aware_info();
+
     Status::type status = track_linepass(
         accelerator,
         orig_post,

python_package/trackcpp/VERSION

@@ -1 +1 @@
-4.10.8
+5.0.0

src/accelerator.cpp

@@ -52,3 +52,47 @@ std::ostream& operator<< (std::ostream &out, const Accelerator& a) {
   out << std::endl << "lattice_version: " << a.lattice_version;
   return out;
 }
+
+void Accelerator::update_time_aware_info(void) const {
+
+  // adjust dl to keep the arrival-time in sync with wall clock
+  this->time_aware_indices.clear();
+  this->time_aware_dl_kicks.clear();
+  std::vector<double> 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; i<this->time_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);
+  }
+
+}

src/auxiliary.cpp

@@ -22,6 +22,10 @@ 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::time_aware_passmethods = {
+  PassMethodsClass::pm_cavity_pass,
+};
+
 void set_random_distribution(unsigned value){
   if (value == Distributions::normal){
     choosen_distribution = value;

src/commands.cpp

@@ -1084,7 +1084,10 @@ 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);
+  // 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");

src/optics.cpp

@@ -81,6 +81,10 @@ Status::type calc_twiss(Accelerator& accelerator,
   std::vector<Pos<double>> 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
   );

src/tests.cpp

@@ -36,6 +36,10 @@ int test_linepass(const Accelerator& accelerator) {
   std::vector<Pos<> > 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<Element
   std::vector<Pos<Tpsa<6,order> > > 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;
 

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<lattice.size(); ++i) {
+    // 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, i, map);
     if (indcs[i]){
       Matrix m (6);
       m[0][0] = map.rx.c[1]; m[0][1] = map.rx.c[2]; m[0][2] = map.rx.c[3];
@@ -137,7 +143,6 @@ Status::type track_findorbit6(
     const Pos<double>& fixed_point_guess) {
 
   const std::vector<Element>& 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<int>    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<Pos<double> > co(7,0);
@@ -160,10 +159,11 @@ Status::type track_findorbit6(
   std::vector<Pos<double> > D(7,0);
   std::vector<Pos<double> > M(6,0);
   Pos<double> dco(1.0,1.0,1.0,1.0,1.0,1.0);
-  Pos<double> 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<double> b = Rf - Ri - theta;
+    Pos<double> b = Rf - Ri;
     std::vector<Pos<double> > 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<double>& fixed_point_guess) {
 
   const std::vector<Element>& 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<double> 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;

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