Merge pull request #28 from lnls-fac/fix-eta-calc

Fix calculation of dispersion function.

Author: fernandohds564

.gitignore

@@ -34,3 +34,4 @@
 
 Debug/
 build/
+.vscode/

include/trackcpp/optics.h

@@ -48,8 +48,7 @@ Status::type calc_twiss(const Accelerator& accelerator,
                         const Pos<double>& fixed_point,
                         Matrix& m66,
                         std::vector<Twiss>& twiss,
-                        Twiss twiss0 = Twiss(),
-                        bool closed_flag = false);
+                        Twiss twiss0 = Twiss());
 
 
 #endif

include/trackcpp/tracking.h

@@ -26,11 +26,11 @@
 #include <limits>
 #include <cmath>
 
-Status::type track_findm66     (const Accelerator& accelerator, std::vector<Pos<double> >& closed_orbit, std::vector<Matrix>& tm, Matrix& m66, Pos<double>& v0);
-Status::type track_findorbit4  (const Accelerator& accelerator, std::vector<Pos<double> >& closed_orbit, const Pos<double>& fixed_point_guess = Pos<double>(0));
-Status::type track_findorbit6  (const Accelerator& accelerator, std::vector<Pos<double> >& closed_orbit, const Pos<double>& fixed_point_guess = Pos<double>(0));
-Pos<double>  linalg_solve4_posvec (const std::vector<Pos<double> >& M, const Pos<double>& b);
-Pos<double>  linalg_solve6_posvec (const std::vector<Pos<double> >& M, const Pos<double>& b);
+Status::type track_findm66(const Accelerator& accelerator, const Pos<double>& fixed_point, std::vector<Matrix>& tm, Matrix& m66, Pos<double>& v0);
+Status::type track_findorbit4(const Accelerator& accelerator, std::vector<Pos<double> >& closed_orbit, const Pos<double>& fixed_point_guess = Pos<double>(0));
+Status::type track_findorbit6(const Accelerator& accelerator, std::vector<Pos<double> >& closed_orbit, const Pos<double>& fixed_point_guess = Pos<double>(0));
+Pos<double>  linalg_solve4_posvec(const std::vector<Pos<double> >& M, const Pos<double>& b);
+Pos<double>  linalg_solve6_posvec(const std::vector<Pos<double> >& M, const Pos<double>& b);
 
 template <typename T>
 Status::type track_elementpass (

src/optics.cpp

@@ -64,8 +64,7 @@ Status::type calc_twiss(const Accelerator& accelerator,
                         const Pos<double>& fixed_point,
                         Matrix& m66,
                         std::vector<Twiss>& twiss,
-                        Twiss twiss0,
-                        bool closed_flag) {
+                        Twiss twiss0) {
 
 #ifdef TIMEIT
   auto start = std::chrono::steady_clock::now();
@@ -84,14 +83,7 @@ Status::type calc_twiss(const Accelerator& accelerator,
   unsigned int element_offset = 0;
   Status::type status = track_linepass(accelerator, fp, closed_orbit, element_offset, lost_plane, true);
   if (status != Status::success) return status;
-  if (not closed_flag) closed_orbit.pop_back();
 
-  // std::vector<Pos<double>> closed_orbit0;
-  // if (not accelerator.cavity_on) {
-  //   Pos<double> fp = fixed_point; fp.de = 0;
-  //   Status::type status = track_linepass(accelerator, fp, closed_orbit0, element_offset, lost_plane, true);
-  //   if (status != Status::success) return status;
-  // }
 
 #ifdef TIMEIT
   end = std::chrono::steady_clock::now(); diff = end - start;
@@ -102,22 +94,10 @@ Status::type calc_twiss(const Accelerator& accelerator,
   start = std::chrono::steady_clock::now();
 #endif
 
-  // finds accumulated transfer matrices
-
-  // std::vector<Matrix> atm0;
-  // if (not accelerator.cavity_on) {
-  //   status = track_findm66 (accelerator, closed_orbit0, atm0, m66);
-  //   if (status != Status::success) return status;
-  //   if (closed_flag) atm.push_back(atm.back());
-  // }
-
   std::vector<Matrix> atm;
   Pos<double> v0;
-  status = track_findm66 (accelerator, closed_orbit, atm, m66, v0);
+  status = track_findm66(accelerator, closed_orbit[0], atm, m66, v0);
   if (status != Status::success) return status;
-  if (closed_flag) atm.push_back(atm.back());
-
-
 
 #ifdef TIMEIT
   end = std::chrono::steady_clock::now(); diff = end - start;
@@ -128,6 +108,11 @@ Status::type calc_twiss(const Accelerator& accelerator,
   start = std::chrono::steady_clock::now();
 #endif
 
+  const double dpp = 1e-8;
+  Pos<double> fpp = fixed_point;
+  fpp.de += dpp;
+  std::vector<Pos<double>> codp;
+
   // initial twiss parameters
   twiss.clear(); twiss.reserve(atm.size());
   if (twiss0.isundef()) { // case of periodic solution
@@ -141,15 +126,32 @@ Status::type calc_twiss(const Accelerator& accelerator,
     twiss0.betay  =  m66[2][3]/sin_muy;
     // --- closed orbit
     twiss0.co = closed_orbit[0];
-    // --- dispersion function based on eta = (1 - M)^(-1) D
-    Vector Dx({m66[0][4], m66[1][4]});
-    Vector Dy({m66[2][4], m66[3][4]});
-    Matrix eye2({{1,0},{0,1}});
-    Matrix mx; m66.getM(mx, 2, 2, 0, 0); mx.linear_combination(1.0,eye2,-1.0,mx); mx.inverse();
-    Matrix my; m66.getM(my, 2, 2, 2, 2); my.linear_combination(1.0,eye2,-1.0,my); my.inverse();
-    twiss0.etax.multiplication(mx, Dx);
-    twiss0.etay.multiplication(my, Dy);
+
+    // // --- dispersion function based on eta = (1 - M)^(-1) D
+    // Vector Dx({m66[0][4], m66[1][4]});
+    // Vector Dy({m66[2][4], m66[3][4]});
+    // Matrix eye2({{1,0},{0,1}});
+    // Matrix mx; m66.getM(mx, 2, 2, 0, 0); mx.linear_combination(1.0,eye2,-1.0,mx); mx.inverse();
+    // Matrix my; m66.getM(my, 2, 2, 2, 2); my.linear_combination(1.0,eye2,-1.0,my); my.inverse();
+    // twiss0.etax.multiplication(mx, Dx);
+    // twiss0.etay.multiplication(my, Dy);
+
+    // Dispersion Function based on tracking:
+    Status::type status = track_findorbit4(accelerator, codp, fpp);
+    if (status != Status::success) return status;
+    twiss0.etax[0] = (codp[0].rx - closed_orbit[0].rx) / dpp;
+    twiss0.etax[1] = (codp[0].px - closed_orbit[0].px) / dpp;
+    twiss0.etay[0] = (codp[0].ry - closed_orbit[0].ry) / dpp;
+    twiss0.etay[1] = (codp[0].py - closed_orbit[0].py) / dpp;
+  } else {
+    fpp.rx += twiss0.etax[0] * dpp;
+    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);
+    if (status != Status::success) return status;
   }
+
   twiss.push_back(twiss0);
 
 #ifdef TIMEIT
@@ -178,20 +180,25 @@ Status::type calc_twiss(const Accelerator& accelerator,
     // --- closed orbit
     tw.co = closed_orbit[i];
 
-    // --- dispersion function
-    Matrix t1(atm[i-1]); t1.inverse();
-    Matrix T; T.multiplication(atm[i],t1);
-    Vector Dx({T[0][4], T[1][4]});
-    Vector Dy({T[2][4], T[3][4]});
-    Matrix mx; T.getMx(mx);
-    Matrix my; T.getMy(my);
-    tw.etax.multiplication(mx, twiss[i-1].etax);
-    tw.etay.multiplication(my, twiss[i-1].etay);
-    tw.etax = Dx + tw.etax.multiplication(mx, twiss[i-1].etax);
-    tw.etay = Dy + tw.etay.multiplication(my, twiss[i-1].etay);
+    // // --- dispersion function based on propagation
+    // Matrix t1(atm[i-1]); t1.inverse();
+    // Matrix T; T.multiplication(atm[i],t1);
+    // Vector Dx({T[0][4], T[1][4]});
+    // Vector Dy({T[2][4], T[3][4]});
+    // Matrix mx; T.getMx(mx);
+    // Matrix my; T.getMy(my);
+    // tw.etax.multiplication(mx, twiss[i-1].etax);
+    // tw.etay.multiplication(my, twiss[i-1].etay);
+    // tw.etax = Dx + tw.etax.multiplication(mx, twiss[i-1].etax);
+    // tw.etay = Dy + tw.etay.multiplication(my, twiss[i-1].etay);
+
+    // Dispersion Function based on tracking
+    tw.etax[0] = (codp[i].rx - closed_orbit[i].rx) / dpp;
+    tw.etax[1] = (codp[i].px - closed_orbit[i].px) / dpp;
+    tw.etay[0] = (codp[i].ry - closed_orbit[i].ry) / dpp;
+    tw.etay[1] = (codp[i].py - closed_orbit[i].py) / dpp;
 
     twiss.push_back(tw);
-
   }
 
   // unwraps betatron phases

src/tracking.cpp

@@ -24,43 +24,39 @@
 // returns a vector with 6-d transfer matrices, one for each element
 //
 // inputs:
-//    accelerator:   structure representing the accelerator
-//    closed_orbit:  Pos vector representing calculated closed orbit.
+//    accelerator:  Structure representing the accelerator
+//    fixed_point:  Pos representing calculated fixed_point.
 //
 // outputs:
 //    tm:     vector of Matrix6 elements. Each component represents the accumulated
-//            transfer matrix from the start of the lattice to the entrance of that element.
+//            transfer matrix from the start of the lattice to the entrance of that element. The last element is the m66 of the line.
 //    m66:    one-turn transfer matrix
 //
 //    v0:     const term of final map
 //
 //    RETURN:      status do tracking (see 'auxiliary.h')
 
 Status::type track_findm66 (const Accelerator& accelerator,
-                            std::vector<Pos<double> >& closed_orbit,
+                            const Pos<double>& fixed_point,
                             std::vector<Matrix>& tm,
                             Matrix& m66,
                             Pos<double>& v0) {
 
   Status::type status  = Status::success;
   const std::vector<Element>& lattice = accelerator.lattice;
 
-  std::vector<double> row0 = {0,0,0,0,0,0};
+  Pos<double> fp = fixed_point;
 
   // case no closed_orbit has been defined
-  if (closed_orbit.size() != lattice.size()) {
-    closed_orbit.clear();
-    for(unsigned int i=0; i<lattice.size(); ++i) {
-      closed_orbit.push_back(Pos<double>(0,0,0,0,0,0));
-    }
-  }
+  if (std::isnan(fp.rx)) fp = Pos<double>(0,0,0,0,0,0);
+
 
   Pos<Tpsa<6,1> > map;
-  map.rx = Tpsa<6,1>(closed_orbit[0].rx, 0); map.px = Tpsa<6,1>(closed_orbit[0].px, 1);
-  map.ry = Tpsa<6,1>(closed_orbit[0].ry, 2); map.py = Tpsa<6,1>(closed_orbit[0].py, 3);
-  map.de = Tpsa<6,1>(closed_orbit[0].de, 4); map.dl = Tpsa<6,1>(closed_orbit[0].dl, 5);
+  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.resize(lattice.size(), Matrix(6));
+  tm.clear(); tm.resize(lattice.size()+1, Matrix(6));
   for(unsigned int i=0; i<lattice.size(); ++i) {
     Matrix& m = tm[i];
     m[0][0] = map.rx.c[1]; m[0][1] = map.rx.c[2]; m[0][2] = map.rx.c[3]; m[0][3] = map.rx.c[4]; m[0][4] = map.rx.c[5]; m[0][5] = map.rx.c[6];
@@ -82,6 +78,8 @@ Status::type track_findm66 (const Accelerator& accelerator,
   m[4][0] = map.de.c[1]; m[4][1] = map.de.c[2]; m[4][2] = map.de.c[3]; m[4][3] = map.de.c[4]; m[4][4] = map.de.c[5]; m[4][5] = map.de.c[6];
   m[5][0] = map.dl.c[1]; m[5][1] = map.dl.c[2]; m[5][2] = map.dl.c[3]; m[5][3] = map.dl.c[4]; m[5][4] = map.dl.c[5]; m[5][5] = map.dl.c[6];
 
+  tm[lattice.size()] = m66;
+
   // 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];