create HYSPEC intensity matrix (#24)

* inital tests setup

* more tests

* change return type to dict

* test calculate different graph options

* test deltaE

* logic and test to handle DeltaE < -Ei

* fixing maths in model.

* chagne plot results to universal key workd "intensity" and update tests

* minor code clean ups

* test -S2

* Clean up some code

* Fix accidental renaming

* Clean up code and add comments


* Clean up code and fix pre-commit

* Clean up test

* Add consistency check
---------

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

scripts/PPTPmodel-Kyle.py

@@ -57,8 +57,8 @@ def polarization_powder(self,Ei, EMin, S2, alpha_p, plot_options="alpha",left=Tr
         Qx= (-1 if left else 1) * kf2d*np.sqrt((1-cos_theta**2))
 
         cos_ang_PQ = (Qx*Px + Qz*Pz)/Q2d/np.sqrt(Px**2+Pz**2)
-        cos_ang_PQ[cos_ang_PQ**2 <0.4] = np.nan 
-        cos_ang_PQ[cos_ang_PQ**2 >0.6] = np.nan
+        # cos_ang_PQ[cos_ang_PQ**2 <0.4] = np.nan 
+        # cos_ang_PQ[cos_ang_PQ**2 >0.6] = np.nan
         ang_PQ = np.degrees(np.arccos(cos_ang_PQ))
 
         kf=np.sqrt(Ei-E)*SE2K
@@ -124,7 +124,7 @@ def polarization_powder(self,Ei, EMin, S2, alpha_p, plot_options="alpha",left=Tr
 
 if __name__ == "__main__":
     obj = Model()
-    output = obj.polarization_powder(20, None, 60, 30,plot_options="alpha",left=True)
+    output = obj.polarization_powder(20, None, 40, 0,plot_options="cos^2(a)",left=True)
     Q_low, Q_hi, E, Q2d, E2d, ang_PQ = output[0], output[1], output[2], output[3], output[4], output[5]
 
     fig, ax = plt.subplots()

src/hyspecppt/hppt/experiment_settings.py

@@ -21,6 +21,11 @@
 DEFAULT_MODE = dict(current_experiment_type="single_crystal")
 # maximum momentum transfer
 MAX_MODQ = 15
+# number of points in the plot
+N_POINTS = 200
+# tank half-width
+TANK_HALF_WIDTH = 30.0
+
 
 # invalid style
 INVALID_QLINEEDIT = """

src/hyspecppt/hppt/hppt_model.py

@@ -11,7 +11,9 @@
     DEFAULT_LATTICE,
     DEFAULT_MODE,
     MAX_MODQ,
+    N_POINTS,
     PLOT_TYPES,
+    TANK_HALF_WIDTH,
 )
 
 logger = logging.getLogger("hyspecppt")
@@ -20,19 +22,19 @@
 class SingleCrystalParameters:
     """Model for single crystal calculations"""
 
-    a: float = DEFAULT_LATTICE["a"]
-    b: float = DEFAULT_LATTICE["b"]
-    c: float = DEFAULT_LATTICE["c"]
-    alpha: float = DEFAULT_LATTICE["alpha"]
-    beta: float = DEFAULT_LATTICE["beta"]
-    gamma: float = DEFAULT_LATTICE["gamma"]
-    h: float = DEFAULT_LATTICE["h"]
-    k: float = DEFAULT_LATTICE["k"]
-    l: float = DEFAULT_LATTICE["l"]
+    a: float
+    b: float
+    c: float
+    alpha: float
+    beta: float
+    gamma: float
+    h: float
+    k: float
+    l: float
 
     def __init__(self) -> None:
         """Constructor"""
-        return
+        self.set_parameters(DEFAULT_LATTICE)
 
     def set_parameters(self, params: dict[str, float]) -> None:
         """Store single crystal parameters
@@ -54,63 +56,58 @@ def set_parameters(self, params: dict[str, float]) -> None:
 
     def get_parameters(self) -> dict[str, float]:
         """Returns all the parameters as a dictionary"""
-        try:
-            return dict(
-                a=self.a,
-                b=self.b,
-                c=self.c,
-                alpha=self.alpha,
-                beta=self.beta,
-                gamma=self.gamma,
-                h=self.h,
-                k=self.k,
-                l=self.l,
-            )
-        except AttributeError:
-            logger.error("The parameters were not initialized-get_parameters-SingleCrystalParameters")
+        return dict(
+            a=self.a,
+            b=self.b,
+            c=self.c,
+            alpha=self.alpha,
+            beta=self.beta,
+            gamma=self.gamma,
+            h=self.h,
+            k=self.k,
+            l=self.l,
+        )
 
     def calculate_modQ(self) -> float:
         """Returns |Q| from lattice parameters and h, k, l"""
-        try:
-            ca = np.cos(np.radians(self.alpha))
-            sa = np.sin(np.radians(self.alpha))
-            cb = np.cos(np.radians(self.beta))
-            sb = np.sin(np.radians(self.beta))
-            cg = np.cos(np.radians(self.gamma))
-            sg = np.sin(np.radians(self.gamma))
-            vabg = np.sqrt(1 - ca**2 - cb**2 - cg**2 + 2 * ca * cb * cg)
-            astar = sa / (self.a * vabg)
-            bstar = sb / (self.b * vabg)
-            cstar = sg / (self.c * vabg)
-            cas = (cb * cg - ca) / (sb * sg)  # noqa: F841
-            cbs = (cg * ca - cb) / (sg * sa)
-            cgs = (ca * cb - cg) / (sa * sb)
-
-            # B matrix
-            B = np.array(
-                [
-                    [astar, bstar * cgs, cstar * cbs],
-                    [0, bstar * np.sqrt(1 - cgs**2), -cstar * np.sqrt(1 - cbs**2) * ca],
-                    [0, 0, 1.0 / self.c],
-                ]
-            )
-
-            modQ = 2 * np.pi * np.linalg.norm(B.dot([self.h, self.k, self.l]))
-            return modQ
-        except AttributeError:
-            logger.error("The parameters were not initialized-calculate_modQ")
+        ca = np.cos(np.radians(self.alpha))
+        sa = np.sin(np.radians(self.alpha))
+        cb = np.cos(np.radians(self.beta))
+        sb = np.sin(np.radians(self.beta))
+        cg = np.cos(np.radians(self.gamma))
+        sg = np.sin(np.radians(self.gamma))
+        vabg = np.sqrt(1 - ca**2 - cb**2 - cg**2 + 2 * ca * cb * cg)
+        astar = sa / (self.a * vabg)
+        bstar = sb / (self.b * vabg)
+        cstar = sg / (self.c * vabg)
+        cas = (cb * cg - ca) / (sb * sg)  # noqa: F841
+        cbs = (cg * ca - cb) / (sg * sa)
+        cgs = (ca * cb - cg) / (sa * sb)
+
+        # B matrix
+        B = np.array(
+            [
+                [astar, bstar * cgs, cstar * cbs],
+                [0, bstar * np.sqrt(1 - cgs**2), -cstar * np.sqrt(1 - cbs**2) * ca],
+                [0, 0, 1.0 / self.c],
+            ]
+        )
+
+        modQ = 2 * np.pi * np.linalg.norm(B.dot([self.h, self.k, self.l]))
+        return modQ
 
 
 class CrosshairParameters:
     """Model for the crosshair parameters"""
 
-    modQ: float = DEFAULT_CROSSHAIR["modQ"]
-    DeltaE: float = DEFAULT_CROSSHAIR["DeltaE"]
-    current_experiment_type: str = DEFAULT_MODE["current_experiment_type"]
+    modQ: float
+    DeltaE: float
+    current_experiment_type: str
     sc_parameters: SingleCrystalParameters
-    experiment_types = ["single_crystal", "powder"]
 
-    def __init__(self):
+    def __init__(self) -> None:
+        """Constructor"""
+        self.set_crosshair(**DEFAULT_MODE, **DEFAULT_CROSSHAIR)
         self.sc_parameters = SingleCrystalParameters()
 
     def set_crosshair(self, current_experiment_type: str, DeltaE: float = None, modQ: float = None) -> None:
@@ -132,18 +129,23 @@ def get_crosshair(self) -> dict[str, float]:
             return dict(DeltaE=self.DeltaE, modQ=modQ)
         return dict(DeltaE=self.DeltaE, modQ=self.modQ)
 
+    def get_experiment_type(self) -> str:
+        """Return experiment type"""
+        return self.current_experiment_type
+
 
 class HyspecPPTModel:
     """Main model"""
 
-    Ei: float = DEFAULT_EXPERIMENT["Ei"]
-    S2: float = DEFAULT_EXPERIMENT["S2"]
-    alpha_p: float = DEFAULT_EXPERIMENT["alpha_p"]
-    plot_type: str = DEFAULT_EXPERIMENT["plot_type"]
+    Ei: float
+    S2: float
+    alpha_p: float
+    plot_type: str
     cp: CrosshairParameters
 
     def __init__(self):
         """Constructor"""
+        self.set_experiment_data(**DEFAULT_EXPERIMENT)
         self.cp = CrosshairParameters()
 
     def set_single_crystal_data(self, params: dict[str, float]) -> None:
@@ -165,70 +167,69 @@ def set_experiment_data(self, Ei: float, S2: float, alpha_p: float, plot_type: s
         self.plot_type = plot_type
 
     def get_experiment_data(self) -> dict[str, float]:
-        data = dict()
-
-        data["Ei"] = self.Ei
-        data["S2"] = self.S2
-        data["alpha_p"] = self.alpha_p
-        data["plot_type"] = self.plot_type
+        data = dict(Ei=self.Ei, S2=self.S2, alpha_p=self.alpha_p, plot_type=self.plot_type)
         return data
 
-    def get_graph_data(self) -> list[float, float, float, list, list, list]:
-        return self.calculate_graph_data()
-
-    def calculate_graph_data(self) -> list[float]:
-        """Returns a list of [Q_low, Q_hi, E, Q2d, E2d, data of plot_types]"""
-        try:
-            SE2K = np.sqrt(2e-3 * e * m_n) * 1e-10 / hbar
-            # def Ei, Emin = - Ei to create Qmin, Qmax to generate plot range
-            if self.cp.DeltaE is not None and self.cp.DeltaE < -self.Ei:
-                EMin = -self.Ei
-            E = np.linspace(EMin, self.Ei * 0.9, 200)
-
-            kfmin = np.sqrt(self.Ei - EMin) * SE2K
-            ki = np.sqrt(self.Ei) * SE2K
-
-            # S2= 60
-            # Create Qmin and Qmax
-            Qmax = np.sqrt(ki**2 + kfmin**2 - 2 * ki * kfmin * np.cos(np.radians(self.S2 + 30)))  # Q=ki or Q=
-            Qmin = 0
-            Q = np.linspace(Qmin, Qmax, 200)
-
-            # Create 2D array
-            E2d, Q2d = np.meshgrid(E, Q)
-
-            Ef2d = self.Ei - E2d
-            kf2d = np.sqrt(Ef2d) * SE2K
-
-            Px = np.cos(np.radians(self.alpha_p))
-            Pz = np.sin(np.radians(self.alpha_p))
-
-            cos_theta = (ki**2 + kf2d**2 - Q2d**2) / (2 * ki * kf2d)
-            cos_theta[cos_theta < np.cos(np.radians(self.S2 + 30))] = np.nan
-            cos_theta[cos_theta > np.cos(np.radians(self.S2 - 30))] = np.nan  # cos(30)
-
-            Qz = ki - kf2d * cos_theta
-
-            if self.S2 >= 30:
-                Qx = (-1) * kf2d * np.sqrt((1 - cos_theta**2))
-            elif self.S2 <= -30:
-                Qx = kf2d * np.sqrt((1 - cos_theta**2))
-
-            cos_ang_PQ = (Qx * Px + Qz * Pz) / Q2d / np.sqrt(Px**2 + Pz**2)
-            ang_PQ = np.degrees(np.arccos(cos_ang_PQ))
-
-            kf = np.sqrt(self.Ei - E) * SE2K
-
-            Q_low = np.sqrt(ki**2 + kf**2 - 2 * ki * kf * np.cos(np.radians(self.S2 - 30)))
-            Q_hi = np.sqrt(ki**2 + kf**2 - 2 * ki * kf * np.cos(np.radians(self.S2 + 30)))
-
-            if self.plot_type == PLOT_TYPES[0]:  # alpha
-                return [Q_low, Q_hi, E, Q2d, E2d, ang_PQ]
-
-            if self.plot_type == PLOT_TYPES[1]:  # cos^2(alpha)
-                return [Q_low, Q_hi, E, Q2d, E2d, np.cos(np.radians(ang_PQ)) ** 2]
-
-            if self.plot_type == PLOT_TYPES[2]:  # "(cos^2(a)+1)/2"
-                return [Q_low, Q_hi, E, Q2d, E2d, (np.cos(np.radians(ang_PQ)) ** 2 + 1) / 2]
-        except AttributeError:
-            logger.error("The parameters were not initialized")
+    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]"""
+        # constant to transform from energy in meV to momentum in Angstrom^-1
+        SE2K = np.sqrt(2e-3 * e * m_n) * 1e-10 / hbar
+
+        # adjust minimum energy
+        if self.cp.DeltaE is not None and self.cp.DeltaE <= -self.Ei:
+            EMin = 1.2 * self.cp.DeltaE
+        else:
+            EMin = -self.Ei
+
+        E = np.linspace(EMin, self.Ei * 0.9, N_POINTS)
+
+        # Calculate lines for the edges of the tank
+        ki = np.sqrt(self.Ei) * SE2K
+        kf = np.sqrt(self.Ei - E) * SE2K
+        Q_low = np.sqrt(ki**2 + kf**2 - 2 * ki * kf * np.cos(np.radians(np.abs(self.S2) - TANK_HALF_WIDTH)))
+        Q_hi = np.sqrt(ki**2 + kf**2 - 2 * ki * kf * np.cos(np.radians(np.abs(self.S2) + TANK_HALF_WIDTH)))
+
+        # Create 2D array
+        Q = np.linspace(0, np.max(Q_hi), N_POINTS)
+        E2d, Q2d = np.meshgrid(E, Q)
+        Ef2d = self.Ei - E2d
+        kf2d = np.sqrt(Ef2d) * SE2K
+
+        # Calculate the angle between Q and z axis. Set to NAN values outside the detector range
+        cos_theta = (ki**2 + kf2d**2 - Q2d**2) / (2 * ki * kf2d)
+        cos_theta[cos_theta < np.cos(np.radians(np.abs(self.S2) + TANK_HALF_WIDTH))] = np.nan
+        cos_theta[cos_theta > np.cos(np.radians(np.abs(self.S2) - TANK_HALF_WIDTH))] = np.nan
+
+        # Calculate Qz
+        Qz = ki - kf2d * cos_theta
+
+        # Calculate Qx. Note that is in the opposite direction as detector position
+        if self.S2 >= TANK_HALF_WIDTH:
+            Qx = (-1) * kf2d * np.sqrt((1 - cos_theta**2))
+        elif self.S2 <= -TANK_HALF_WIDTH:
+            Qx = kf2d * np.sqrt((1 - cos_theta**2))
+
+        # Transform polarization angle in the lab frame to vector
+        Px = np.cos(np.radians(self.alpha_p))
+        Pz = np.sin(np.radians(self.alpha_p))
+
+        # Calculate angle between polarization vector and momentum transfer
+        cos_ang_PQ = (Qx * Px + Qz * Pz) / Q2d / np.sqrt(Px**2 + Pz**2)
+
+        # Select return value for intensity
+        if self.plot_type == PLOT_TYPES[0]:  # alpha
+            ang_PQ = np.arccos(cos_ang_PQ)
+            intensity = np.degrees(ang_PQ)
+        elif self.plot_type == PLOT_TYPES[1]:  # cos^2(alpha)
+            intensity = cos_ang_PQ**2
+        elif self.plot_type == PLOT_TYPES[2]:  # "(cos^2(a)+1)/2"
+            intensity = (cos_ang_PQ**2 + 1) / 2
+
+        return dict(
+            Q_low=Q_low,
+            Q_hi=Q_hi,
+            E=E,
+            Q2d=Q2d,
+            E2d=E2d,
+            intensity=intensity,
+        )

src/hyspecppt/hppt/hppt_presenter.py

@@ -1,6 +1,6 @@
 """Presenter for the Main tab"""
 
-from .experiment_settings import DEFAULT_CROSSHAIR, DEFAULT_EXPERIMENT, DEFAULT_LATTICE, DEFAULT_MODE, PLOT_TYPES
+from .experiment_settings import PLOT_TYPES
 
 
 class HyspecPPTPresenter:
@@ -20,20 +20,13 @@ def __init__(self, view: any, model: any):
         self.view.connect_sc_mode_switch(self.handle_switch_to_sc)
 
         # populate fields
-        self.view.sc_widget.set_values(DEFAULT_LATTICE)
+        self.view.sc_widget.set_values(self.model.get_single_crystal_data())
         self.view.experiment_widget.initializeCombo(PLOT_TYPES)
-        self.view.experiment_widget.set_values(DEFAULT_EXPERIMENT)
-        self.view.crosshair_widget.set_values(DEFAULT_CROSSHAIR)
-
-        # model init
-        # to be removed needs to happen in the model
-        self.model.set_experiment_data(**DEFAULT_EXPERIMENT)
-        self.model.set_crosshair_data(**DEFAULT_CROSSHAIR, **DEFAULT_MODE)
-        self.model.set_single_crystal_data(params=DEFAULT_LATTICE)
+        self.view.experiment_widget.set_values(self.model.get_experiment_data())
 
         # set default selection mode
         experiment_type = self.view.selection_widget.powder_label
-        if DEFAULT_MODE["current_experiment_type"].startswith("single"):
+        if self.model.cp.get_experiment_type().startswith("single"):
             experiment_type = self.view.selection_widget.sc_label
         self.view.selection_widget.selector_init(experiment_type)  # pass the default mode from experiment type