Calculate Q-beam (Q-ki) (#44)

* calculate angle between Q and beam(z)

* calculate kiQ

* calculate Q-beam angle by hand

* test use scenarios

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* test use scenarios

* resolve merge conflicts

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

src/hyspecppt/hppt/hppt_model.py

@@ -247,6 +247,18 @@ def check_plot_update(self, deltaE) -> bool:
         # check if it is the same
         return self.Emin != Emin
 
+    def get_ang_Q_beam(self) -> float:
+        """Returns the angle between Q and the beam"""
+        SE2K = np.sqrt(2e-3 * e * m_n) * 1e-10 / hbar
+        crosshair_data = self.get_crosshair_data()
+        deltaE = crosshair_data["DeltaE"]
+        modQ = crosshair_data["modQ"]
+        ki = np.sqrt(self.Ei) * SE2K
+        kf = np.sqrt(self.Ei - deltaE) * SE2K
+        with np.errstate(all="ignore"):  # ignore the state when momentum energy not conserved
+            cos_kiQ = (ki**2 + modQ**2 - kf**2) / (2 * ki * modQ)
+            return np.degrees(np.arccos(cos_kiQ)) if self.S2 < 0 else -np.degrees(np.arccos(cos_kiQ))
+
     def calculate_graph_data(self) -> dict[str, np.array]:
         """Returns a dictionary of arrays [Q_low, Q_hi, E, Q2d, E2d, data of plot_types]
 

tests/hppt_model/test_hyspecpptmodel.py

@@ -221,6 +221,67 @@ def test_zero_alpha():
     assert np.isclose(model.calculate_graph_data()["intensity"][94][105], 0.209092)
 
 
+def test_Q_beam_angle():
+    """Test the angle between Q and beam based on hyspecppt/scripts/test_Q-beam_angle.pdf"""
+    model = HyspecPPTModel()
+    # Only Ei matters as it is used in calculating ki
+    assert model.Ei == 20.0
+
+    assert np.isnan(model.get_ang_Q_beam())
+    model.set_crosshair_data("powder", DeltaE=0, modQ=3.107)
+    assert np.isclose(model.get_ang_Q_beam(), -60, 0.01)
+
+    model.set_crosshair_data("powder", DeltaE=0, modQ=2.0)
+    assert np.isclose(model.get_ang_Q_beam(), -71.22, 0.01)
+
+
+def test_Q_beam_angle_opposite_to_S2():
+    """Test that changing S2 flip signs of Q-beam angle"""
+    model = HyspecPPTModel()
+    model.set_experiment_data(Ei=20.0, S2=-30, alpha_p=0, plot_type=PLOT_TYPES[0])
+    model.set_crosshair_data("powder", DeltaE=0, modQ=3.107)
+    assert np.isclose(model.get_ang_Q_beam(), 60, 0.01)
+    model.set_crosshair_data("powder", DeltaE=0, modQ=2.0)
+    assert np.isclose(model.get_ang_Q_beam(), 71.22, 0.01)
+
+
+def test_scattering_triagnle_not_closed():
+    """Test intermediate state when scattering triangle not closed"""
+    model = HyspecPPTModel()
+    model.set_crosshair_data("powder", DeltaE=5, modQ=0)
+    assert np.isnan(model.get_ang_Q_beam())
+
+
+def test_maxq():
+    """Test intermediate state when scattering triangle not closed"""
+    model = HyspecPPTModel()
+    model.set_crosshair_data("powder", DeltaE=0, modQ=6.22)
+    assert np.isnan(model.get_ang_Q_beam())
+
+    model.set_crosshair_data("powder", DeltaE=0, modQ=6.21)
+    assert np.isclose(model.get_ang_Q_beam(), -1.926, 0.001)
+
+
+def test_inelastic_position():
+    """Test intermediate state when scattering triangle not closed"""
+    model = HyspecPPTModel()
+    model.set_crosshair_data("powder", DeltaE=10, modQ=3)
+    assert np.isclose(model.get_ang_Q_beam(), -42.122, 0.001)
+
+    model.set_crosshair_data("powder", DeltaE=10, modQ=2)
+    assert np.isclose(model.get_ang_Q_beam(), -44.747, 0.001)
+
+
+def test_single_crystal_q_beam_mode():
+    """Test changing single crystal data will also calculate correct Q-beam angle"""
+    model = HyspecPPTModel()
+    model.set_crosshair_data("single_crystal", DeltaE=10)
+    params = dict(a=5.0, b=6.0, c=7.0, alpha=90.0, beta=90.0, gamma=120.0, h=1.0, k=2.0, l=3.0)
+    model.set_single_crystal_data(params)
+    assert np.isclose(model.get_crosshair_data()["modQ"], 4.326)
+    assert np.isclose(model.get_ang_Q_beam(), -28.864, 0.001)
+
+
 def test_cos2_sin2():
     """Test calculating different graph data"""
     model = HyspecPPTModel()