Add Twiss alpha to BeamRelevant and Docs (#3993)

2 files changed, 27 insertions, 15 deletions

diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index 1424e518a..005b1ffaa 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -2748,11 +2748,15 @@ Reduced Diagnostics
         :math:`\epsilon_z = \dfrac{1}{mc} \sqrt{\delta_z^2 \delta_{pz}^2 -
         \Big\langle (z-\langle z \rangle) (p_z-\langle p_z \rangle) \Big\rangle^2}`.
 
-        [19], [20]: beta function for the transverse directions (m)
+        [19], [20]: Twiss alpha for the transverse directions
+        :math:`\alpha_x = - \Big\langle (x-\langle x \rangle) (p_x-\langle p_x \rangle) \Big\rangle \Big/ \epsilon_x`,
+        :math:`\alpha_y = - \Big\langle (y-\langle y \rangle) (p_y-\langle p_y \rangle) \Big\rangle \Big/ \epsilon_y`.
+
+        [21], [22]: beta function for the transverse directions (m)
         :math:`\beta_x = \dfrac{{\delta_x}^2}{\epsilon_x}`,
         :math:`\beta_y = \dfrac{{\delta_y}^2}{\epsilon_y}`.
 
-        [21]: The charge of the beam (C).
+        [23]: The charge of the beam (C).
 
         For 2D-XZ,
         :math:`\langle y \rangle`,
diff --git a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp
index 4c2afede5..a99d476fa 100644
--- a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp
+++ b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp
@@ -47,9 +47,10 @@ BeamRelevant::BeamRelevant (std::string rd_name)
     // 10,11,12: rms px,py,pz
     //       13: rms gamma
     // 14,15,16: emittance x,y,z
-    //    17,18: beta-function x,y
-    //       19: charge
-    m_data.resize(20, 0.0_rt);
+    //    17,18: Twiss-alpha x,y
+    //    19,20: beta-function x,y
+    //       21: charge
+    m_data.resize(22, 0.0_rt);
 #elif (defined WARPX_DIM_XZ)
     //     0, 1: mean x,z
     //  2, 3, 4: mean px,py,pz
@@ -58,9 +59,10 @@ BeamRelevant::BeamRelevant (std::string rd_name)
     //  8, 9,10: rms px,py,pz
     //       11: rms gamma
     //    12,13: emittance x,z
-    //       14: beta-function x
-    //       15: charge
-    m_data.resize(16, 0.0_rt);
+    //       14: Twiss-alpha x
+    //       15: beta-function x
+    //       16: charge
+    m_data.resize(17, 0.0_rt);
 #elif (defined WARPX_DIM_1D_Z)
     //       0 : mean z
     //   1,2,3 : mean px,py,pz
@@ -102,6 +104,8 @@ BeamRelevant::BeamRelevant (std::string rd_name)
             ofs << "[" << c++ << "]emittance_x(m)";   ofs << m_sep;
             ofs << "[" << c++ << "]emittance_y(m)";   ofs << m_sep;
             ofs << "[" << c++ << "]emittance_z(m)";   ofs << m_sep;
+            ofs << "[" << c++ << "]alpha_x()";        ofs << m_sep;
+            ofs << "[" << c++ << "]alpha_y()";        ofs << m_sep;
             ofs << "[" << c++ << "]beta_x(m)";        ofs << m_sep;
             ofs << "[" << c++ << "]beta_y(m)";        ofs << m_sep;
             ofs << "[" << c++ << "]charge(C)";        ofs << std::endl;
@@ -124,6 +128,7 @@ BeamRelevant::BeamRelevant (std::string rd_name)
             ofs << "[" << c++ << "]gamma_rms()";      ofs << m_sep;
             ofs << "[" << c++ << "]emittance_x(m)";   ofs << m_sep;
             ofs << "[" << c++ << "]emittance_z(m)";   ofs << m_sep;
+            ofs << "[" << c++ << "]alpha_x()";        ofs << m_sep;
             ofs << "[" << c++ << "]beta_x(m)";        ofs << m_sep;
             ofs << "[" << c++ << "]charge(C)";        ofs << std::endl;
 #elif (defined WARPX_DIM_1D_Z)
@@ -373,9 +378,11 @@ void BeamRelevant::ComputeDiags (int step)
         m_data[14] = std::sqrt(x_ms*ux_ms-xux*xux) / PhysConst::c;
         m_data[15] = std::sqrt(y_ms*uy_ms-yuy*yuy) / PhysConst::c;
         m_data[16] = std::sqrt(z_ms*uz_ms-zuz*zuz) / PhysConst::c;
-        m_data[17] = (PhysConst::c * x_ms) / std::sqrt(x_ms*ux_ms-xux*xux);
-        m_data[18] = (PhysConst::c * y_ms) / std::sqrt(y_ms*uy_ms-yuy*yuy);
-        m_data[19] = charge;
+        m_data[17] = - (PhysConst::c * xux) / std::sqrt(x_ms*ux_ms-xux*xux);
+        m_data[18] = - (PhysConst::c * yuy) / std::sqrt(y_ms*uy_ms-yuy*yuy);
+        m_data[19] = (PhysConst::c * x_ms) / std::sqrt(x_ms*ux_ms-xux*xux);
+        m_data[20] = (PhysConst::c * y_ms) / std::sqrt(y_ms*uy_ms-yuy*yuy);
+        m_data[21] = charge;
 #elif (defined WARPX_DIM_XZ)
         m_data[0]  = x_mean;
         m_data[1]  = z_mean;
@@ -391,8 +398,9 @@ void BeamRelevant::ComputeDiags (int step)
         m_data[11] = std::sqrt(gm_ms);
         m_data[12] = std::sqrt(x_ms*ux_ms-xux*xux) / PhysConst::c;
         m_data[13] = std::sqrt(z_ms*uz_ms-zuz*zuz) / PhysConst::c;
-        m_data[14] = (PhysConst::c * x_ms) / std::sqrt(x_ms*ux_ms-xux*xux);
-        m_data[15] = charge;
+        m_data[14] = - (PhysConst::c * xux) / std::sqrt(x_ms*ux_ms-xux*xux);
+        m_data[15] = (PhysConst::c * x_ms) / std::sqrt(x_ms*ux_ms-xux*xux);
+        m_data[16] = charge;
         amrex::ignore_unused(y_mean, y_ms, yuy);
 #elif (defined WARPX_DIM_1D_Z)
         m_data[0]  = z_mean;
@@ -403,8 +411,8 @@ void BeamRelevant::ComputeDiags (int step)
         m_data[5]  = std::sqrt(z_ms);
         m_data[6]  = std::sqrt(ux_ms) * m;
         m_data[7]  = std::sqrt(uy_ms) * m;
-        m_data[8] = std::sqrt(uz_ms) * m;
-        m_data[9] = std::sqrt(gm_ms);
+        m_data[8]  = std::sqrt(uz_ms) * m;
+        m_data[9]  = std::sqrt(gm_ms);
         m_data[10] = std::sqrt(z_ms*uz_ms-zuz*zuz) / PhysConst::c;
         m_data[11] = charge;
         amrex::ignore_unused(x_mean, x_ms, xux, y_mean, y_ms, yuy);