diff --git a/src/tavi/instrument/components/collimators.py b/src/tavi/instrument/components/collimators.py new file mode 100644 index 00000000..e097e737 --- /dev/null +++ b/src/tavi/instrument/components/collimators.py @@ -0,0 +1,67 @@ +# -*- coding: utf-8 -*- +from typing import Optional + +import numpy as np + +from tavi.instrument.components.component_base import TASComponent + + +class Collimators(TASComponent): + """collimitor divergences, in mins of arc""" + + def __init__( + self, + param_dict: Optional[dict] = None, + component_name: str = "collimitors", + ): + # defalut values + self.h_pre_mono: float = 30.0 # mins of arc + self.h_pre_sample: float = 30.0 + self.h_post_sample: float = 30.0 + self.h_post_ana: float = 30.0 + self.v_pre_mono: float = 30.0 + self.v_pre_sample: float = 30.0 + self.v_post_sample: float = 30.0 + self.v_post_ana: float = 30.0 + + super().__init__(param_dict, component_name) + + @property + def horizontal_divergence(self) -> list: + """list of horizontal divergence in minitus of arc""" + return [ + self.h_pre_mono, + self.h_pre_sample, + self.h_post_sample, + self.h_post_ana, + ] + + @property + def _horizontal_divergence(self) -> list: + """list of horizontal divergence in radian""" + return [ + np.deg2rad(self.h_pre_mono / 60), + np.deg2rad(self.h_pre_sample / 60), + np.deg2rad(self.h_post_sample / 60), + np.deg2rad(self.h_post_ana / 60), + ] + + @property + def vertical_divergence(self) -> list: + """list of vertical divergence in minutes of arcs""" + return [ + self.v_pre_mono, + self.v_pre_sample, + self.v_post_sample, + self.v_post_ana, + ] + + @property + def _vertical_divergence(self) -> list: + """list of vertical divergence in radian""" + return [ + np.deg2rad(self.v_pre_mono / 60), + np.deg2rad(self.v_pre_sample / 60), + np.deg2rad(self.v_post_sample / 60), + np.deg2rad(self.v_post_ana / 60), + ] diff --git a/src/tavi/instrument/components/component_base.py b/src/tavi/instrument/components/component_base.py new file mode 100644 index 00000000..a6b58e44 --- /dev/null +++ b/src/tavi/instrument/components/component_base.py @@ -0,0 +1,127 @@ +# -*- coding: utf-8 -*- +from typing import Optional, Union + +import numpy as np + +from tavi.utilities import cm2angstrom + + +class TASComponent(object): + """ + A component of the triple-axis spectrometer + + Attributes: + type (str): identify which component + Methods: + update(param, value): update a parameter + """ + + def __init__( + self, + param_dict: Optional[dict] = None, + component_name: str = "", + ): + """Load parameters if provided.""" + self.component_name = component_name + if param_dict is not None: + for key, val in param_dict.items(): + setattr(self, key, val) + + def update( + self, + param: str, + value: Union[str, float, int], + ) -> None: + """update a paramter if exist""" + if hasattr(self, param): + setattr(self, param, value) + print(f"Setting {self.component_name} parameter {param} to {value}") + else: + print(f"{self.component_name} has no attribute {param}") + + @staticmethod + def _min2rad( + angle_in_minutes: Optional[float], + param_name: str = "Angle", + ) -> Optional[float]: + """Convert from minutes to radian is angle is not None""" + if angle_in_minutes is None: + print(f"{param_name} is None.") + return None + return np.deg2rad(angle_in_minutes / 60) + + @staticmethod + def _cm2angstrom_given_shape( + length_in_cm: Optional[float], + shape: str, + param_name: str = "Length", + ) -> Optional[float]: + """ + Convert length from centimeter to angstrom. + + Apply correction based on shape, for resolution calculation. + Divide length by np.sqrt(12) if rectangular, + Divive diameter D by 4 if circular or spherical. + """ + if length_in_cm is None: + print(f"{param_name} is None.") + return None + match shape: + case "rectangular": + return length_in_cm / np.sqrt(12) * cm2angstrom + case "circular" | "spherical": + return length_in_cm / 4 * cm2angstrom + case _: + print("Unrecognized shape. Needs to be rectangular or circular/spherical.") + return None + + @staticmethod + def _cm2angstrom( + length_in_cm: Optional[float], + param_name: str = "Length", + ) -> Optional[float]: + """ + Convert length from centimeter to angstrom. + """ + if length_in_cm is None: + print(f"{param_name} is None.") + return None + + return length_in_cm * cm2angstrom + + +class Distances(TASComponent): + """Distance between components, in units of centimeters""" + + def __init__( + self, + param_dict: Optional[dict] = None, + component_name: str = "distances", + ): + self.src_mono: Optional[float] = None + self.mono_sample: Optional[float] = None + self.sample_ana: Optional[float] = None + self.ana_det: Optional[float] = None + self.mono_monitor: Optional[float] = None + + super().__init__(param_dict, component_name) + + @property + def _src_mono(self): + return TASComponent._cm2angstrom(self.src_mono, "Source to monochromator distance") + + @property + def _mono_sample(self): + return TASComponent._cm2angstrom(self.src_mono, "Monochromator to sample distance") + + @property + def _sample_ana(self): + return TASComponent._cm2angstrom(self.src_mono, "Sample to analyzer distance") + + @property + def _ana_det(self): + return TASComponent._cm2angstrom(self.src_mono, "Analyzer to detector distance") + + @property + def _mono_monitor(self): + return TASComponent._cm2angstrom(self.src_mono, "Monochromator to monitor distance") diff --git a/src/tavi/instrument/components/detector.py b/src/tavi/instrument/components/detector.py new file mode 100644 index 00000000..12bec74a --- /dev/null +++ b/src/tavi/instrument/components/detector.py @@ -0,0 +1,35 @@ +# -*- coding: utf-8 -*- +from typing import Literal, Optional + +from tavi.instrument.components.component_base import TASComponent + + +# TODO not implemented in resolution calculation +class Detector(TASComponent): + """ + Neutron detector + + Note: + Need to think about detector efficiency, saturation + """ + + def __init__( + self, + param_dict: Optional[dict] = None, + component_name: str = "detector", + ): + self.shape: Literal["rectangular", "spherical"] = "rectangular" + self.width: Optional[float] = None + self.height: Optional[float] = None + + super().__init__(param_dict, component_name) + + @property + def _width(self): + """width in angstrom, with correction based on shape, for resolution calculation""" + return TASComponent._cm2angstrom_given_shape(self.width, self.shape, "Detector width") + + @property + def _height(self): + """Height in angstrom, with correction based on shape, for resolution calculation""" + return TASComponent._cm2angstrom_given_shape(self.height, self.shape, "Detector height") diff --git a/src/tavi/instrument/goni.py b/src/tavi/instrument/components/goni.py similarity index 98% rename from src/tavi/instrument/goni.py rename to src/tavi/instrument/components/goni.py index 892a227e..c41c7adf 100644 --- a/src/tavi/instrument/goni.py +++ b/src/tavi/instrument/components/goni.py @@ -2,7 +2,7 @@ import numpy as np -from tavi.instrument.tas_cmponents import TASComponent +from tavi.instrument.components.component_base import TASComponent from tavi.utilities import MotorAngles diff --git a/src/tavi/instrument/components/guide.py b/src/tavi/instrument/components/guide.py new file mode 100644 index 00000000..537bcd8f --- /dev/null +++ b/src/tavi/instrument/components/guide.py @@ -0,0 +1,29 @@ +# -*- coding: utf-8 -*- +from typing import Optional + +from tavi.instrument.components.component_base import TASComponent + + +class Guide(TASComponent): + """Neutron guide in the monochromator drum""" + + def __init__( + self, + param_dict: Optional[dict] = None, + component_name: str = "guide", + ): + self.in_use: bool = False + self.div_h: Optional[float] = None + self.div_v: Optional[float] = None + + super().__init__(param_dict, component_name) + + @property + def _div_h(self): + """Horizontal divergence in radian""" + return TASComponent._min2rad(self.div_h, "Guide horizontal divergence") + + @property + def _div_v(self): + """Vertical divergence in radian""" + return TASComponent._min2rad(self.div_v, "Guide vertical divergence") diff --git a/src/tavi/instrument/components/monitor.py b/src/tavi/instrument/components/monitor.py new file mode 100644 index 00000000..efcfa2fb --- /dev/null +++ b/src/tavi/instrument/components/monitor.py @@ -0,0 +1,35 @@ +# -*- coding: utf-8 -*- +from typing import Literal, Optional + +from tavi.instrument.components.component_base import TASComponent + + +# TODO not implemented in resolution calculation +class Monitor(TASComponent): + """ + Neutron Monitor + + Note: + Needs to think about monitor efficiency + """ + + def __init__( + self, + param_dict: Optional[dict] = None, + component_name: str = "monitor", + ): + self.shape: Literal["rectangular", "spherical"] = "rectangular" + self.width: Optional[float] = None + self.height: Optional[float] = None + + super().__init__(param_dict, component_name) + + @property + def _width(self): + """width in angstrom, with correction based on shape, for resolution calculation""" + return TASComponent._cm2angstrom_given_shape(self.width, self.shape, "Monitor width") + + @property + def _height(self): + """Height in angstrom, with correction based on shape, for resolution calculation""" + return TASComponent._cm2angstrom_given_shape(self.height, self.shape, "Monitor height") diff --git a/src/tavi/instrument/mono_ana.py b/src/tavi/instrument/components/mono_ana.py similarity index 97% rename from src/tavi/instrument/mono_ana.py rename to src/tavi/instrument/components/mono_ana.py index 46a355c8..72d00a46 100644 --- a/src/tavi/instrument/mono_ana.py +++ b/src/tavi/instrument/components/mono_ana.py @@ -2,7 +2,7 @@ import numpy as np -from tavi.instrument.tas_cmponents import TASComponent +from tavi.instrument.components.component_base import TASComponent # --------------------------------------------------------------- # d_spacing table from Shirane Appendix 3, in units of Angstrom diff --git a/src/tavi/instrument/components/source.py b/src/tavi/instrument/components/source.py new file mode 100644 index 00000000..dbcd3f38 --- /dev/null +++ b/src/tavi/instrument/components/source.py @@ -0,0 +1,30 @@ +# -*- coding: utf-8 -*- +from typing import Literal, Optional + +from tavi.instrument.components.component_base import TASComponent + + +class Source(TASComponent): + """Neutron source""" + + def __init__( + self, + param_dict: Optional[dict] = None, + component_name="source", + ): + self.name: str = "HFIR" + self.shape: Literal["rectangular", "circular"] = "rectangular" + self.width: Optional[float] = None + self.height: Optional[float] = None + + super().__init__(param_dict, component_name) + + @property + def _width(self): + """width in angstrom, with correction based on shape, for resolution calculation""" + return TASComponent._cm2angstrom_given_shape(self.width, self.shape, "Source width") + + @property + def _height(self): + """Height in angstrom, with correction based on shape, for resolution calculation""" + return TASComponent._cm2angstrom_given_shape(self.height, self.shape, "Source height") diff --git a/src/tavi/instrument/resolution/cooper_nathans.py b/src/tavi/instrument/resolution/cooper_nathans.py index 10b9e291..f380922a 100755 --- a/src/tavi/instrument/resolution/cooper_nathans.py +++ b/src/tavi/instrument/resolution/cooper_nathans.py @@ -1,43 +1,40 @@ -from typing import Optional +from typing import Optional, Union import numpy as np -from tavi.instrument.mono_ana import MonoAna +from tavi.instrument.components.collimators import Collimators +from tavi.instrument.components.mono_ana import MonoAna from tavi.instrument.resolution.ellipsoid import ResoEllipsoid from tavi.instrument.tas import TAS -from tavi.instrument.tas_cmponents import Collimators -from tavi.utilities import get_angle_bragg, get_angle_from_triangle, ksq2eng, rotation_matrix_2d, sig2fwhm +from tavi.utilities import ( + get_angle_bragg, + get_angle_from_triangle, + ksq2eng, + rotation_matrix_2d, + sig2fwhm, + spice_to_mantid, +) class CN(TAS): """Copper-Nathans method - Attibutes: - _mat_f - _mat_g - Methods: + validate_instrument_parameters cooper_nathans """ # 4 soller slits collimators NUM_COLLS = 4 - IDX_COLL0_H = 0 - IDX_COLL0_V = 2 - IDX_COLL1_H = 1 - IDX_COLL1_V = 3 - IDX_COLL2_H = 4 - IDX_COLL2_V = 6 - IDX_COLL3_H = 5 - IDX_COLL3_V = 7 + IDX_COLL0_H, IDX_COLL0_V = 0, 2 + IDX_COLL1_H, IDX_COLL1_V = 1, 3 + IDX_COLL2_H, IDX_COLL2_V = 4, 6 + IDX_COLL3_H, IDX_COLL3_V = 5, 7 # 1 monochromator and 1 analyzer - NUM_MONOS = 1 - NUM_ANAS = 1 - IDX_MONO0_H = 0 - IDX_MONO0_V = 1 - IDX_ANA0_H = 2 - IDX_ANA0_V = 3 + NUM_MONOS, NUM_ANAS = 1, 1 + IDX_MONO0_H, IDX_MONO0_V = 0, 1 + IDX_ANA0_H, IDX_ANA0_V = 2, 3 def __init__(self, SPICE_CONVENTION: bool = True) -> None: """Load instrument configuration from json if provided""" @@ -50,8 +47,7 @@ def __init__(self, SPICE_CONVENTION: bool = True) -> None: def validate_instrument_parameters(self): """Check if enough instrument parameters are provided for Cooper-Nathans mehtod""" - # monochromator - try: + try: # monochromator mono = self.monochromator except AttributeError: print("Monochromator info are missing.") @@ -61,8 +57,7 @@ def validate_instrument_parameters(self): elif not all(val > 0 for val in mono_mosaic): raise ValueError("Mosaic of monochromator cannot be negative.") - # analyzer - try: + try: # analyzer ana = self.analyzer except AttributeError: print("Analyzer info are missing.") @@ -165,106 +160,116 @@ def calc_mat_c(theta_m, theta_a): def cooper_nathans( self, - hkl: tuple[float, float, float], - ei: float, - ef: float, + hkl_list: Union[tuple[float, float, float], list[tuple[float, float, float]]], + ei: Union[float, list[float]], + ef: Union[float, list[float]], projection: tuple = ((1, 0, 0), (0, 1, 0), (0, 0, 1)), R0: bool = False, - ): + ) -> ResoEllipsoid: """Calculate resolution using Cooper-Nathans method Args: - hkl : momentum transfer, miller indices in reciprocal lattice + hkl: momentum transfer, miller indices in reciprocal lattice ei: incident energy, in units of meV - ef : final energy, in units of meV - R0 : calculate normalization factor if True + ef: final energy, in units of meV + R0: calculate normalization factor if True projection (tuple): three non-coplaner vectors. If projection is None, the calculation is done in local Q frame """ self.validate_instrument_parameters() - if not (isinstance(hkl, tuple | list) and len(hkl) == 3): - raise ValueError("q needs to be a tupe of size 3.") + self._mat_f = CN.calc_mat_f(self.monochromator, self.analyzer) if self._mat_f is None else self._mat_f + self._mat_g = CN.calc_mat_g(self.collimators) if self._mat_g is None else self._mat_g - if self._mat_f is None: - self._mat_f = CN.calc_mat_f(self.monochromator, self.analyzer) + if not isinstance(hkl_list, list): + hkl_list = [hkl_list] + # TODO + else: # generate ei, ef list + pass - if self._mat_g is None: - self._mat_g = CN.calc_mat_g(self.collimators) + rez_list = [] + for hlk in hkl_list: - # q_lab = conv_mat @ hkl - # q_mod = np.linalg.norm(q_lab) - q_mod = np.linalg.norm(self.sample.b_mat @ hkl) * 2 * np.pi + # q_lab = conv_mat @ hkl + # q_mod = np.linalg.norm(q_lab) + q_mod = np.linalg.norm(self.sample.b_mat @ hkl) * 2 * np.pi - ki = np.sqrt(ei / ksq2eng) - kf = np.sqrt(ef / ksq2eng) + ki = np.sqrt(ei / ksq2eng) + kf = np.sqrt(ef / ksq2eng) - rez = ResoEllipsoid(hkl=hkl, en=ei - ef, projection=projection, sample=self.sample) + rez = ResoEllipsoid(hkle=hkl + (ei - ef,), projection=projection, sample=self.sample) - try: - two_theta = get_angle_from_triangle(ki, kf, q_mod) * self.goniometer.sense - except TypeError: - rez.STATUS = False - print(f"Cannot close triangle for ei={ei}, ef={ef}, hkl={np.round(hkl,3)}.") - return rez + try: + two_theta = get_angle_from_triangle(ki, kf, q_mod) * self.goniometer.sense + except TypeError: + rez.STATUS = False + print(f"Cannot close triangle for ei={ei}, ef={ef}, hkl={np.round(hkl,3)}.") + return rez - # phi = , always has the oppositie sign of s2 - phi = get_angle_from_triangle(ki, q_mod, kf) * self.goniometer.sense * (-1) + # phi = , always has the oppositie sign of s2 + phi = get_angle_from_triangle(ki, q_mod, kf) * self.goniometer.sense * (-1) - theta_m = get_angle_bragg(ki, self.monochromator.d_spacing) * self.monochromator.sense - theta_a = get_angle_bragg(kf, self.analyzer.d_spacing) * self.analyzer.sense + theta_m = get_angle_bragg(ki, self.monochromator.d_spacing) * self.monochromator.sense + theta_a = get_angle_bragg(kf, self.analyzer.d_spacing) * self.analyzer.sense - # TODO - # curved monochromator and analyzer + # TODO + # curved monochromator and analyzer - # TODO - # reflection efficiency - - mat_a = CN.calc_mat_a(ki, kf, theta_m, theta_a) - mat_b = CN.calc_mat_b(ki, kf, phi, two_theta) - mat_c = CN.calc_mat_c(theta_m, theta_a) - - mat_h = mat_c.T @ self._mat_f @ mat_c + self._mat_g - mat_h_inv = np.linalg.inv(mat_h) - mat_ba = mat_b @ mat_a - mat_cov = mat_ba @ mat_h_inv @ mat_ba.T - - # TODO how to add smaple mosaic in cooper-nathans? - mat_cov[1, 1] += q_mod**2 * self.sample._mosaic_h**2 - mat_cov[2, 2] += q_mod**2 * self.sample._mosaic_v**2 - - mat_reso = np.linalg.inv(mat_cov) * sig2fwhm**2 - - motor_angles = self.calculate_motor_angles(peak=hkl, ei=ei, ef=ef) - if motor_angles is None: - rez.STATUS = False - return rez - - r_mat = self.goniometer.r_mat(motor_angles) - ub_mat = self.sample.ub_mat - if self.SPICE_CONVENTION: - ub_mat = TAS.spice_to_mantid(ub_mat) - conv_mat = 2 * np.pi * r_mat @ ub_mat - rez.project_to_frame(mat_reso, phi, conv_mat) - - # TODO check normalization factor - # ------------------------------------------------------------------------- - # - if the instruments works in kf=const mode and the scans are counted for - # or normalised to monitor counts no ki^3 or kf^3 factor is needed. - # - if the instrument works in ki=const mode the kf^3 factor is needed. - - if R0: # calculate - r0 = np.pi**2 / 4 / np.sin(theta_m) / np.sin(theta_a) - r0 *= np.sqrt(np.linalg.det(self._mat_f) / np.linalg.det(mat_h)) - else: - r0 = 0 - - rez.r0 = r0 - - if np.isnan(rez.r0) or np.isinf(rez.r0) or np.isnan(rez.mat.any()) or np.isinf(rez.mat.any()): - rez.STATUS = False - else: - rez.STATUS = True - - rez.set_labels() + # TODO + # reflection efficiency + + mat_a = CN.calc_mat_a(ki, kf, theta_m, theta_a) + mat_b = CN.calc_mat_b(ki, kf, phi, two_theta) + mat_c = CN.calc_mat_c(theta_m, theta_a) + + mat_h = mat_c.T @ self._mat_f @ mat_c + self._mat_g + mat_h_inv = np.linalg.inv(mat_h) + mat_ba = mat_b @ mat_a + mat_cov = mat_ba @ mat_h_inv @ mat_ba.T + + # TODO how to add smaple mosaic in cooper-nathans? + mat_cov[1, 1] += q_mod**2 * self.sample._mosaic_h**2 + mat_cov[2, 2] += q_mod**2 * self.sample._mosaic_v**2 + + mat_reso = np.linalg.inv(mat_cov) * sig2fwhm**2 + + motor_angles = self.calculate_motor_angles(peak=hkl, ei=ei, ef=ef) + if motor_angles is None: + rez.STATUS = False + return rez + + r_mat = self.goniometer.r_mat(motor_angles) + ub_mat = self.sample.ub_mat + + if self.SPICE_CONVENTION: + ub_mat = spice_to_mantid(ub_mat) + + conv_mat = 2 * np.pi * np.matmul(r_mat, ub_mat) + rez._project_to_frame(mat_reso, phi, conv_mat) + + # TODO check normalization factor + # ------------------------------------------------------------------------- + # - if the instruments works in kf=const mode and the scans are counted for + # or normalised to monitor counts no ki^3 or kf^3 factor is needed. + # - if the instrument works in ki=const mode the kf^3 factor is needed. + + if R0: # calculate + r0 = np.pi**2 / 4 / np.sin(theta_m) / np.sin(theta_a) + r0 *= np.sqrt(np.linalg.det(self._mat_f) / np.linalg.det(mat_h)) + else: + r0 = 0 + + rez.r0 = r0 + + if np.isnan(rez.r0) or np.isinf(rez.r0) or np.isnan(rez.mat.any()) or np.isinf(rez.mat.any()): + rez.STATUS = False + else: + rez.STATUS = True + + rez.set_labels() + + rez_list.append(rez) + + if len(rez_list) == 1: + return rez[0] return rez diff --git a/src/tavi/instrument/resolution/curve.py b/src/tavi/instrument/resolution/curve.py new file mode 100644 index 00000000..59bf5992 --- /dev/null +++ b/src/tavi/instrument/resolution/curve.py @@ -0,0 +1,63 @@ +import numpy as np + +from tavi.utilities import sig2fwhm + +np.set_printoptions(floatmode="fixed", precision=4) + + +class ResoCurve(object): + """1D Gaussian curve + + Attributes: + cen (float) + fwhm (float) + xlabel + ylabel + title + legend + + Methods: + generate_curve + generate_plot + """ + + def __init__(self) -> None: + self.center: float = 0.0 + self.fwhm = None + self.r0 = None + self.x = None + self.y = None + self.xlabel = None + self.ylabel = None + self.title = None + self.legend = None + + def generate_curve(self, num_of_sigmas=3, points=100): + """Generate points on the Gaussian curve""" + + sigma = self.fwhm / sig2fwhm + cen = self.cen + self.x = np.linspace( + -num_of_sigmas * sigma + cen, + num_of_sigmas * sigma + cen, + points, + ) + amp = self.r0 / np.sqrt(2 * np.pi) / sigma + self.y = amp * np.exp(-((self.x - cen) ** 2) / sigma**2 / 2) + + def generate_plot(self, ax, c="black", linestyle="solid"): + """Generate the plot""" + self.generate_curve() + + s = ax.plot( + self.x, + self.y, + c=c, + linestyle=linestyle, + label=self.legend, + ) + ax.set_xlabel(self.xlabel) + ax.set_ylabel(self.ylabel) + ax.set_title(self.title) + ax.grid(alpha=0.6) + ax.legend() diff --git a/src/tavi/instrument/resolution/ellipse_curve.py b/src/tavi/instrument/resolution/ellipse.py similarity index 55% rename from src/tavi/instrument/resolution/ellipse_curve.py rename to src/tavi/instrument/resolution/ellipse.py index 77f538c0..63201b41 100644 --- a/src/tavi/instrument/resolution/ellipse_curve.py +++ b/src/tavi/instrument/resolution/ellipse.py @@ -1,72 +1,16 @@ import matplotlib.pyplot as plt import numpy as np +import numpy.linalg as la from mpl_toolkits.axisartist import Subplot from mpl_toolkits.axisartist.grid_finder import MaxNLocator from mpl_toolkits.axisartist.grid_helper_curvelinear import GridHelperCurveLinear +from tavi.data.scan_data import ScanData1D from tavi.utilities import sig2fwhm np.set_printoptions(floatmode="fixed", precision=4) -class ResoCurve(object): - """1D Gaussian curve - - Attributes: - cen (float) - fwhm (float) - xlabel - ylabel - title - legend - - Methods: - generate_curve - generate_plot - """ - - def __init__(self) -> None: - self.center: float = 0.0 - self.fwhm = None - self.r0 = None - self.x = None - self.y = None - self.xlabel = None - self.ylabel = None - self.title = None - self.legend = None - - def generate_curve(self, num_of_sigmas=3, points=100): - """Generate points on the Gaussian curve""" - - sigma = self.fwhm / sig2fwhm - cen = self.cen - self.x = np.linspace( - -num_of_sigmas * sigma + cen, - num_of_sigmas * sigma + cen, - points, - ) - amp = self.r0 / np.sqrt(2 * np.pi) / sigma - self.y = amp * np.exp(-((self.x - cen) ** 2) / sigma**2 / 2) - - def generate_plot(self, ax, c="black", linestyle="solid"): - """Generate the plot""" - self.generate_curve() - - s = ax.plot( - self.x, - self.y, - c=c, - linestyle=linestyle, - label=self.legend, - ) - ax.set_xlabel(self.xlabel) - ax.set_ylabel(self.ylabel) - ax.set_title(self.title) - ax.grid(alpha=0.6) - ax.legend() - - class ResoEllipse(object): """2D ellipses @@ -77,7 +21,7 @@ class ResoEllipse(object): vecs: eigen-vectors angle: angle between the two plotting axes, in degrees axes_labels (str) - ORIGIN (bool|None): shift the origin if True + Methods: generate_ellipse(ellipse_points=128) @@ -87,26 +31,26 @@ class ResoEllipse(object): """ - g_esp = 1e-8 + def __init__(self, mat, centers, angle, axes_labels): + self.mat = mat + self.centers = centers + self.angle = angle + self.axes_labels = axes_labels - def __init__(self): - self.mat = None - self.centers = None - self.fwhms = None - self.vecs = None - self.angle = None - self.axes_labels = None - - self.ORIGIN = None self.grid_helper = None - def generate_ellipse(self, ellipse_points=128): + evals, self.evecs = la.eig(mat) # evecs normalized already + self.fwhms = 1.0 / np.sqrt(np.abs(evals)) * sig2fwhm + + self.generate_axes() + + def get_points(self, num_points=128): """Generate points on a ellipse""" - phi = np.linspace(0, 2.0 * np.pi, ellipse_points) + phi = np.linspace(0, 2.0 * np.pi, num_points) pts = np.dot( - self.vecs, + self.evecs, np.array( [ self.fwhms[0] / 2 * np.cos(phi), @@ -114,10 +58,10 @@ def generate_ellipse(self, ellipse_points=128): ], ), ) - if self.ORIGIN: - pts[0] += self.centers[0] - pts[1] += self.centers[1] - return pts + + pts[0] += self.centers[0] + pts[1] += self.centers[1] + return ScanData1D(pts[0], pts[1]) def _tr(self, x, y): x, y = np.asarray(x), np.asarray(y) @@ -140,7 +84,7 @@ def generate_axes(self): def generate_plot(self, ax, c="black", linestyle="solid"): """Gnerate the ellipse for plotting""" - pts = self.generate_ellipse() + pts = self.get_points() if self.grid_helper is None: @@ -164,7 +108,7 @@ def generate_plot(self, ax, c="black", linestyle="solid"): return None def plot(self): - """Plot the ellipsis.""" + """Plot the ellipses.""" fig = plt.figure() ax = Subplot(fig, 1, 1, 1, grid_helper=self.grid_helper) diff --git a/src/tavi/instrument/resolution/ellipsoid.py b/src/tavi/instrument/resolution/ellipsoid.py index baa7576e..5132239b 100755 --- a/src/tavi/instrument/resolution/ellipsoid.py +++ b/src/tavi/instrument/resolution/ellipsoid.py @@ -2,10 +2,10 @@ import matplotlib.pyplot as plt import numpy as np -import numpy.linalg as la from mpl_toolkits.axisartist import Axes -from tavi.instrument.resolution.ellipse_curve import ResoCurve, ResoEllipse +from tavi.instrument.resolution.curve import ResoCurve +from tavi.instrument.resolution.ellipse import ResoEllipse from tavi.sample.xtal import Xtal from tavi.utilities import get_angle_vec, sig2fwhm @@ -20,8 +20,7 @@ class ResoEllipsoid(object): frame (str): "q", "hkl", or "proj" projection (tuple): three non-coplanar vectors q (tuple): momentum transfer (h', k', l') in the coordinate system specified by projection - hkl (tuple): momentum transfer (h, k, l) - en (float): energy transfer + hkle (tuple): momentum transfer (h, k, l) and energy transfer agnles (tuple): angles between plotting axes mat (np.ndarray): 4 by 4 resolution matrix @@ -38,34 +37,35 @@ class ResoEllipsoid(object): def __init__( self, - hkl: tuple[float, float, float], - en: float, + hkle: tuple[float, float, float, float], sample: Xtal, projection: Optional[tuple] = ((1, 0, 0), (0, 1, 0), (0, 0, 1)), ) -> None: - self.STATUS: Optional[bool] = None - self.q: Optional[tuple[float, float, float]] = None - self.hkl: tuple[float, float, float] = hkl + self.STATUS: bool + self.q: tuple[float, float, float] + + *hkl, en = hkle + self.hkl: tuple = tuple(hkl) self.en: float = en self.projection: Optional[tuple] = projection self.angles: tuple[float, float, float] = (90.0, 90.0, 90.0) self.axes_labels = None - self.mat = None - self.r0 = None + self.mat: np.ndarray + self.r0: Optional[float] = None match self.projection: case None: # Local Q frame self.frame = "q" self.angles = (90.0, 90.0, 90.0) - q_norm = np.linalg.norm(sample.b_mat @ hkl) * 2 * np.pi + q_norm = np.linalg.norm(sample.b_mat @ self.hkl) * 2 * np.pi self.q = (q_norm, 0.0, 0.0) case ((1, 0, 0), (0, 1, 0), (0, 0, 1)): # HKL self.frame = "hkl" - self.q = hkl + self.q = self.hkl self.angles = (sample.gamma_star, sample.alpha_star, sample.beta_star) case _: # customized projection p1, p2, p3 = self.projection @@ -89,7 +89,7 @@ def __init__( self.set_labels() - def project_to_frame(self, mat_reso, phi, conv_mat): + def _project_to_frame(self, mat_reso, phi, conv_mat): """determinate the frame from the projection vectors""" mat_lab_to_local = np.array( @@ -190,7 +190,7 @@ def quadric_proj(quadric, idx): # print("\nProjected row/column %d:\n%s\n->\n%s.\n" % (idx, str(quadric), str(ortho_proj))) return np.delete(np.delete(ortho_proj, idx, axis=0), idx, axis=1) - def generate_ellipse(self, axes=(0, 1), PROJECTION=False, ORIGIN=True): + def get_ellipse(self, axes=(0, 1), PROJECTION=False, ORIGIN=True) -> ResoEllipse: """Gnerate a 2D ellipse by either making a cut or projection Arguments: @@ -199,48 +199,41 @@ def generate_ellipse(self, axes=(0, 1), PROJECTION=False, ORIGIN=True): ORIGIN: shift the center if True """ - ellipse = ResoEllipse() - ellipse.ORIGIN = ORIGIN + qe_list = np.concatenate((self.q, self.en), axis=None) # axes = np.sort(axes) # match tuple(np.sort(axes)): match tuple(axes): case (0, 1) | (1, 0): - ellipse.angle = np.round(self.angles[0], 2) + angle = np.round(self.angles[0], 2) case (1, 2) | (2, 1): - ellipse.angle = np.round(self.angles[1], 2) + angle = np.round(self.angles[1], 2) case (0, 2) | (2, 0): - ellipse.angle = np.round(self.angles[2], 2) + angle = np.round(self.angles[2], 2) case _: - ellipse.angle = 90.0 + angle = 90.0 + axes_labels = (self.axes_labels[axes[0]], self.axes_labels[axes[1]]) + + if ORIGIN: + centers = (qe_list[axes[0]], qe_list[axes[1]]) + else: + centers = (0.0, 0.0) + q_res = self.mat if PROJECTION: # make a projection - q_res = self.mat - ellipse.centers = (qe_list[axes[0]], qe_list[axes[1]]) - ellipse.axes_labels = (self.axes_labels[axes[0]], self.axes_labels[axes[1]]) # Qres_QxE_proj = np.delete(np.delete(self.mat, 2, axis=0), 2, axis=1) for i in (3, 2, 1, 0): if i not in axes: q_res = ResoEllipsoid.quadric_proj(q_res, i) - ellipse.mat = q_res + mat = q_res else: # make a slice - q_res = self.mat - ellipse.centers = (qe_list[axes[0]], qe_list[axes[1]]) - ellipse.axes_labels = (self.axes_labels[axes[0]], self.axes_labels[axes[1]]) for i in (3, 2, 1, 0): if i not in axes: q_res = np.delete(np.delete(q_res, i, axis=0), i, axis=1) - ellipse.mat = q_res + mat = q_res - evals, evecs = la.eig(ellipse.mat) # evecs normalized already - fwhms = 1.0 / np.sqrt(np.abs(evals)) * sig2fwhm - # angles = np.array([np.arctan2(evecs[1][0], evecs[0][0])]) - ellipse.fwhms = fwhms - ellipse.vecs = evecs - ellipse.generate_axes() - - return ellipse + return ResoEllipse(mat, centers, angle, axes_labels) def set_labels(self): """Set axes labels based on the frame""" @@ -257,77 +250,45 @@ def set_labels(self): "E (meV)", ) - # TODO - def calc_ellipses(self, verbose=True): - """Calculate FWHMs and rotation angles""" - - # 2d sliced ellipses - - elps_qx_en = self.generate_ellipse(axes=(0, 3), PROJECTION=False) - # elps_qx_en.plot() - - elps_qy_en = self.generate_ellipse(axes=(1, 3), PROJECTION=False) - elps_qz_en = self.generate_ellipse(axes=(2, 3), PROJECTION=False) - - elps_qx_qy = self.generate_ellipse(axes=(0, 1), PROJECTION=False) - elps_qx_qy.plot() - elps_qy_qz = self.generate_ellipse(axes=(1, 2), PROJECTION=False) - elps_qx_qz = self.generate_ellipse(axes=(0, 2), PROJECTION=False) - elps_qx_qz.plot() - plt.show() - - # 2d projected ellipses - # Qres_QxE_proj = np.delete(np.delete(self.mat, 2, axis=0), 2, axis=1) - - elps_proj_qx_en = self.generate_ellipse(axes=(0, 3), PROJECTION=True) - elps_proj_qy_en = self.generate_ellipse(axes=(1, 3), PROJECTION=True) - elps_proj_qz_en = self.generate_ellipse(axes=(2, 3), PROJECTION=True) - - elps_proj_qx_qy = self.generate_ellipse(axes=(0, 1), PROJECTION=True) - elps_proj_qy_qz = self.generate_ellipse(axes=(1, 2), PROJECTION=True) - elps_proj_qx_qz = self.generate_ellipse(axes=(0, 2), PROJECTION=True) - - return None - def plot(self): """Plot all 2D ellipses""" # fig = plt.figure() fig = plt.figure(figsize=(10, 6)) - elps_qx_en = self.generate_ellipse(axes=(0, 3), PROJECTION=False) + elps_qx_en = self.get_ellipse(axes=(0, 3), PROJECTION=False) ax = fig.add_subplot(231, axes_class=Axes, grid_helper=elps_qx_en.grid_helper) elps_qx_en.generate_plot(ax, c="black", linestyle="solid") - elps_proj_qx_en = self.generate_ellipse(axes=(0, 3), PROJECTION=True) + elps_proj_qx_en = self.get_ellipse(axes=(0, 3), PROJECTION=True) elps_proj_qx_en.generate_plot(ax, c="black", linestyle="dashed") - elps_qy_en = self.generate_ellipse(axes=(1, 3), PROJECTION=False) + elps_qy_en = self.get_ellipse(axes=(1, 3), PROJECTION=False) ax = fig.add_subplot(232, axes_class=Axes, grid_helper=elps_qy_en.grid_helper) elps_qy_en.generate_plot(ax, c="black", linestyle="solid") - elps_proj_qy_en = self.generate_ellipse(axes=(1, 3), PROJECTION=True) + elps_proj_qy_en = self.get_ellipse(axes=(1, 3), PROJECTION=True) elps_proj_qy_en.generate_plot(ax, c="black", linestyle="dashed") - elps_qz_en = self.generate_ellipse(axes=(2, 3), PROJECTION=False) + elps_qz_en = self.get_ellipse(axes=(2, 3), PROJECTION=False) ax = fig.add_subplot(233, axes_class=Axes, grid_helper=elps_qz_en.grid_helper) elps_qz_en.generate_plot(ax, c="black", linestyle="solid") - elps_proj_qz_en = self.generate_ellipse(axes=(2, 3), PROJECTION=True) + elps_proj_qz_en = self.get_ellipse(axes=(2, 3), PROJECTION=True) elps_proj_qz_en.generate_plot(ax, c="black", linestyle="dashed") - elps_qx_qy = self.generate_ellipse(axes=(0, 1), PROJECTION=False) + elps_qx_qy = self.get_ellipse(axes=(0, 1), PROJECTION=False) ax = fig.add_subplot(234, axes_class=Axes, grid_helper=elps_qx_qy.grid_helper) elps_qx_qy.generate_plot(ax, c="black", linestyle="solid") - elps_proj_qx_qy = self.generate_ellipse(axes=(0, 1), PROJECTION=True) + elps_proj_qx_qy = self.get_ellipse(axes=(0, 1), PROJECTION=True) elps_proj_qx_qy.generate_plot(ax, c="black", linestyle="dashed") - elps_qy_qz = self.generate_ellipse(axes=(1, 2), PROJECTION=False) + elps_qy_qz = self.get_ellipse(axes=(1, 2), PROJECTION=False) ax = fig.add_subplot(235, axes_class=Axes, grid_helper=elps_qy_qz.grid_helper) elps_qy_qz.generate_plot(ax, c="black", linestyle="solid") - elps_proj_qy_qz = self.generate_ellipse(axes=(1, 2), PROJECTION=True) + elps_proj_qy_qz = self.get_ellipse(axes=(1, 2), PROJECTION=True) elps_proj_qy_qz.generate_plot(ax, c="black", linestyle="dashed") - elps_qx_qz = self.generate_ellipse(axes=(0, 2), PROJECTION=False) + elps_qx_qz = self.get_ellipse(axes=(0, 2), PROJECTION=False) ax = fig.add_subplot(236, axes_class=Axes, grid_helper=elps_qx_qz.grid_helper) elps_qx_qz.generate_plot(ax, c="black", linestyle="solid") - elps_proj_qx_qz = self.generate_ellipse(axes=(0, 2), PROJECTION=True) + elps_proj_qx_qz = self.get_ellipse(axes=(0, 2), PROJECTION=True) elps_proj_qx_qz.generate_plot(ax, c="black", linestyle="dashed") fig.tight_layout(pad=2) diff --git a/src/tavi/instrument/tas.py b/src/tavi/instrument/tas.py index 0f7ddcee..d0824a2b 100644 --- a/src/tavi/instrument/tas.py +++ b/src/tavi/instrument/tas.py @@ -23,11 +23,7 @@ def __init__(self, SPICE_CONVENTION: bool = True): self.SPICE_CONVENTION = SPICE_CONVENTION # use coordination system defined in SPICE @staticmethod - def q_lab( - two_theta_deg: float, - ei: float, - ef: float, - ): + def q_lab(two_theta_deg: float, ei: float, ef: float): """ Reutrn momentum transfer q in lab frame, using Mantid convention diff --git a/src/tavi/instrument/tas_base.py b/src/tavi/instrument/tas_base.py index 4cc5e14e..34d2c912 100644 --- a/src/tavi/instrument/tas_base.py +++ b/src/tavi/instrument/tas_base.py @@ -4,9 +4,14 @@ from typing import Optional, Union from tavi.data.scan import Scan -from tavi.instrument.goni import Goniometer -from tavi.instrument.mono_ana import MonoAna -from tavi.instrument.tas_cmponents import Collimators, Detector, Distances, Guide, Monitor, Source +from tavi.instrument.components.collimators import Collimators +from tavi.instrument.components.component_base import Distances +from tavi.instrument.components.detector import Detector +from tavi.instrument.components.goni import Goniometer +from tavi.instrument.components.guide import Guide +from tavi.instrument.components.monitor import Monitor +from tavi.instrument.components.mono_ana import MonoAna +from tavi.instrument.components.source import Source from tavi.sample.powder import Powder from tavi.sample.sample import Sample from tavi.sample.xtal import Xtal @@ -32,16 +37,18 @@ class TASBase(object): """ def __init__(self) -> None: - """Load instrument configuration from json if provided, otherwise leace as None""" + """Load instrument configuration from json if provided, otherwise leave as None""" self.monochromator: MonoAna + self.analyzer: MonoAna self.goniometer: Goniometer + self.collimators: Collimators self.sample: Union[Sample, Xtal, Powder] self.source: Optional[Source] = None - self.collimators: Optional[Collimators] = None + self.guide: Optional[Guide] = None self.monitor: Optional[Monitor] = None - self.analyzer: Optional[MonoAna] = None + self.detector: Optional[Detector] = None self.arms: Optional[Distances] = None diff --git a/src/tavi/instrument/tas_cmponents.py b/src/tavi/instrument/tas_cmponents.py deleted file mode 100644 index 9d2ed8fb..00000000 --- a/src/tavi/instrument/tas_cmponents.py +++ /dev/null @@ -1,301 +0,0 @@ -# -*- coding: utf-8 -*- -from typing import Literal, Optional, Union - -import numpy as np - -from tavi.utilities import cm2angstrom - - -class TASComponent(object): - """ - A component of the triple-axis spectrometer - - Attributes: - type (str): identify which component - Methods: - update(param, value): update a parameter - """ - - def __init__( - self, - param_dict: Optional[dict] = None, - component_name: str = "", - ): - """Load parameters if provided.""" - self.component_name = component_name - if param_dict is not None: - for key, val in param_dict.items(): - setattr(self, key, val) - - def update( - self, - param: str, - value: Union[str, float, int], - ) -> None: - """update a paramter if exist""" - if hasattr(self, param): - setattr(self, param, value) - print(f"Setting {self.component_name} parameter {param} to {value}") - else: - print(f"{self.component_name} has no attribute {param}") - - @staticmethod - def _min2rad( - angle_in_minutes: Optional[float], - param_name: str = "Angle", - ) -> Optional[float]: - """Convert from minutes to radian is angle is not None""" - if angle_in_minutes is None: - print(f"{param_name} is None.") - return None - return np.deg2rad(angle_in_minutes / 60) - - @staticmethod - def _cm2angstrom_given_shape( - length_in_cm: Optional[float], - shape: str, - param_name: str = "Length", - ) -> Optional[float]: - """ - Convert length from centimeter to angstrom. - - Apply correction based on shape, for resolution calculation. - Divide length by np.sqrt(12) if rectangular, - Divive diameter D by 4 if circular or spherical. - """ - if length_in_cm is None: - print(f"{param_name} is None.") - return None - match shape: - case "rectangular": - return length_in_cm / np.sqrt(12) * cm2angstrom - case "circular" | "spherical": - return length_in_cm / 4 * cm2angstrom - case _: - print("Unrecognized shape. Needs to be rectangular or circular/spherical.") - return None - - @staticmethod - def _cm2angstrom( - length_in_cm: Optional[float], - param_name: str = "Length", - ) -> Optional[float]: - """ - Convert length from centimeter to angstrom. - """ - if length_in_cm is None: - print(f"{param_name} is None.") - return None - - return length_in_cm * cm2angstrom - - -class Source(TASComponent): - """Neutron source""" - - def __init__( - self, - param_dict: Optional[dict] = None, - component_name="source", - ): - self.name: str = "HFIR" - self.shape: Literal["rectangular", "circular"] = "rectangular" - self.width: Optional[float] = None - self.height: Optional[float] = None - - super().__init__(param_dict, component_name) - - @property - def _width(self): - """width in angstrom, with correction based on shape, for resolution calculation""" - return TASComponent._cm2angstrom_given_shape(self.width, self.shape, "Source width") - - @property - def _height(self): - """Height in angstrom, with correction based on shape, for resolution calculation""" - return TASComponent._cm2angstrom_given_shape(self.height, self.shape, "Source height") - - -class Guide(TASComponent): - """Neutron guide in the monochromator drum""" - - def __init__( - self, - param_dict: Optional[dict] = None, - component_name: str = "guide", - ): - self.in_use: bool = False - self.div_h: Optional[float] = None - self.div_v: Optional[float] = None - - super().__init__(param_dict, component_name) - - @property - def _div_h(self): - """Horizontal divergence in radian""" - return TASComponent._min2rad(self.div_h, "Guide horizontal divergence") - - @property - def _div_v(self): - """Vertical divergence in radian""" - return TASComponent._min2rad(self.div_v, "Guide vertical divergence") - - -class Collimators(TASComponent): - """collimitor divergences, in mins of arc""" - - def __init__( - self, - param_dict: Optional[dict] = None, - component_name: str = "collimitors", - ): - # defalut values - self.h_pre_mono: float = 30.0 # mins of arc - self.h_pre_sample: float = 30.0 - self.h_post_sample: float = 30.0 - self.h_post_ana: float = 30.0 - self.v_pre_mono: float = 30.0 - self.v_pre_sample: float = 30.0 - self.v_post_sample: float = 30.0 - self.v_post_ana: float = 30.0 - - super().__init__(param_dict, component_name) - - @property - def horizontal_divergence(self) -> list: - """list of horizontal divergence in minitus of arc""" - return [ - self.h_pre_mono, - self.h_pre_sample, - self.h_post_sample, - self.h_post_ana, - ] - - @property - def _horizontal_divergence(self) -> list: - """list of horizontal divergence in radian""" - return [ - np.deg2rad(self.h_pre_mono / 60), - np.deg2rad(self.h_pre_sample / 60), - np.deg2rad(self.h_post_sample / 60), - np.deg2rad(self.h_post_ana / 60), - ] - - @property - def vertical_divergence(self) -> list: - """list of vertical divergence in minutes of arcs""" - return [ - self.v_pre_mono, - self.v_pre_sample, - self.v_post_sample, - self.v_post_ana, - ] - - @property - def _vertical_divergence(self) -> list: - """list of vertical divergence in radian""" - return [ - np.deg2rad(self.v_pre_mono / 60), - np.deg2rad(self.v_pre_sample / 60), - np.deg2rad(self.v_post_sample / 60), - np.deg2rad(self.v_post_ana / 60), - ] - - -# TODO not implemented in resolution calculation -class Monitor(TASComponent): - """ - Neutron Monitor - - Note: - Needs to think about monitor efficiency - """ - - def __init__( - self, - param_dict: Optional[dict] = None, - component_name: str = "monitor", - ): - self.shape: Literal["rectangular", "spherical"] = "rectangular" - self.width: Optional[float] = None - self.height: Optional[float] = None - - super().__init__(param_dict, component_name) - - @property - def _width(self): - """width in angstrom, with correction based on shape, for resolution calculation""" - return TASComponent._cm2angstrom_given_shape(self.width, self.shape, "Monitor width") - - @property - def _height(self): - """Height in angstrom, with correction based on shape, for resolution calculation""" - return TASComponent._cm2angstrom_given_shape(self.height, self.shape, "Monitor height") - - -# TODO not implemented in resolution calculation -class Detector(TASComponent): - """ - Neutron detector - - Note: - Need to think about detector efficiency, saturation - """ - - def __init__( - self, - param_dict: Optional[dict] = None, - component_name: str = "detector", - ): - self.shape: Literal["rectangular", "spherical"] = "rectangular" - self.width: Optional[float] = None - self.height: Optional[float] = None - - super().__init__(param_dict, component_name) - - @property - def _width(self): - """width in angstrom, with correction based on shape, for resolution calculation""" - return TASComponent._cm2angstrom_given_shape(self.width, self.shape, "Detector width") - - @property - def _height(self): - """Height in angstrom, with correction based on shape, for resolution calculation""" - return TASComponent._cm2angstrom_given_shape(self.height, self.shape, "Detector height") - - -class Distances(TASComponent): - """Distance between components, in units of centimeters""" - - def __init__( - self, - param_dict: Optional[dict] = None, - component_name: str = "distances", - ): - self.src_mono: Optional[float] = None - self.mono_sample: Optional[float] = None - self.sample_ana: Optional[float] = None - self.ana_det: Optional[float] = None - self.mono_monitor: Optional[float] = None - - super().__init__(param_dict, component_name) - - @property - def _src_mono(self): - return TASComponent._cm2angstrom(self.src_mono, "Source to monochromator distance") - - @property - def _mono_sample(self): - return TASComponent._cm2angstrom(self.src_mono, "Monochromator to sample distance") - - @property - def _sample_ana(self): - return TASComponent._cm2angstrom(self.src_mono, "Sample to analyzer distance") - - @property - def _ana_det(self): - return TASComponent._cm2angstrom(self.src_mono, "Analyzer to detector distance") - - @property - def _mono_monitor(self): - return TASComponent._cm2angstrom(self.src_mono, "Monochromator to monitor distance") diff --git a/tests/test_cooper_nathans.py b/tests/test_cooper_nathans.py index 504b9993..f2908831 100755 --- a/tests/test_cooper_nathans.py +++ b/tests/test_cooper_nathans.py @@ -10,7 +10,7 @@ def test_copper_nathans_localQ(tas_params): tas, ei, ef, hkl, _, r0 = tas_params - rez = tas.cooper_nathans(ei=ei, ef=ef, hkl=hkl, projection=None, R0=r0) + rez = tas.cooper_nathans(hkl_list=hkl, ei=ei, ef=ef, projection=None, R0=r0) mat = np.array( [ [9583.2881, -4671.0614, -0.0000, 986.5610], @@ -25,7 +25,7 @@ def test_copper_nathans_localQ(tas_params): def test_copper_nathans_hkl(tas_params): tas, ei, ef, hkl, _, r0 = tas_params - rez = tas.cooper_nathans(ei=ei, ef=ef, hkl=hkl, R0=r0) + rez = tas.cooper_nathans(hkl_list=hkl, ei=ei, ef=ef, R0=r0) mat = np.array( [ [33305.0843, 33224.4963, -2651.8290, -5152.9962], @@ -40,7 +40,7 @@ def test_copper_nathans_hkl(tas_params): def test_copper_nathans_projection(tas_params): tas, ei, ef, hkl, projection, r0 = tas_params - rez = tas.cooper_nathans(ei=ei, ef=ef, hkl=hkl, projection=projection, R0=r0) + rez = tas.cooper_nathans(hkl_list=hkl, ei=ei, ef=ef, projection=projection, R0=r0) mat = np.array( [ [1.3306e05, -5.3037e03, -1.7660e-01, -1.0306e04], @@ -52,6 +52,22 @@ def test_copper_nathans_projection(tas_params): assert np.allclose(rez.mat, mat, atol=1e-1) +def test_copper_nathans_list(tas_params): + tas, ei, ef, _, _, r0 = tas_params + + rez = tas.cooper_nathans(hkl_list=[(0, 0, 3), (0, 0, -3)], ei=ei, ef=ef, projection=None, R0=r0) + mat = np.array( + [ + [9583.2881, -4671.0614, -0.0000, 986.5610], + [-4671.0614, 21359.2992, 0.0000, -4129.1553], + [0.0000, 0.0000, 77.7036, 0.0000], + [986.5610, -4129.1553, -0.0000, 864.3494], + ] + ) + assert len(rez) == 2 + assert np.allclose(rez[0].mat, mat, atol=1e-1) + + @pytest.fixture def tas_params(): # cooper_nathans_CTAX diff --git a/tests/test_ellipse.py b/tests/test_ellipse.py new file mode 100644 index 00000000..dac52759 --- /dev/null +++ b/tests/test_ellipse.py @@ -0,0 +1,61 @@ +import matplotlib.pyplot as plt +import numpy as np +import pytest +from mpl_toolkits.axisartist import Subplot + +from tavi.instrument.resolution.cooper_nathans import CN +from tavi.plotter import Plot2D +from tavi.sample.xtal import Xtal + +np.set_printoptions(floatmode="fixed", precision=4) + + +def test_local_q(tas_params): + tas, ei, ef, hkl, _, R0 = tas_params + rez = tas.cooper_nathans(ei=ei, ef=ef, hkl=hkl, projection=None, R0=R0) + ellipse = rez.get_ellipse(axes=(0, 3), PROJECTION=False) + + assert np.allclose(ellipse.angle, 90) + assert ellipse.axes_labels == ("Q_para (1/A)", "E (meV)") + + +def test_hkl(tas_params): + tas, ei, ef, hkl, _, R0 = tas_params + rez = tas.cooper_nathans(ei=ei, ef=ef, hkl=hkl, R0=R0) + ellipse = rez.get_ellipse(axes=(0, 1), PROJECTION=False) + + assert np.allclose(ellipse.angle, 60) + assert ellipse.axes_labels == ("H (r.l.u.)", "K (r.l.u.)") + + p1 = Plot2D() + p1.add_curve(ellipse.get_points(), fmt="-k") + + fig = plt.figure() + ax = Subplot(fig, 1, 1, 1) + fig.add_subplot(ax) + p1.plot(ax) + plt.show() + + +@pytest.fixture +def tas_params(): + # cooper_nathans_CTAX + + instrument_config_json_path = "./src/tavi/instrument/instrument_params/cg4c.json" + tas = CN(SPICE_CONVENTION=False) + tas.load_instrument_params_from_json(instrument_config_json_path) + + sample_json_path = "./test_data/test_samples/nitio3.json" + sample = Xtal.from_json(sample_json_path) + tas.mount_sample(sample) + + ei = 4.8 + ef = 4.8 + hkl = (0, 0, 3) + + projection = ((1, 1, 0), (0, 0, 1), (1, -1, 0)) + R0 = False + + tas_params = (tas, ei, ef, hkl, projection, R0) + + return tas_params diff --git a/tests/test_ellipsoid.py b/tests/test_ellipsoid.py new file mode 100644 index 00000000..97681a2b --- /dev/null +++ b/tests/test_ellipsoid.py @@ -0,0 +1,73 @@ +import matplotlib.pyplot as plt +import numpy as np +import pytest + +from tavi.instrument.resolution.cooper_nathans import CN +from tavi.sample.xtal import Xtal + +np.set_printoptions(floatmode="fixed", precision=4) + + +def test_local_q(tas_params): + tas, ei, ef, hkl, _, R0 = tas_params + rez = tas.cooper_nathans(ei=ei, ef=ef, hkl=hkl, projection=None, R0=R0) + + assert np.allclose(rez.hkl, (0, 0, 3)) + assert np.allclose(rez.q, (3 * 2 * np.pi / tas.sample.c, 0, 0)) + assert rez.projection is None + assert np.allclose(rez.angles, (90, 90, 90)) + assert rez.axes_labels == ("Q_para (1/A)", "Q_perp (1/A)", "Q_up (1/A)", "E (meV)") + + +def test_hkl(tas_params): + tas, ei, ef, hkl, _, R0 = tas_params + rez = tas.cooper_nathans(ei=ei, ef=ef, hkl=hkl, R0=R0) + + assert np.allclose(rez.hkl, (0, 0, 3)) + assert np.allclose(rez.q, (0, 0, 3)) + assert rez.projection == ((1, 0, 0), (0, 1, 0), (0, 0, 1)) + assert np.allclose(rez.angles, (60, 90, 90)) + assert rez.axes_labels == ("H (r.l.u.)", "K (r.l.u.)", "L (r.l.u.)", "E (meV)") + + +def test_projection(tas_params): + tas, ei, ef, hkl, projection, R0 = tas_params + rez = tas.cooper_nathans(ei=ei, ef=ef, hkl=hkl, projection=projection, R0=R0) + + assert np.allclose(rez.hkl, (0, 0, 3)) + assert np.allclose(rez.q, (0, 3, 0)) + assert rez.projection == ((1, 1, 0), (0, 0, 1), (1, -1, 0)) + assert np.allclose(rez.angles, (90, 90, 90)) + assert rez.axes_labels == ("(1, 1, 0)", "(0, 0, 1)", "(1, -1, 0)", "E (meV)") + + +def test_plot(tas_params): + tas, ei, ef, hkl, _, R0 = tas_params + rez = tas.cooper_nathans(ei=ei, ef=ef, hkl=hkl, R0=R0) + rez.plot() + + plt.show() + + +@pytest.fixture +def tas_params(): + # cooper_nathans_CTAX + + instrument_config_json_path = "./src/tavi/instrument/instrument_params/cg4c.json" + tas = CN(SPICE_CONVENTION=False) + tas.load_instrument_params_from_json(instrument_config_json_path) + + sample_json_path = "./test_data/test_samples/nitio3.json" + nitio3 = Xtal.from_json(sample_json_path) + tas.mount_sample(nitio3) + + ei = 4.8 + ef = 4.8 + hkl = (0, 0, 3) + + projection = ((1, 1, 0), (0, 0, 1), (1, -1, 0)) + R0 = False + + tas_params = (tas, ei, ef, hkl, projection, R0) + + return tas_params