From 1b32e010e03a9441ce797a26276b7b36aacbac42 Mon Sep 17 00:00:00 2001 From: Bing Li <104096768+bingli621@users.noreply.github.com> Date: Thu, 22 Aug 2024 09:34:21 -0400 Subject: [PATCH] modified tas component classes --- .DS_Store | Bin 10244 -> 10244 bytes scripts/example_cooper_nathans.py | 5 +- scripts/test_plotter.py | 2 +- scripts/test_rez.py | 2 +- src/tavi/instrument/goni.py | 158 +++++++ .../instrument/instrument_params/cg4c.json | 4 +- .../instrument/instrument_params/hb1a.json | 72 ++++ .../instrument/instrument_params/hb3.json | 4 +- .../instrument_params/python_dicts/cg4c.py | 105 ----- .../instrument_params/python_dicts/hb3.py | 106 ----- .../python_dicts/takin_test.py | 112 ----- .../instrument_params/takin_test.json | 4 +- src/tavi/instrument/mono_ana.py | 145 +++++++ .../resolution/backups/cooper_nathans_bak.py | 196 --------- .../resolution/backups/reso_ellipses_bak.py | 374 ----------------- .../instrument/resolution/cooper_nathans.py | 23 +- src/tavi/instrument/tas.py | 120 ++++-- src/tavi/instrument/tas_cmponents.py | 395 ++++++++---------- src/tavi/sample/powder.py | 16 - src/tavi/sample/sample.py | 56 ++- src/tavi/sample/xtal.py | 69 +-- src/tavi/utilities.py | 87 ---- test_data/nexus_exp1031.h5 | Bin 10594272 -> 10594272 bytes test_data/test_samples/sample.json | 4 +- tests/test_cooper_nathans.py | 6 +- tests/test_sample.py | 16 +- tests/test_tas.py | 19 +- 27 files changed, 743 insertions(+), 1357 deletions(-) create mode 100644 src/tavi/instrument/goni.py create mode 100644 src/tavi/instrument/instrument_params/hb1a.json delete mode 100644 src/tavi/instrument/instrument_params/python_dicts/cg4c.py delete mode 100644 src/tavi/instrument/instrument_params/python_dicts/hb3.py delete mode 100644 src/tavi/instrument/instrument_params/python_dicts/takin_test.py create mode 100644 src/tavi/instrument/mono_ana.py delete mode 100755 src/tavi/instrument/resolution/backups/cooper_nathans_bak.py delete mode 100755 src/tavi/instrument/resolution/backups/reso_ellipses_bak.py diff --git a/.DS_Store b/.DS_Store index 9b6a3a630946a21b63916b413a90eb7fec1bf222..528ecb2450932f92f6e95f7e93b9ea2d1a4fa09d 100644 GIT binary patch delta 1150 zcmdVZTS!zv7zgn0Ki>A>q{lkD9xt()m87e!Y32YA=ebhC>u z5mI{SLR1hvbdO+Dg0d(`qRT^|K}7~Vgh4@6=t~bxd)5bg=_P{ZW#<3QJbW`hrr*); zxH+Ok<>nRGQgZY4QRM1^5U6QdpudsAYWsPGzR*sv`6|QGns=)P~7W;aht}c^gx5?@E?$V5X8_$+x zzvOK5`aHpkn!|HVvPY_KceiNmMk${cnB;1yp;znndpmtbv6$;KmXm#D$wnV(Wf`Hx zY+p{FyF*JWvm03=6D_0_R8BS2Ox?7b2Iv%>r(5)h#%O{j=>vVDX#flnh(t6Lm=O;Z zsmMSkvXG5pl%oP`QHRZNVGFjx(}xcDupM3KMIZL#AP(U$hM?m(PT(9a;u5an8isKX z_wfJ^@dB^$25&KmkNApj_|76(ER$IrOJa-JDptsfSTQSM6|9yuFc;H|{ux5|a?Rc^ zIEN8y)l(FUenpK&AY^B)tXfxHGi$f#&vtrSf>#(#BxE&ls~jPo-!!f!Bql9Q7co*i zy%95QR;$e9VUrQ41Pbf&m=GF%04~&fshi z>v9n5I!16CcW@V@7{g;c!BafLbG*bmOyM)8e z>aiw%i1~$+kw8 zYswQu5ARldJ7m3YWlIE65!_pRfv!l!RK$vTu}fIz_jqN$UMgeDc(Flf3F%Mi#T87A zPZ+8NHM5q8AGKu^W4^ z7sEJ&Q5?l_oWN;};{q<>GH&Am9^w(E@B(k}7Vj{FpB!=E7$wIrl(TaV z&Z&>j!B}HsbS>U!w5b*;d7fVh1-<@mFA+CpNSbTwYU>-Ef02~_B&iOk`ZTdmsf8jl zi*XtTKW$1)OV7z=YGGPdn29kOq7bQ-TFp#^kyI>aOpMhy^`5;_ODO>>joeXVHjBJ? zhv)h1==^hF3(t1DPfzJNy`~SE+YiKOS}Cv~2UbnVjw(1%gH;-o3r$!X!RpYcyy%1< z0gbCi1KWuKgfWO+7{LJ?jG!Hhpq<1SjNvTK;XEdB&4%l^ft$F6dw7g#yfO%)Hmjq% l=@B1)(@cz{55yiGOrD6hw8s5^-z*f}X%1eqx&HZ^{Q`Q*1#18R diff --git a/scripts/example_cooper_nathans.py b/scripts/example_cooper_nathans.py index c64e2048..04befb56 100755 --- a/scripts/example_cooper_nathans.py +++ b/scripts/example_cooper_nathans.py @@ -1,5 +1,6 @@ import matplotlib.pylab as plt import numpy as np + from tavi.instrument.resolution.cooper_nathans import CN np.set_printoptions(floatmode="fixed", precision=4) @@ -33,7 +34,7 @@ def test_cooper_nathans_compare_3(): instrument_config_json_path = "./src/tavi/instrument/instrument_params/takin_test.json" sample_json_path = "./test_data/test_samples/sample_test.json" - tas.load_instrument_from_json(instrument_config_json_path) + tas.load_instrument_params_from_json(instrument_config_json_path) tas.load_sample_from_json(sample_json_path) if tas.sample.ub_matrix is None: @@ -64,7 +65,7 @@ def test_cooper_nathans_CTAX(): instrument_config_json_path = "./src/tavi/instrument/instrument_params/cg4c.json" sample_json_path = "./test_data/test_samples/nitio3.json" - tas.load_instrument_from_json(instrument_config_json_path) + tas.load_instrument_params_from_json(instrument_config_json_path) tas.load_sample_from_json(sample_json_path) # ----------- reset UB ----------------------------------- diff --git a/scripts/test_plotter.py b/scripts/test_plotter.py index 0a1cab55..ecbd7ba8 100644 --- a/scripts/test_plotter.py +++ b/scripts/test_plotter.py @@ -68,7 +68,7 @@ def test_2d_plot(tavi, tas): instrument_config_json_path = "./src/tavi/instrument/instrument_params/cg4c.json" sample_json_path = "./test_data/test_samples/nitio3.json" - tas.load_instrument_from_json(instrument_config_json_path) + tas.load_instrument_params_from_json(instrument_config_json_path) tas.load_sample_from_json(sample_json_path) test_2d_plot(tavi, tas) diff --git a/scripts/test_rez.py b/scripts/test_rez.py index 09e096d2..467761c1 100644 --- a/scripts/test_rez.py +++ b/scripts/test_rez.py @@ -57,7 +57,7 @@ def test_1D(tas): instrument_config_json_path = "./src/tavi/instrument/instrument_params/takin_test.json" sample_json_path = "./test_data/test_samples/nitio3.json" - tas.load_instrument_from_json(instrument_config_json_path) + tas.load_instrument_params_from_json(instrument_config_json_path) tas.load_sample_from_json(sample_json_path) # test_1D(tas) diff --git a/src/tavi/instrument/goni.py b/src/tavi/instrument/goni.py new file mode 100644 index 00000000..ad1edacb --- /dev/null +++ b/src/tavi/instrument/goni.py @@ -0,0 +1,158 @@ +from typing import Literal, Optional + +import numpy as np + +from tavi.instrument.tas_cmponents import TASComponent + + +class Goniometer(TASComponent): + """Goniometer table, type = Y-ZX or YZ-X""" + + def __init__( + self, + param_dict: Optional[dict] = None, + component_name: str = "goniometer", + ): + self.type: str = "Y-ZX" # Y-mZ-X for Huber stage at HB1A and HB3 + self.sense: Literal[-1, +1] = -1 + + super().__init__(param_dict) + self.component_name = component_name + + @staticmethod + def rot_x(nu): + """rotation matrix about y-axis by angle nu + + Args: + nu (float): angle in degrees + + Note: + Using Mantid convention, beam along z, y is up, x in plane + """ + + angle = np.deg2rad(nu) + c = np.cos(angle) + s = np.sin(angle) + mat = np.array( + [ + [1, 0, 0], + [0, c, -s], + [0, s, c], + ] + ) + return mat + + @staticmethod + def rot_y(omega): + """rotation matrix about y-axis by angle omega + + Args: + omega (float): angle in degrees + + Note: + Using Mantid convention, beam along z, y is up, x in plane + """ + + angle = np.deg2rad(omega) + c = np.cos(angle) + s = np.sin(angle) + mat = np.array( + [ + [c, 0, s], + [0, 1, 0], + [-s, 0, c], + ] + ) + return mat + + @staticmethod + def rot_z(mu): + """rotation matrix about z-axis by angle mu + + Args: + mu (float): angle in degrees + + Note: + Using Mantid convention, beam along z, y is up, x in plane + """ + + angle = np.deg2rad(mu) + c = np.cos(angle) + s = np.sin(angle) + mat = np.array( + [ + [c, -s, 0], + [s, c, 0], + [0, 0, 1], + ] + ) + return mat + + def r_mat( + self, + angles_deg: tuple[float, float, float], + ): + "Goniometer rotation matrix R" + + omega, sgl, sgu = angles_deg # s2, s1, sgl, sgu + match self.type: + case "Y-ZX": # HB3 + r_mat = Goniometer.rot_y(omega) @ Goniometer.rot_z(-1 * sgl) @ Goniometer.rot_x(sgu) + case "YZ-X": # CG4C ?? + r_mat = Goniometer.rot_y(omega) @ Goniometer.rot_z(sgl) @ Goniometer.rot_x(-1 * sgu) + case _: + r_mat = None + print("Unknow goniometer type. Curruntly support Y-ZX and YZ-X") + + return r_mat + + def r_mat_inv( + self, + angles: tuple[float, float, float], + ): + """inverse of rotation matrix""" + # return np.linalg.inv(self.r_mat(angles)) + return self.r_mat(angles).T + + def angles_from_r_mat(self, r_mat): + """Calculate goniometer angles from the R matrix + + Note: + range of np.arcsin is -pi/2 to pi/2 + range of np.atan2 is -pi to pi + """ + + match self.type: + case "Y-ZX" | "YZ-X": # Y-mZ-X (s1, sgl, sgu) for HB1A and HB3, Y-Z-mX (s1, sgl, sgu) for CG4C + # sgl1 = np.arcsin(r_mat[1, 0]) * rad2deg + # sgl2 = np.arccos(np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg + sgl_rad = np.arctan2(r_mat[1, 0], np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) + sgl = np.rad2deg(sgl_rad) + + # sgu1 = np.arcsin(-r_mat[1, 2] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg + # sgu2 = np.arccos(r_mat[1, 1] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg + sgu_rad = np.arctan2( + -r_mat[1, 2] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2), + r_mat[1, 1] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2), + ) + sgu = np.rad2deg(sgu_rad) + + # omega1 = np.arcsin(-r_mat[2, 0] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg + # omega2 = np.arccos(r_mat[0, 0] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg + omega_rad = np.arctan2( + -r_mat[2, 0] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2), + r_mat[0, 0] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2), + ) + omega = np.rad2deg(omega_rad) + + match self.type: + case "Y-ZX": + angles = (omega, -1 * sgl, sgu) + case "YZ-X": + angles = (omega, sgl, -1 * sgu) + + case _: + angles = None + print("Unknow goniometer type. Curruntly support Y-ZX and YZ-X.") + + return angles diff --git a/src/tavi/instrument/instrument_params/cg4c.json b/src/tavi/instrument/instrument_params/cg4c.json index 0252decc..644b5677 100644 --- a/src/tavi/instrument/instrument_params/cg4c.json +++ b/src/tavi/instrument/instrument_params/cg4c.json @@ -11,7 +11,7 @@ }, "monochromator": { "type": "PG002", - "mosaic": 30, + "mosaic_h": 30, "mosaic_v": 30, "sense": -1, "shape": "rectangular", @@ -35,7 +35,7 @@ "analyzer": { "type": "Pg002", "d_spacing": 3.35416, - "mosaic": 90, + "mosaic_h": 90, "mosaic_v": 90, "sense": -1, "shape": "rectangular", diff --git a/src/tavi/instrument/instrument_params/hb1a.json b/src/tavi/instrument/instrument_params/hb1a.json new file mode 100644 index 00000000..eed0e724 --- /dev/null +++ b/src/tavi/instrument/instrument_params/hb1a.json @@ -0,0 +1,72 @@ +{ + "source": { + "shape": "rectangular", + "width": 0.0, + "height": 0.0 + }, + "guide": { + "in_use": false, + "div_h": 0.0, + "div_v": 0.0 + }, + "monochromator": { + "type": "PG002", + "mosaic_h": 30, + "mosaic_v": 30, + "sense": 1, + "shape": "rectangular", + "width": 10.0, + "height": 18.0, + "depth": 0.0, + "curved_h": false, + "curvh": 0.0, + "optimally_curved_h": false, + "curved_v": true, + "curvv": 0.0, + "optimally_curved_v": false + }, + "monitor": {}, + "goniometer": { + "sense": -1, + "type": "Y-ZX" + }, + "analyzer": { + "type": "Pg002", + "mosaic_h": 30, + "mosaic_v": 30, + "sense": 1, + "shape": "rectangular", + "width": 15.0, + "height": 15.6, + "depth": 0.0, + "curved_h": false, + "curvh": 0.0, + "optimally_curved_h": false, + "curved_v": true, + "curvv": 0.0, + "optimally_curved_v": false + }, + "detector": { + "shape": "rectangular", + "width": 5.08, + "height": 5.08 + }, + "distances": { + "src_mono1": 474.5, + "mono1_mono2": 227.6, + "src_mono": 474.5, + "mono_sample": 143.13, + "sample_ana": 74.0, + "ana_det": 50.0 + }, + "collimators": { + "h_pre_mono": 40, + "h_pre_sample": 40, + "h_post_sample": 40, + "h_post_ana": 80, + "v_pre_mono": 600, + "v_pre_sample": 600, + "v_post_sample": 600, + "v_post_ana": 600 + } +} \ No newline at end of file diff --git a/src/tavi/instrument/instrument_params/hb3.json b/src/tavi/instrument/instrument_params/hb3.json index 1f9a1f00..618c3e33 100644 --- a/src/tavi/instrument/instrument_params/hb3.json +++ b/src/tavi/instrument/instrument_params/hb3.json @@ -11,7 +11,7 @@ }, "monochromator": { "type": "PG002", - "mosaic": 30, + "mosaic_h": 30, "mosaic_v": 30, "sense": -1, "shape": "rectangular", @@ -37,7 +37,7 @@ "analyzer": { "type": "Pg002", "d_spacing": 3.35416, - "mosaic": 40, + "mosaic_h": 40, "mosaic_v": 40, "sense": -1, "shape": "rectangular", diff --git a/src/tavi/instrument/instrument_params/python_dicts/cg4c.py b/src/tavi/instrument/instrument_params/python_dicts/cg4c.py deleted file mode 100644 index be8518e1..00000000 --- a/src/tavi/instrument/instrument_params/python_dicts/cg4c.py +++ /dev/null @@ -1,105 +0,0 @@ -from tavi.utilities import * - -source = { - "shape": "rectangular", # rectangular or circular - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 7, # / np.sqrt(12), # in cm - "height": 15, # / np.sqrt(12), # in cm -} - - -# guide before monochromator -# guide = { -# "in_use": False, -# "div_h": 0.0, -# "div_v": 0.0, -# } - -monochromator = { - "type": "PG002", - "shape": "rectangular", - "d_spacing": mono_ana_xtal["PG002"], - "mosaic": 30, # horizontal mosaic - "mosaic_v": 30, # vertical mosaic, if anisotropic - "sense": -1, - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 7.0, # / np.sqrt(12), # in cm - "height": 15.0, # / np.sqrt(12), # in cm - "depth": 0.2, # / np.sqrt(12), # in cm - # horizontal focusing - "curvh": 0.0, # in cm^-1 - # vertical focusing - "curvv": 60.4, # at Ei = 4 meV, in cm^-1 -} - -monitor = { - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - # "width": 5 / np.sqrt(12), - # "height": 12 / np.sqrt(12), -} - -goniometer = { - "sense": +1, - "type": "Y-ZX", # s1-sgl-sgu - # CTAX's angle convention is actually YZ-X, - # but the UB calculation in SPICE is done using the convetion Y-ZX -} - -analyzer = { - "type": "Pg002", - "d_spacing": mono_ana_xtal["Pg002"], - "mosaic": 90, # horizontal mosaic - "mosaic_v": 90, # vertical mosaic, if anisotropic - "sense": -1, - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 20.0 / np.sqrt(12), - "height": 15.0 / np.sqrt(12), - "depth": 0.2 / np.sqrt(12), - # horizontal focusing - "curvh": 0.0, - # vertical focusing - "curvv": 163.2, # in cm^-1 -} -detector = { - "shape": "rectangular", # rectangular or circular - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 5 / np.sqrt(12), - "height": 10 / np.sqrt(12), -} - -distances = { - "src_mono": 530.0, # in cm - "mono_sample": 160.0, - "sample_ana": 106.0, - "ana_det": 50.0, - # "mono_monitor": 86.0, -} - -collimators = { # in units of mins of arc - "h_pre_mono": 40, - "h_pre_sample": 100, - "h_post_sample": 80, - "h_post_ana": 120, - "v_pre_mono": 600, - "v_pre_sample": 600, - "v_post_sample": 600, - "v_post_ana": 600, -} - - -cg4c_config_params = { - "source": source, - # "guide": guide, - "monochromator": monochromator, - "goniometer": goniometer, - "monitor": monitor, - "analyzer": analyzer, - "detector": detector, - "distances": distances, - "collimators": collimators, -} diff --git a/src/tavi/instrument/instrument_params/python_dicts/hb3.py b/src/tavi/instrument/instrument_params/python_dicts/hb3.py deleted file mode 100644 index 093ca322..00000000 --- a/src/tavi/instrument/instrument_params/python_dicts/hb3.py +++ /dev/null @@ -1,106 +0,0 @@ -import numpy as np - -from tavi.utilities import * - -source = { - "shape": "rectangular", # rectangular or circular - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 15, # / np.sqrt(12), # in cm - "height": 15, # / np.sqrt(12), # in cm -} - - -# guide before monochromator -guide = { - "in_use": False, - "div_h": 0.0, - "div_v": 0.0, -} - -monochromator = { - "type": "PG002", - "shape": "rectangular", - "d_spacing": mono_ana_xtal["PG002"], - "mosaic": 30, # horizontal mosaic - "mosaic_v": 30, # vertical mosaic, if anisotropic - "sense": -1, - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 7.62, # / np.sqrt(12), - "height": 10.16, # / np.sqrt(12), - "depth": 0.25, # / np.sqrt(12), - # horizontal focusing - "curvh": 0.0, - # vertical focusing - "curvv": 0.0, -} - -monitor = { - "shape": "rectangular", - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 5 / np.sqrt(12), - "height": 12 / np.sqrt(12), -} - -goniometer = { - "sense": +1, - "type": "Y-ZX", # s1-sgl-sgu -} - -analyzer = { - "type": "Pg002", - "d_spacing": mono_ana_xtal["Pg002"], - "mosaic": 40, # horizontal mosaic - "mosaic_v": 40, # vertical mosaic, if anisotropic - "sense": -1, - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 7.62, # / np.sqrt(12), - "height": 7, # / np.sqrt(12), - "depth": 0.2, # / np.sqrt(12), - # horizontal focusing - "curvh": 0.0, - # vertical focusing - "curvv": 0.0, -} -detector = { - "shape": "rectangular", # rectangular or circular - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 4, # / np.sqrt(12), - "height": 12, # / np.sqrt(12), -} - -distances = { - "src_mono": 650.0, # in cm - "mono_sample": 190.0, - "sample_ana": 160.0, - "ana_det": 60.0, - "mono_monitor": 86.0, -} - -collimators = { # in units of mins of arc - "h_pre_mono": 48, - "h_pre_sample": 40, - "h_post_sample": 40, - "h_post_ana": 120, - "v_pre_mono": 150, - "v_pre_sample": 270, - "v_post_sample": 300, - "v_post_ana": 600, -} - - -config_params = { - "source": source, - "guide": guide, - "monochromator": monochromator, - "goniometer": goniometer, - "monitor": monitor, - "analyzer": analyzer, - "detector": detector, - "distances": distances, - "collimators": collimators, -} diff --git a/src/tavi/instrument/instrument_params/python_dicts/takin_test.py b/src/tavi/instrument/instrument_params/python_dicts/takin_test.py deleted file mode 100644 index fa8fef7c..00000000 --- a/src/tavi/instrument/instrument_params/python_dicts/takin_test.py +++ /dev/null @@ -1,112 +0,0 @@ -from tavi.utilities import * - -source = { - "shape": "rectangular", # rectangular or circular - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 6.0, # in cm - "height": 12.0, # in cm -} - - -# guide before monochromator -guide = { - "in_use": False, - "div_h": 15.0, # min of arc - "div_v": 15.0, # min of arc -} - -monochromator = { - "type": "PG002", - "shape": "rectangular", - # "d_spacing": mono_ana_xtal["PG002"], - "mosaic": 45, # horizontal mosaic - "mosaic_v": 45, # vertical mosaic, if anisotropic - "sense": +1, - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 12.0, - "height": 8.0, - "depth": 0.15, - # horizontal focusing - "curved_h": False, - "curvh": 0.0, - "optimally_curved_h": False, - # vertical focusing - "curved_v": False, - "curvv": 0.0, - "optimally_curved_v": False, -} - -monitor = { - "shape": "rectangular", - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 5, # / np.sqrt(12), - "height": 12, # / np.sqrt(12), -} - -goniometer = { - "sense": -1, - "type": "Y-ZX", # s1-sgl-sgu -} - -analyzer = { - "type": "Pg002", - "d_spacing": mono_ana_xtal["Pg002"], - "mosaic": 45, # horizontal mosaic - "mosaic_v": 45, # vertical mosaic, if anisotropic - "sense": +1, - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 12.0, - "height": 8.0, - "depth": 0.3, - # horizontal focusing - "curved_h": False, - "curvh": 0.0, - "optimally_curved_h": False, - # vertical focusing - "curved_v": False, - "curvv": 0.0, - "optimally_curved_v": False, -} -detector = { - "shape": "rectangular", # rectangular or circular - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - "width": 1.5, # cm - "height": 5.0, # cm -} - -distances = { # in units of cm - "src_mono": 10.0, - "mono_sample": 200.0, - "sample_ana": 115.0, - "ana_det": 85.0, - # "mono_monitor": 86.0 , -} - -collimators = { # in units of mins of arc - "h_pre_mono": 30, - "h_pre_sample": 30, - "h_post_sample": 30, - "h_post_ana": 30, - "v_pre_mono": 30, - "v_pre_sample": 30, - "v_post_sample": 30, - "v_post_ana": 30, -} - - -instrument_params = { - "source": source, - "guide": guide, - "monochromator": monochromator, - "monitor": monitor, - "goniometer": goniometer, - "analyzer": analyzer, - "detector": detector, - "distances": distances, - "collimators": collimators, -} diff --git a/src/tavi/instrument/instrument_params/takin_test.json b/src/tavi/instrument/instrument_params/takin_test.json index 95638b4f..ec33e9d9 100644 --- a/src/tavi/instrument/instrument_params/takin_test.json +++ b/src/tavi/instrument/instrument_params/takin_test.json @@ -11,7 +11,7 @@ }, "monochromator": { "type": "PG002", - "mosaic": 45, + "mosaic_h": 45, "mosaic_v": 45, "sense": 1, "shape": "rectangular", @@ -37,7 +37,7 @@ "analyzer": { "type": "Pg002", "d_spacing": 3.35416, - "mosaic": 45, + "mosaic_h": 45, "mosaic_v": 45, "sense": 1, "shape": "rectangular", diff --git a/src/tavi/instrument/mono_ana.py b/src/tavi/instrument/mono_ana.py new file mode 100644 index 00000000..7bb0a433 --- /dev/null +++ b/src/tavi/instrument/mono_ana.py @@ -0,0 +1,145 @@ +from typing import Literal, Optional + +import numpy as np + +from tavi.instrument.tas_cmponents import TASComponent +from tavi.utilities import cm2angstrom + +# --------------------------------------------------------------- +# d_spacing table from Shirane Appendix 3, in units of Angstrom +# --------------------------------------------------------------- +mono_ana_xtal = { + "PG002": 3.35416, + "Pg002": 3.35416, + "PG004": 1.67708, + "Cu111": 2.08717, + "Cu220": 1.27813, + "Ge111": 3.26627, + "Ge220": 2.00018, + "Ge311": 1.70576, + "Ge331": 1.29789, + "Be002": 1.79160, + "Be110": 1.14280, + "Heusler": 3.435, # Cu2MnAl(111) +} + + +class MonoAna(TASComponent): + """Monochromator and analyzer class""" + + def __init__( + self, + param_dict: Optional[dict] = None, + component_name: str = "", + ): + # defalut values + self.type: str = "PG002" + self.d_spacing: float = mono_ana_xtal["PG002"] + self.mosaic_h: float = 45.0 # horizontal mosaic, min of arc + self.mosaic_v: float = 45.0 # vertical mosaic, if anisotropic + self.sense: Literal[-1, 1] = -1 # +1 for counter-clockwise, -1 for clockwise + + # divide by np.sqrt(12) if rectangular + # Diameter D/4 if spherical + self.shape: Literal["rectangular", "spherical"] = "rectangular" + self.width: float = 12.0 + self.height: float = 8.0 + self.depth: float = 0.15 + # horizontal focusing + self.curved_h: bool = False + self.curvh: float = 0.0 + self.optimally_curved_h: bool = False + # vertical focusing + self.curved_v: bool = False + self.curvv: float = 0.0 + self.optimally_curved_v: bool = False + + super().__init__(param_dict, component_name) + self.d_spacing = mono_ana_xtal[self.type] + + @property + def _mosaic_h(self): + return np.deg2rad(self.mosaic_h / 60) + + @property + def _mosaic_v(self): + return np.deg2rad(self.mosaic_v / 60) + + @property + def _width(self): + """width in angstrom, with correction based on shape, for resolution calculation""" + match self.shape: + case "rectangular": + return self.width / np.sqrt(12) * cm2angstrom + case "spherical": + return self.width / 4 * cm2angstrom + case _: + print("Unrecognized monochromator shape. Needs to be rectangular or spherical.") + return None + + @property + def _height(self): + """height in angstrom, with correction based on shape, for resolution calculation""" + match self.shape: + case "rectangular": + return self.height / np.sqrt(12) * cm2angstrom + case "spherical": + return self.height / 4 * cm2angstrom + case _: + print("Unrecognized monochromator shape. Needs to be rectangular or spherical.") + return None + + @property + def depth_ang(self): + """depth in angstrom, with correction based on shape, for resolution calculation""" + match self.shape: + case "rectangular": + return self.depth / np.sqrt(12) * cm2angstrom + case "spherical": + return self.depth / 4 * cm2angstrom + case _: + print("Unrecognized monochromator shape. Needs to be rectangular or spherical.") + return None + + +# # TODO implement curvature +# class Analyzer(object): +# def __init__(self, param_dict): +# self.type = "Pg002" +# self.d_spacing = mono_ana_xtal["Pg002"] +# self.mosaic = 45 # horizontal mosaic, min of arc +# self.mosaic_v = 45 # vertical mosaic, if anisotropic +# self.sense = -1 +# # divide by np.sqrt(12) if rectangular +# # Diameter D/4 if spherical +# self.shape = "rectangular" +# self.width = 12.0 +# self.height = 8.0 +# self.depth = 0.3 +# # horizontal focusing +# self.curved_h = False +# self.curvh = 0.0 +# self.optimally_curved_h = False +# # vertical focusing +# self.curved_v = False +# self.curvv = 0.0 +# self.optimally_curved_v = False + +# for key, val in param_dict.items(): +# match key: +# # case "mosaic" | "mosaic_v" | "curvh" | "curvv": +# # setattr(self, key, val * min2rad) +# case "width" | "height" | "depth": +# # divide by np.sqrt(12) if rectangular +# # Diameter D/4 if spherical +# if param_dict["shape"] == "rectangular": +# setattr(self, key, val / np.sqrt(12) * cm2angstrom) +# elif param_dict["shape"] == "spherical": +# setattr(self, key, val / 4 * cm2angstrom) +# else: +# print("Analyzer shape needs to be either rectangular or spherical.") + +# setattr(self, key, val * cm2angstrom) +# case _: +# setattr(self, key, val) +# self.d_spacing = mono_ana_xtal[self.type] diff --git a/src/tavi/instrument/resolution/backups/cooper_nathans_bak.py b/src/tavi/instrument/resolution/backups/cooper_nathans_bak.py deleted file mode 100755 index 53b4c756..00000000 --- a/src/tavi/instrument/resolution/backups/cooper_nathans_bak.py +++ /dev/null @@ -1,196 +0,0 @@ -import numpy as np -import numpy.linalg as la -from tavi.utilities import * -from tavi.instrument.tas import TAS -from tavi.instrument.resolution.reso_ellipses import ResoEllipsoid - - -class CN(TAS): - """Copper-Nathans method - - Attibutes: - _mat_f - _mat_g - - - Methods: - cooper_nathans: resolution in Q-local frame - cooper_nathans_hkle: resolution in (H,K,L,E) frame - - """ - - # 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 - # 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 - - def __init__(self): - super().__init__() - - # constants independent of q and eng - self._mat_f = None - self._mat_g = None - - def cooper_nathans_hkle(self, ei, ef, hkl, R0=False): - """Calculate Cooper-Nathans resolution fucntion in hkle frame""" - - angles = self.find_angles(hkl, ei, ef) - r_mat = self.goniometer.r_mat(angles[1:]) - ub_mat = self.sample.ub_matrix - conv_mat = 2 * np.pi * r_mat @ ub_mat - q_lab = conv_mat @ hkl - q = np.linalg.norm(q_lab) - rez_hkle = self.cooper_nathans(ei, ef, q, R0) - - ki = np.sqrt(ei / ksq2eng) - kf = np.sqrt(ef / ksq2eng) - - # opposite sign of s2 - phi = get_angle(ki, q, kf) * np.sign(angles[0]) * (-1) - # phi = np.deg2rad(0) - - conv_mat_4d = np.eye(4) - # conv_mat_4d[0:3, 0:3] = rot_y(phi * rad2deg) @ rot_x(90) @ conv_mat - conv_mat_4d[0:3, 0:3] = ( - np.array( - [ - [np.sin(phi), 0, np.cos(phi)], - [np.cos(phi), 0, -np.sin(phi)], - [0, 1, 0], - ] - ) - @ conv_mat - ) - - rez_hkle.mat = conv_mat_4d.T @ rez_hkle.mat @ conv_mat_4d - self.frame = "HKLE" - return rez_hkle - - def cooper_nathans(self, ei, ef, q, R0=False): - """Calculate resolution using Cooper-Nathans method - - Args: - ei (float): incident energy, in units of meV - ef (float): final energy, in units of meV - q (float | tuple of floats): momentum transfer, in units of inverse angstrom - R0 (bool): calculate normalization factor if True - """ - self.frame = "Qlocal" - if self._mat_f is None: - # matrix F, divergence of monochromator and analyzer, [pop75] Appendix 1 - mat_f = np.zeros(((CN.NUM_MONOS + CN.NUM_ANAS) * 2, (CN.NUM_MONOS + CN.NUM_ANAS) * 2)) - mat_f[CN.IDX_MONO0_H, CN.IDX_MONO0_H] = 1.0 / self.monochromator.mosaic**2 - mat_f[CN.IDX_MONO0_V, CN.IDX_MONO0_V] = 1.0 / self.monochromator.mosaic_v**2 - mat_f[CN.IDX_ANA0_H, CN.IDX_ANA0_H] = 1.0 / self.analyzer.mosaic**2 - mat_f[CN.IDX_ANA0_V, CN.IDX_ANA0_V] = 1.0 / self.analyzer.mosaic_v**2 - self._mat_f = mat_f - - if self._mat_g is None: - # matrix G, divergence of collimators, [pop75] Appendix 1 - mat_g = np.zeros((CN.NUM_COLLS * 2, CN.NUM_COLLS * 2)) - mat_g[CN.IDX_COLL0_H, CN.IDX_COLL0_H] = 1.0 / self.collimators.h_pre_mono**2 - mat_g[CN.IDX_COLL0_V, CN.IDX_COLL0_V] = 1.0 / self.collimators.v_pre_mono**2 - mat_g[CN.IDX_COLL1_H, CN.IDX_COLL1_H] = 1.0 / self.collimators.h_pre_sample**2 - mat_g[CN.IDX_COLL1_V, CN.IDX_COLL1_V] = 1.0 / self.collimators.v_pre_sample**2 - mat_g[CN.IDX_COLL2_H, CN.IDX_COLL2_H] = 1.0 / self.collimators.h_post_sample**2 - mat_g[CN.IDX_COLL2_V, CN.IDX_COLL2_V] = 1.0 / self.collimators.v_post_sample**2 - mat_g[CN.IDX_COLL3_H, CN.IDX_COLL3_H] = 1.0 / self.collimators.h_post_ana**2 - mat_g[CN.IDX_COLL3_V, CN.IDX_COLL3_V] = 1.0 / self.collimators.v_post_ana**2 - - self._mat_g = mat_g - - if isinstance(q, tuple | list): - if len(q) == 3: - h, k, l = q - q = self.sample.hkl2q((h, k, l)) - else: - print("Length of q is not 3.") - elif isinstance(q, float): - pass - else: - print("q is not float or tupe or list.") - - ki = np.sqrt(ei / ksq2eng) - kf = np.sqrt(ef / ksq2eng) - en = ei - ef - - two_theta = get_angle(ki, kf, q) * self.goniometer.sense - # theta_s = two_theta / 2 - - # phi = , always has the oppositie sign of theta_s - phi = get_angle(ki, q, 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 - - # TODO - # curved monochromator and analyzer - - # TODO - # reflection efficiency - - # TODO - # matrix A,Y=AU, tranform from collimators angular divergence to ki-kf frame - mat_a = np.zeros((6, 2 * CN.NUM_COLLS)) - mat_a[0, CN.IDX_COLL0_H] = 0.5 * ki / np.tan(theta_m) - mat_a[0, CN.IDX_COLL1_H] = -0.5 * ki / np.tan(theta_m) - mat_a[1, CN.IDX_COLL1_H] = ki - mat_a[2, CN.IDX_COLL1_V] = -ki - - mat_a[3, CN.IDX_COLL2_H] = 0.5 * kf / np.tan(theta_a) - mat_a[3, CN.IDX_COLL3_H] = -0.5 * kf / np.tan(theta_a) - mat_a[4, CN.IDX_COLL2_H] = kf - mat_a[5, CN.IDX_COLL2_V] = kf - - # matrix B, X=BY, transform from ki-kf frame to momentum transfer q-frame - mat_b = np.zeros((4, 6)) - mat_b[0:3, 0:3] = rotation_matrix_2d(phi) - mat_b[0:3, 3:6] = rotation_matrix_2d(phi - two_theta) * (-1) - mat_b[3, 0] = 2 * ksq2eng * ki - mat_b[3, 3] = -2 * ksq2eng * kf - - # matrix C, constrinat between mono/ana mosaic and collimator divergence - mat_c = np.zeros(((CN.NUM_MONOS + CN.NUM_ANAS) * 2, CN.NUM_COLLS * 2)) - mat_c[CN.IDX_MONO0_H, CN.IDX_COLL0_H] = 0.5 - mat_c[CN.IDX_MONO0_H, CN.IDX_COLL1_H] = 0.5 - mat_c[CN.IDX_MONO0_V, CN.IDX_COLL0_V] = 0.5 / np.sin(theta_m) - mat_c[CN.IDX_MONO0_V, CN.IDX_COLL1_V] = -0.5 / np.sin(theta_m) - - mat_c[CN.IDX_ANA0_H, CN.IDX_COLL2_H] = 0.5 - mat_c[CN.IDX_ANA0_H, CN.IDX_COLL3_H] = 0.5 - mat_c[CN.IDX_ANA0_V, CN.IDX_COLL2_V] = 0.5 / np.sin(theta_a) - mat_c[CN.IDX_ANA0_V, CN.IDX_COLL3_V] = -0.5 / np.sin(theta_a) - - mat_h = mat_c.T @ self._mat_f @ mat_c - mat_hg_inv = la.inv(mat_h + self._mat_g) - mat_ba = mat_b @ mat_a - mat_cov = mat_ba @ mat_hg_inv @ mat_ba.T - - # TODO smaple mosaic???? - mat_cov[1, 1] += q**2 * self.sample.mosaic**2 - mat_cov[2, 2] += q**2 * self.sample.mosaic_v**2 - - mat_reso = la.inv(mat_cov) * sig2fwhm**2 - - # TODO normalization factor - if R0: # calculate - r0 = 1 - else: - r0 = 0 - - rez = ResoEllipsoid(mat_reso, r0) - - return rez diff --git a/src/tavi/instrument/resolution/backups/reso_ellipses_bak.py b/src/tavi/instrument/resolution/backups/reso_ellipses_bak.py deleted file mode 100755 index 5a4423f8..00000000 --- a/src/tavi/instrument/resolution/backups/reso_ellipses_bak.py +++ /dev/null @@ -1,374 +0,0 @@ -import numpy as np -import numpy.linalg as la -import mpl_toolkits.mplot3d as mplot3d -import matplotlib -import matplotlib.pyplot as plt -from tavi.utilities import * - - -np.set_printoptions(floatmode="fixed", precision=4) - - -class ResoEllipsoid(object): - """Manage the resolution ellipoid - - Attributs: - frame (str): "q", "hkl", or "proj" - STATUS (None | bool): True if resolution calculation is successful - mat (float): 4 by 4 resolution matrix - r0 (float | None): Normalization factor - - Methods: - - """ - - g_eps = 1e-8 - - def __init__(self): - - self.q = None - self.en = None - self.frame = None - self.projection = None - - self.STATUS = None - self.mat = None - self.r0 = None - self.coh_fwhms = None - self.incoh_fwhms = None - self.principal_fwhms = None - - # def __init__(self, mat, r0): - - # self.mat = mat - # self.r0 = r0 - # self.coh_fwhms = None - # self.incoh_fwhms = None - # self.principal_fwhms = None - - # if np.isnan(r0) or np.isinf(r0) or np.isnan(mat.any()) or np.isinf(mat.any()): - # self.STATUS = False - # else: - # self.STATUS = True - - def ellipsoid_volume(self): - """volume of the ellipsoid""" - pass - - @staticmethod - def quadric_proj(quadric, idx): - """projects along one axis of the quadric""" - - if np.abs(quadric[idx, idx]) < ResoEllipsoid.g_eps: - return np.delete(np.delete(quadric, idx, axis=0), idx, axis=1) - - # row/column along which to perform the orthogonal projection - vec = 0.5 * (quadric[idx, :] + quadric[:, idx]) # symmetrise if not symmetric - vec /= np.sqrt(quadric[idx, idx]) # normalise to indexed component - proj_op = np.outer(vec, vec) # projection operator - ortho_proj = quadric - proj_op # projected quadric - - # comparing with simple projection - # rank = len(quadric) - # vec /= np.sqrt(np.dot(vec, vec)) - # proj_op = np.outer(vec, vec) - # ortho_proj = np.dot((np.identity(rank) - proj_op), quadric) - - # remove row/column that was projected out - # 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 coh_fwhms_calc(self): - """Coherent FWHMs""" - reso = self.mat - fwhms = [] - - for i in range(len(reso)): - fwhms.append(sig2fwhm / np.sqrt(reso[i, i])) - fwhms = np.array(fwhms) - self.coh_fwhms = fwhms - - def incoh_fwhms_calc(self): - """Incoherent FWHMs""" - - reso = self.mat - proj_q_para = ResoEllipsoid.quadric_proj(reso, 3) - proj_q_para = ResoEllipsoid.quadric_proj(proj_q_para, 2) - proj_q_para = ResoEllipsoid.quadric_proj(proj_q_para, 1) - - proj_q_perp = ResoEllipsoid.quadric_proj(reso, 3) - proj_q_perp = ResoEllipsoid.quadric_proj(proj_q_perp, 2) - proj_q_perp = ResoEllipsoid.quadric_proj(proj_q_perp, 0) - - proj_q_up = ResoEllipsoid.quadric_proj(reso, 3) - proj_q_up = ResoEllipsoid.quadric_proj(proj_q_up, 1) - proj_q_up = ResoEllipsoid.quadric_proj(proj_q_up, 0) - - proj_en = ResoEllipsoid.quadric_proj(reso, 2) - proj_en = ResoEllipsoid.quadric_proj(proj_en, 1) - proj_en = ResoEllipsoid.quadric_proj(proj_en, 0) - - fwhms = np.array( - [ - 1.0 / np.sqrt(np.abs(proj_q_para[0, 0])) * sig2fwhm, - 1.0 / np.sqrt(np.abs(proj_q_perp[0, 0])) * sig2fwhm, - 1.0 / np.sqrt(np.abs(proj_q_up[0, 0])) * sig2fwhm, - 1.0 / np.sqrt(np.abs(proj_en[0, 0])) * sig2fwhm, - ] - ) - self.incoh_fwhms = fwhms - - def principal_fwhms_calc(self): - """FWHMs of principal axes in 4D""" - - evals, _ = la.eig(self.mat) - fwhms = np.array(1.0 / np.sqrt(np.abs(evals)) * sig2fwhm) - self.principal_fwhms = fwhms - - @staticmethod - def descr_ellipse(quadric): - """2D ellipse""" - - evals, evecs = la.eig(quadric) # evecs normalized already - fwhms = 1.0 / np.sqrt(np.abs(evals)) * sig2fwhm - # angles = np.array([np.arctan2(evecs[1][0], evecs[0][0])]) - - return (fwhms, evecs) - - @staticmethod - def _gen_ellipse_cut(mat, axes=(0, 1)): - """generate a 2D ellipsis by removing other axes""" - q_res = mat - 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) - - fwhms, vec = ResoEllipsoid.descr_ellipse(q_res) - return (fwhms, vec) - - @staticmethod - def _gen_ellipse_project(mat, axes=(0, 1)): - """generate a 2D ellipsis by projecting out other axes""" - q_res = mat - # 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) - - fwhms, rot = ResoEllipsoid.descr_ellipse(q_res) - return (fwhms, rot) - - def calc_ellipses(self, verbose=True): - """Calculate FWHMs and rotation angles""" - self.coh_fwhms_calc() - self.incoh_fwhms_calc() - self.principal_fwhms_calc() - - if verbose: - print(f"Coherent-elastic fwhms:{self.coh_fwhms}") - print(f"Incoherent-elastic fwhms: {self.incoh_fwhms}") - print(f"Principal axes fwhms: {self.principal_fwhms}") - - # 2d sliced ellipses - - fwhms_QxE, rot_QxE = ResoEllipsoid._gen_ellipse_cut(self.mat, axes=(0, 3)) - fwhms_QyE, rot_QyE = ResoEllipsoid._gen_ellipse_cut(self.mat, axes=(1, 3)) - fwhms_QzE, rot_QzE = ResoEllipsoid._gen_ellipse_cut(self.mat, axes=(2, 3)) - - fwhms_QxQy, rot_QxQy = ResoEllipsoid._gen_ellipse_cut(self.mat, axes=(0, 1)) - fwhms_QyQz, rot_QyQz = ResoEllipsoid._gen_ellipse_cut(self.mat, axes=(1, 2)) - fwhms_QxQz, rot_QxQz = ResoEllipsoid._gen_ellipse_cut(self.mat, axes=(0, 2)) - - if verbose: - print(f"2d Qx,E slice fwhms and vectors: {fwhms_QxE}, {rot_QxE}") - print(f"2d Qy,E slice fwhms and vectors:{fwhms_QyE},{rot_QyE}") - print(f"2d Qz,E slice fwhms and vectors:{fwhms_QzE},{rot_QzE}") - print(f"2d Qx,Qy slice fwhms and vectors:{fwhms_QxQy},{rot_QxQy}") - print(f"2d Qy,Qz slice fwhms and vectors:{fwhms_QyQz},{rot_QyQz}") - print(f"2d Qx,Qz slice fwhms and vectors:{fwhms_QxQz},{rot_QxQz}") - - # 2d projected ellipses - # Qres_QxE_proj = np.delete(np.delete(self.mat, 2, axis=0), 2, axis=1) - - (fwhms_QxE_proj, rot_QxE_proj) = ResoEllipsoid._gen_ellipse_project(self.mat, axes=(0, 3)) - (fwhms_QyE_proj, rot_QyE_proj) = ResoEllipsoid._gen_ellipse_project(self.mat, axes=(1, 3)) - (fwhms_QzE_proj, rot_QzE_proj) = ResoEllipsoid._gen_ellipse_project(self.mat, axes=(2, 3)) - (fwhms_QxQy_proj, rot_QxQy_proj) = ResoEllipsoid._gen_ellipse_project(self.mat, axes=(0, 1)) - (fwhms_QyQz_proj, rot_QyQz_proj) = ResoEllipsoid._gen_ellipse_project(self.mat, axes=(1, 2)) - (fwhms_QxQz_proj, rot_QxQz_proj) = ResoEllipsoid._gen_ellipse_project(self.mat, axes=(0, 2)) - - if verbose: - print(f"2d Qx,E projection fwhms and vectors: {fwhms_QxE_proj}, {rot_QxE_proj}") - print(f"2d Qy,E projection fwhms and vectors: {fwhms_QyE_proj}, {rot_QyE_proj}") - print(f"2d Qz,E projection fwhms and vectors: {fwhms_QzE_proj}, {rot_QzE_proj}") - print(f"2d Qx,Qy projection fwhms and vectors: {fwhms_QxQy_proj}, {rot_QxQy_proj}") - print(f"2d Qy,Qz projection fwhms and vectors: {fwhms_QyQz_proj}, {rot_QyQz_proj}") - print(f"2d Qx,Qz projection fwhms and vectors: {fwhms_QxQz_proj}, {rot_QxQz_proj}") - - results = { - "fwhms_QxE": fwhms_QxE, - "rot_QxE": rot_QxE, - "fwhms_QyE": fwhms_QyE, - "rot_QyE": rot_QyE, - "fwhms_QzE": fwhms_QzE, - "rot_QzE": rot_QzE, - "fwhms_QxQy": fwhms_QxQy, - "rot_QxQy": rot_QxQy, - "fwhms_QyQz": fwhms_QyQz, - "rot_QyQz": rot_QyQz, - "fwhms_QxQz": fwhms_QxQz, - "rot_QxQz": rot_QxQz, - "fwhms_QxE_proj": fwhms_QxE_proj, - "rot_QxE_proj": rot_QxE_proj, - "fwhms_QyE_proj": fwhms_QyE_proj, - "rot_QyE_proj": rot_QyE_proj, - "fwhms_QzE_proj": fwhms_QzE_proj, - "rot_QzE_proj": rot_QzE_proj, - "fwhms_QxQy_proj": fwhms_QxQy_proj, - "rot_QxQy_proj": rot_QxQy_proj, - "fwhms_QyQz_proj": fwhms_QyQz_proj, - "rot_QyQz_proj": rot_QyQz_proj, - "fwhms_QxQz_proj": fwhms_QxQz_proj, - "rot_QxQz_proj": rot_QxQz_proj, - } - - return results - - def plot_ellipses( - self, - ellis, - verbose=True, - plot_results=True, - file="", - dpi=600, - ellipse_points=128, - use_tex=False, - ): - """Plot 2D and 3D ellipses""" - matplotlib.rc("text", usetex=use_tex) - - def ellfkt(rad, vecs, phi): - return np.dot( - vecs, - np.array( - [ - rad[0] * np.cos(phi), - rad[1] * np.sin(phi), - ], - ), - ) - - phi = np.linspace(0, 2.0 * np.pi, ellipse_points) - - ell_QxE = ellfkt(ellis["fwhms_QxE"] * 0.5, ellis["rot_QxE"], phi) - ell_QyE = ellfkt(ellis["fwhms_QyE"] * 0.5, ellis["rot_QyE"], phi) - ell_QzE = ellfkt(ellis["fwhms_QzE"] * 0.5, ellis["rot_QzE"], phi) - ell_QxQy = ellfkt(ellis["fwhms_QxQy"] * 0.5, ellis["rot_QxQy"], phi) - ell_QyQz = ellfkt(ellis["fwhms_QyQz"] * 0.5, ellis["rot_QyQz"], phi) - ell_QxQz = ellfkt(ellis["fwhms_QxQz"] * 0.5, ellis["rot_QxQz"], phi) - - ell_QxE_proj = ellfkt(ellis["fwhms_QxE_proj"] * 0.5, ellis["rot_QxE_proj"], phi) - ell_QyE_proj = ellfkt(ellis["fwhms_QyE_proj"] * 0.5, ellis["rot_QyE_proj"], phi) - ell_QzE_proj = ellfkt(ellis["fwhms_QzE_proj"] * 0.5, ellis["rot_QzE_proj"], phi) - ell_QxQy_proj = ellfkt(ellis["fwhms_QxQy_proj"] * 0.5, ellis["rot_QxQy_proj"], phi) - ell_QyQz_proj = ellfkt(ellis["fwhms_QyQz_proj"] * 0.5, ellis["rot_QyQz_proj"], phi) - ell_QxQz_proj = ellfkt(ellis["fwhms_QxQz_proj"] * 0.5, ellis["rot_QxQz_proj"], phi) - - labelQpara = "Qpara (1/A)" - labelQperp = "Qperp (1/A)" - labelQup = "Qup (1/A)" - - if use_tex: - labelQpara = "$Q_{\parallel}$ (\AA$^{-1}$)" - labelQperp = "$Q_{\perp}$ (\AA$^{-1}$)" - labelQup = "$Q_{up}$ (\AA$^{-1}$)" - - # Qpara, E axis - - fig = plt.figure() - subplot_QxE = fig.add_subplot(231) - subplot_QxE.set_xlabel(labelQpara) - subplot_QxE.set_ylabel("E (meV)") - subplot_QxE.plot(ell_QxE[0], ell_QxE[1], c="black", linestyle="dashed") - subplot_QxE.plot(ell_QxE_proj[0], ell_QxE_proj[1], c="black", linestyle="solid") - subplot_QxE.grid(alpha=0.6) - - # Qperp, E axis - subplot_QyE = fig.add_subplot(232) - subplot_QyE.set_xlabel(labelQperp) - subplot_QyE.set_ylabel("E (meV)") - subplot_QyE.plot(ell_QyE[0], ell_QyE[1], c="black", linestyle="dashed") - subplot_QyE.plot(ell_QyE_proj[0], ell_QyE_proj[1], c="black", linestyle="solid") - subplot_QyE.grid(alpha=0.6) - # Qup, E axis - subplot_QzE = fig.add_subplot(233) - subplot_QzE.set_xlabel(labelQup) - subplot_QzE.set_ylabel("E (meV)") - subplot_QzE.plot(ell_QzE[0], ell_QzE[1], c="black", linestyle="dashed") - subplot_QzE.plot(ell_QzE_proj[0], ell_QzE_proj[1], c="black", linestyle="solid") - subplot_QzE.grid(alpha=0.6) - # Qpara, Qperp axis - subplot_QxQy = fig.add_subplot(234) - subplot_QxQy.set_xlabel(labelQpara) - subplot_QxQy.set_ylabel(labelQperp) - subplot_QxQy.plot(ell_QxQy[0], ell_QxQy[1], c="black", linestyle="dashed") - subplot_QxQy.plot(ell_QxQy_proj[0], ell_QxQy_proj[1], c="black", linestyle="solid") - subplot_QxQy.grid(alpha=0.6) - subplot_QxQy.set_aspect("equal") - # Qperp, Qup axis - subplot_QyQz = fig.add_subplot(235) - subplot_QyQz.set_xlabel(labelQperp) - subplot_QyQz.set_ylabel(labelQup) - subplot_QyQz.plot(ell_QyQz[0], ell_QyQz[1], c="black", linestyle="dashed") - subplot_QyQz.plot(ell_QyQz_proj[0], ell_QyQz_proj[1], c="black", linestyle="solid") - subplot_QyQz.grid(alpha=0.6) - subplot_QyQz.set_aspect("equal") - # Qpara, Qup axis - subplot_QxQz = fig.add_subplot(236) - subplot_QxQz.set_xlabel(labelQpara) - subplot_QxQz.set_ylabel(labelQup) - subplot_QxQz.plot(ell_QxQz[0], ell_QxQz[1], c="black", linestyle="dashed") - subplot_QxQz.plot(ell_QxQz_proj[0], ell_QxQz_proj[1], c="black", linestyle="solid") - subplot_QxQz.grid(alpha=0.6) - subplot_QxQz.set_aspect("equal") - - plt.tight_layout() - - # 3d plot - fig3d = plt.figure() - subplot3d = fig3d.add_subplot(111, projection="3d") - - subplot3d.set_xlabel(labelQpara) - subplot3d.set_ylabel(labelQperp) - # subplot3d.set_ylabel(labelQup) - subplot3d.set_zlabel("E (meV)") - - # xE - subplot3d.plot(ell_QxE[0], ell_QxE[1], zs=0.0, zdir="y", c="black", linestyle="dashed") - subplot3d.plot(ell_QxE_proj[0], ell_QxE_proj[1], zs=0.0, zdir="y", c="black", linestyle="solid") - # yE - subplot3d.plot(ell_QyE[0], ell_QyE[1], zs=0.0, zdir="x", c="black", linestyle="dashed") - subplot3d.plot(ell_QyE_proj[0], ell_QyE_proj[1], zs=0.0, zdir="x", c="black", linestyle="solid") - # zE - # subplot3d.plot(ell_QzE[0], ell_QzE[1], zs=0., zdir="x", c="black", linestyle="dashed") - # subplot3d.plot(ell_QzE_proj[0], ell_QzE_proj[1], zs=0., zdir="x", c="black", linestyle="solid") - # xy - subplot3d.plot(ell_QxQy[0], ell_QxQy[1], zs=0.0, zdir="z", c="black", linestyle="dashed") - subplot3d.plot(ell_QxQy_proj[0], ell_QxQy_proj[1], zs=0.0, zdir="z", c="black", linestyle="solid") - - # save plots to files - if file != "": - import os - - splitext = os.path.splitext(file) - file3d = splitext[0] + "_3d" + splitext[1] - - if verbose: - print('Saving 2d plot to "%s".' % file) - print('Saving 3d plot to "%s".' % file3d) - fig.savefig(file, dpi=dpi) - fig3d.savefig(file3d, dpi=dpi) - - # show plots - if plot_results: - plt.show() diff --git a/src/tavi/instrument/resolution/cooper_nathans.py b/src/tavi/instrument/resolution/cooper_nathans.py index f08a19cd..8dc33646 100755 --- a/src/tavi/instrument/resolution/cooper_nathans.py +++ b/src/tavi/instrument/resolution/cooper_nathans.py @@ -2,11 +2,11 @@ import numpy.linalg as la from tavi.instrument.resolution.reso_ellipses import ResoEllipsoid -from tavi.instrument.tas import TAS +from tavi.instrument.tas import TripleAxisSpectrometer from tavi.utilities import * -class CN(TAS): +class CN(TripleAxisSpectrometer): """Copper-Nathans method Attibutes: @@ -38,8 +38,9 @@ class CN(TAS): IDX_ANA0_H = 2 IDX_ANA0_V = 3 - def __init__(self): - super().__init__() + def __init__(self, path_to_json=None) -> None: + """Load instrument configuration from json if provided""" + super().__init__(path_to_json) # constants independent of q and eng self._mat_f = None @@ -47,11 +48,11 @@ def __init__(self): def cooper_nathans( self, - ei, - ef, - hkl, - projection=((1, 0, 0), (0, 1, 0), (0, 0, 1)), - R0=False, + ei: float, + ef: float, + hkl: tuple[float], + projection: tuple[tuple] = ((1, 0, 0), (0, 1, 0), (0, 0, 1)), + R0: bool = False, ): """Calculate resolution using Cooper-Nathans method @@ -232,9 +233,9 @@ def cooper_nathans( mat_cov = mat_ba @ mat_h_inv @ mat_ba.T # TODO smaple mosaic???? - sample_mosaic = np.deg2rad(self.sample.mosaic / 60) + sample_mosaic_h = np.deg2rad(self.sample.mosaic_h / 60) sample_mosaic_v = np.deg2rad(self.sample.mosaic_v / 60) - mat_cov[1, 1] += q_mod**2 * sample_mosaic**2 + mat_cov[1, 1] += q_mod**2 * sample_mosaic_h**2 mat_cov[2, 2] += q_mod**2 * sample_mosaic_v**2 mat_reso = la.inv(mat_cov) * sig2fwhm**2 diff --git a/src/tavi/instrument/tas.py b/src/tavi/instrument/tas.py index d469e2d0..8dd42489 100644 --- a/src/tavi/instrument/tas.py +++ b/src/tavi/instrument/tas.py @@ -1,63 +1,95 @@ import json +from pathlib import Path +from typing import Optional, Union import numpy as np -from tavi.instrument.tas_cmponents import * +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 Arms, Collimators, Detector, Guide, Monitor, Source from tavi.sample.powder import Powder from tavi.sample.sample import Sample from tavi.sample.xtal import Xtal from tavi.utilities import * -class TAS(object): - """Triple-axis instrument class. Manage instrument congigutarion parameters +class TripleAxisSpectrometer(object): + """ + Triple-axis instrument class. Manage instrument congigutarion parameters Attibutes: + source (Source): + collimators (Collimators): + guide (Guide): + monichromator (MonoAna): + goniometer (Goniometer): + analyzer (MonoAna): + detector (Detector): + arms (Arms): + sample (Sample | Xtal | Powder): Methods: - load_instrument(config_params): - load_sample(sample_params): - - + load_instrument_params_from_json """ - def __init__(self, path_to_json=None): - """load default instrument configuration""" - # TODO - self.source = None - self.collimators = None - self.guide = None - self.monochromator = None - self.monitor = None - self.goniometer = None - self.analyzer = None - self.detector = None - self.arms = None - self.sample = None + def __init__( + self, + path_to_json: str = None, + ): + """Load instrument configuration from json if provided""" + self.source: Optional[Source] = None + self.collimators: Optional[Collimators] = None + self.guide: Optional[Guide] = None + self.monochromator: Optional[MonoAna] = None + self.monitor: Optional[Monitor] = None + self.goniometer: Optional[Goniometer] = None + self.analyzer: Optional[MonoAna] = None + self.detector: Optional[Detector] = None + self.arms: Optional[Arms] = None + self.sample: Union[Sample, Xtal, Powder, None] = None if path_to_json is not None: - self.load_instrument_from_json(path_to_json) - - def load_instrument_from_dicts(self, config_params): + self.load_instrument_params_from_json(path_to_json) + + def _load_instrument_parameters(self, config_params: dict[dict]): + components = { + "source": Source, + "collimators": Collimators, + "guide": Guide, + "monochromator": MonoAna, + "monitor": Monitor, + "goniometer": Goniometer, + "analyzer": MonoAna, + "detector": Detector, + "distances": Arms, + } + + for component_name, component_class in components.items(): + param = config_params.get(component_name) + if param is not None: + setattr(self, component_name, component_class(param, component_name)) + + def load_instrument_params_from_json( + self, + path_to_json: str, + ): """Load instrument configuration from a json file""" - self.source = Source(config_params["source"]) - self.collimators = Collimators(config_params["collimators"]) - self.guide = Guide(config_params["guide"]) - self.monochromator = Monochromator(config_params["monochromator"]) - self.monitor = Monitor(config_params["monitor"]) - self.goniometer = Goniometer(config_params["goniometer"]) - self.analyzer = Analyzer(config_params["analyzer"]) - self.detector = Detector(config_params["detector"]) - self.arms = Arms(config_params["distances"]) - - def load_instrument_from_json(self, path_to_json): - """Load instrument configuration from a json file""" + json_file = Path(path_to_json) + if json_file.is_file(): + with open(json_file, "r", encoding="utf-8") as file: + config_params = json.load(file) - with open(path_to_json, "r", encoding="utf-8") as file: - config_params = json.load(file) + self._load_instrument_parameters(config_params) + else: + print("Invalid path for instrument configuration json file.") - self.load_instrument_from_dicts(config_params) + def load_insrument_params_from_scan( + self, + scan: Scan, + ): + pass # def save_instrument(self): # """Save configuration into a dictionary""" @@ -78,11 +110,11 @@ def load_sample_from_json(self, path_to_json): try: # is type xtal ot powder? sample_type = sample_params["type"] if sample_type == "xtal": - sample = Xtal.from_json(sample_params) + sample = Xtal.from_json(path_to_json) elif sample_type == "powder": - sample = Powder.from_json(sample_params) + sample = Powder.from_json(path_to_json) except KeyError: # sample type is not given - sample = Sample.from_json(sample_params) + sample = Sample.from_json(path_to_json) self.load_sample(sample) @@ -139,8 +171,8 @@ def _find_u_from_2peaks(self, peaks, angles, ki, kf): two_theta1, _, _, _ = angles1 two_theta2, _, _, _ = angles2 - q_lab1 = TAS.q_lab(two_theta1, ki=ki, kf=kf) - q_lab2 = TAS.q_lab(two_theta2, ki=ki, kf=kf) + q_lab1 = TripleAxisSpectrometer.q_lab(two_theta1, ki=ki, kf=kf) + q_lab2 = TripleAxisSpectrometer.q_lab(two_theta2, ki=ki, kf=kf) # Goniometer angles all zeros in q_sample frame q_sample1 = self.goniometer.r_mat_inv(angles1[1:]) @ q_lab1 @@ -309,7 +341,7 @@ def find_angles(self, peak, ei=13.5, ef=None): t_mat_inv = np.linalg.inv(t_mat) - q_lab1 = TAS.q_lab(two_theta, ki, kf) / q_norm + q_lab1 = TripleAxisSpectrometer.q_lab(two_theta, ki, kf) / q_norm q_lab2 = np.array([q_lab1[2], 0, -q_lab1[0]]) q_lab3 = np.array([0, 1, 0]) diff --git a/src/tavi/instrument/tas_cmponents.py b/src/tavi/instrument/tas_cmponents.py index 8f954092..f7b69ca6 100644 --- a/src/tavi/instrument/tas_cmponents.py +++ b/src/tavi/instrument/tas_cmponents.py @@ -1,123 +1,197 @@ +from typing import Literal, Optional, Union + import numpy as np -from tavi.utilities import * +from tavi.utilities import cm2angstrom -class Source(object): - """Neutron source""" +class TASComponent(object): + """ + A component of the triple-axis spectrometer - def __init__(self, param_dict): - self.name = None - self.shape = "rectangular" - self.width = None - self.height = None + Attributes: + type (str): identify which component + Methods: + update(param, value): update a parameter + """ - for key, val in param_dict.items(): - match key: - case "width" | "height": - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - if param_dict["shape"] == "rectangular": - setattr(self, key, val / np.sqrt(12) * cm2angstrom) - elif param_dict["shape"] == "circular": - setattr(self, key, val / 4 * cm2angstrom) - else: - print("Source shape needs to be either rectangular or circular.") - case _: - setattr(self, key, val) + 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", + ) -> 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( + length_in_cm: Optional[float], + shape: str, + param_name: str = "Length", + ) -> float: + """ + Convert length from centimeter to angstrom. -class Guide(object): - """Guide""" + 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 - def __init__(self, param_dict): - self.in_use = False - self.div_h = None - self.div_v = None - for key, val in param_dict.items(): - match key: - case "div_h" | "div_v": - setattr(self, key, val * min2rad) - case _: - setattr(self, key, val) +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[str] = None + self.height: Optional[str] = None -class Monochromator(object): - """Monochromator""" + super().__init__(param_dict, component_name) + pass - def __init__(self, param_dict): - self.type = "PG002" - self.d_spacing = mono_ana_xtal["PG002"] - self.mosaic = 45 # horizontal mosaic, min of arc - self.mosaic_v = 45 # vertical mosaic, if anisotropic - self.sense = -1 + @property + def _width(self): + """width in angstrom, with correction based on shape, for resolution calculation""" + return TASComponent._cm2angstrom(self.width, self.shape, "Source width") - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - self.shape = "rectangular" # or spherical - self.width = 12.0 - self.height = 8.0 - self.depth = 0.15 - # horizontal focusing - self.curved_h = False - self.curvh = 0.0 - self.optimally_curved_h = False - # vertical focusing - self.curved_v = False - self.curvv = 0.0 - self.optimally_curved_v = False + @property + def _height(self): + """Height in angstrom, with correction based on shape, for resolution calculation""" + return TASComponent._cm2angstrom(self.height, self.shape, "Source height") - for key, val in param_dict.items(): - match key: - # case "mosaic" | "mosaic_v" | "curvh" | "curvv": - # setattr(self, key, val * min2rad) - case "width" | "height" | "depth": - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - if param_dict["shape"] == "rectangular": - setattr(self, key, val / np.sqrt(12) * cm2angstrom) - elif param_dict["shape"] == "spherical": - setattr(self, key, val / 4 * cm2angstrom) - else: - print("Monochromator shape needs to be either rectangular or spherical.") - case _: - setattr(self, key, val) - self.d_spacing = mono_ana_xtal[self.type] +class Guide(TASComponent): + """Neutron guide in the monochromator drum""" -class Collimators(object): - """collimitor divergences, in mins of arc""" + def __init__( + self, + param_dict: Optional[dict] = None, + component_name: str = "guide", + ): + self.in_use: bool = False + self.div_h: Optional[str] = None + self.div_v: Optional[str] = None - def __init__(self, param_dict): - self.h_pre_mono = 30 # mins of arc - self.h_pre_sample = 30 - self.h_post_sample = 30 - self.h_post_ana = 30 - self.v_pre_mono = 30 - self.v_pre_sample = 30 - self.v_post_sample = 30 - self.v_post_ana = 30 + super().__init__(param_dict, component_name) - for key, val in param_dict.items(): - # if key in ( - # "h_pre_mono", - # "h_pre_sample", - # "h_post_sample", - # "h_post_ana", - # "v_pre_mono", - # "v_pre_sample", - # "v_post_sample", - # "v_post_ana", - # ): - # setattr(self, key, val * min2rad) - # else: - setattr(self, key, val) + @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") -# TODO + +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 +# TODO ========================================================= class Monitor(object): - def __init__(self, param_dict): + def __init__(self, param_dict, component_name): self.shape = "rectangular" # or spherical self.width = 5 * cm2angstrom self.height = 12 * cm2angstrom @@ -137,51 +211,9 @@ def __init__(self, param_dict): setattr(self, key, val) -class Analyzer(object): - def __init__(self, param_dict): - self.type = "Pg002" - self.d_spacing = mono_ana_xtal["Pg002"] - self.mosaic = 45 # horizontal mosaic, min of arc - self.mosaic_v = 45 # vertical mosaic, if anisotropic - self.sense = -1 - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - self.shape = "rectangular" - self.width = 12.0 - self.height = 8.0 - self.depth = 0.3 - # horizontal focusing - self.curved_h = False - self.curvh = 0.0 - self.optimally_curved_h = False - # vertical focusing - self.curved_v = False - self.curvv = 0.0 - self.optimally_curved_v = False - - for key, val in param_dict.items(): - match key: - # case "mosaic" | "mosaic_v" | "curvh" | "curvv": - # setattr(self, key, val * min2rad) - case "width" | "height" | "depth": - # divide by np.sqrt(12) if rectangular - # Diameter D/4 if spherical - if param_dict["shape"] == "rectangular": - setattr(self, key, val / np.sqrt(12) * cm2angstrom) - elif param_dict["shape"] == "spherical": - setattr(self, key, val / 4 * cm2angstrom) - else: - print("Analyzer shape needs to be either rectangular or spherical.") - - setattr(self, key, val * cm2angstrom) - case _: - setattr(self, key, val) - self.d_spacing = mono_ana_xtal[self.type] - - # TODO class Detector(object): - def __init__(self, param_dict): + def __init__(self, param_dict, component_name): self.shape = "rectangular" self.width = 1.5 # in cm self.height = 5.0 # in cm @@ -204,7 +236,7 @@ def __init__(self, param_dict): class Arms(object): """lengths of arms, in units of cm""" - def __init__(self, param_dict): + def __init__(self, param_dict, component_name): # in units of cm self.src_mono = 0 self.mono_sample = 0 @@ -214,82 +246,3 @@ def __init__(self, param_dict): for key, val in param_dict.items(): setattr(self, key, val * cm2angstrom) - - -class Goniometer(object): - """Goniometer table, type = Y-ZX or YZ-X""" - - def __init__(self, param_dict): - self.type = "Y-ZX" # Y-mZ-X for Huber stage at HB1A and HB3 - self.sense = -1 - - for key, val in param_dict.items(): - setattr(self, key, val) - - def r_mat(self, angles): - "rotation matrix" - - omega, sgl, sgu = angles # s2, s1, sgl, sgu - match self.type: - case "Y-ZX": # HB3 - r_mat = rot_y(omega) @ rot_z(-1 * sgl) @ rot_x(sgu) - case "YZ-X": # CG4C - r_mat = rot_y(omega) @ rot_z(sgl) @ rot_x(-1 * sgu) - case _: - r_mat = None - print("Unknow goniometer type. Curruntly support Y-ZX and YZ-X") - - return r_mat - - def r_mat_inv( - self, - angles, - ): - """inverse of rotation matrix""" - # return np.linalg.inv(self.r_mat(angles)) - return self.r_mat(angles).T - - def angles_from_r_mat(self, r_mat): - """Calculate goniometer angles from the R matrix - - Note: - range of np.arcsin is -pi/2 to pi/2 - range of np.atan2 is -pi to pi - """ - - match self.type: - case "Y-ZX" | "YZ-X": # Y-mZ-X (s1, sgl, sgu) for HB1A and HB3, Y-Z-mX (s1, sgl, sgu) for CG4C - # sgl1 = np.arcsin(r_mat[1, 0]) * rad2deg - # sgl2 = np.arccos(np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg - sgl = np.arctan2(r_mat[1, 0], np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg - - # sgu1 = np.arcsin(-r_mat[1, 2] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg - # sgu2 = np.arccos(r_mat[1, 1] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg - sgu = ( - np.arctan2( - -r_mat[1, 2] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2), - r_mat[1, 1] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2), - ) - * rad2deg - ) - - # omega1 = np.arcsin(-r_mat[2, 0] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg - # omega2 = np.arccos(r_mat[0, 0] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2)) * rad2deg - omega = ( - np.arctan2( - -r_mat[2, 0] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2), - r_mat[0, 0] / np.sqrt(r_mat[0, 0] ** 2 + r_mat[2, 0] ** 2), - ) - * rad2deg - ) - match self.type: - case "Y-ZX": - angles = (omega, -1 * sgl, sgu) - case "YZ-X": - angles = (omega, sgl, -1 * sgu) - - case _: - angles = None - print("Unknow goniometer type. Curruntly support Y-ZX and YZ-X.") - - return angles diff --git a/src/tavi/sample/powder.py b/src/tavi/sample/powder.py index 268d3d7c..7db9e35f 100755 --- a/src/tavi/sample/powder.py +++ b/src/tavi/sample/powder.py @@ -13,19 +13,3 @@ class Powder(Sample): def __init__(self, lattice_params): super().__init__(lattice_params) self.type = "powder" - - @classmethod - def from_json(cls, sample_params): - """Alternate constructor from json""" - lattice_params = ( - sample_params["a"], - sample_params["b"], - sample_params["c"], - sample_params["alpha"], - sample_params["beta"], - sample_params["gamma"], - ) - - sample = cls(lattice_params=lattice_params) - - return sample diff --git a/src/tavi/sample/sample.py b/src/tavi/sample/sample.py index 2ca2644b..4af36e58 100755 --- a/src/tavi/sample/sample.py +++ b/src/tavi/sample/sample.py @@ -1,18 +1,20 @@ # -*- coding: utf-8 -*- -import numpy as np +import json -# from tavi.utilities import * +import numpy as np class Sample(object): """ + Sample class + Attributes: shape (str): "cuboid" or "cylindrical" width (float): in units of cm height (float): in units of cm depth (float): in units of cm - mosaic (fload): in units of minutes of arc + mosaic_h (fload): in units of minutes of arc mosaic_v (fload): verital mosaic if anisotropic, in units of minutes of arc a, b, c lattice constants in Angstrom @@ -23,9 +25,6 @@ class Sample(object): a_star_vec, b_star_vec, c_star_vec reciprocal lattice vector i_star, j_star, k_star bases for the reciprocal space lattice vectors - - - Methods: real_vec_cart reciprocal_vec_cart @@ -33,10 +32,12 @@ class Sample(object): reciprocal_basis b_mat - Static Methods: v_alpha_beta_gamma_calc(alpha,m beta, gamma) + Class Methods: + from_json(path_to_json): construct sample from a json file + """ @@ -55,9 +56,14 @@ def __init__( self.set_mosaic() # with defalt values self.set_shape() # with defalt values + # TODO @classmethod - def from_json(cls, sample_params): + def from_json(cls, path_to_json): """Alternate constructor from json""" + + with open(path_to_json, "r", encoding="utf-8") as file: + sample_params = json.load(file) + lattice_params = ( sample_params["a"], sample_params["b"], @@ -66,7 +72,29 @@ def from_json(cls, sample_params): sample_params["beta"], sample_params["gamma"], ) - return cls(lattice_params=lattice_params) + sample = cls(lattice_params=lattice_params) + + ub_matrix = sample_params.get("ub_matrix") + if ub_matrix is not None: + sample.ub_matrix = np.array(ub_matrix).reshape(3, 3) + + plane_normal = sample_params.get("plane_normal") + if plane_normal is not None: + sample.plane_normal = np.array(plane_normal) + + shape = sample_params.get("shape") + width = sample_params.get("width") + height = sample_params.get("height") + depth = sample_params.get("depth") + if all([shape, width, height, depth]): + sample.set_shape(shape, width, height, depth) + + mosaic_h = sample_params.get("mosaic_h") + mosaic_v = sample_params.get("mosaic_v") + if all([mosaic_h, mosaic_v]): + sample.set_mosaic(mosaic_h, mosaic_v) + + return sample def set_shape( self, @@ -83,11 +111,11 @@ def set_shape( def set_mosaic( self, - mosaic: float = 30, # horizontal mosaic + mosaic_h: float = 30, # horizontal mosaic mosaic_v: float = 30, # vertical mosaic ) -> None: """Set horizontal and vertical mosaic in units of minitues of arc""" - self.mosaic_h = mosaic # * min2rad + self.mosaic_h = mosaic_h # * min2rad self.mosaic_v = mosaic_v # * min2rad def update_lattice_parameters( @@ -124,6 +152,12 @@ def update_lattice_parameters( self.c_star_vec, ) = self._reciprocal_space_vectors() + ( + self.i_star, + self.j_star, + self.k_star, + ) = self._reciprocal_basis() + @staticmethod def v_alpha_beta_gamma_calc(alpha, beta, gamma) -> float: """ diff --git a/src/tavi/sample/xtal.py b/src/tavi/sample/xtal.py index 7b4658cc..9de3c801 100755 --- a/src/tavi/sample/xtal.py +++ b/src/tavi/sample/xtal.py @@ -1,25 +1,31 @@ import numpy as np from tavi.sample.sample import Sample -from tavi.utilities import * class Xtal(Sample): - """Singel crystal sample + """ + Singel crystal class Attibutes: type (str): "xtal" - inv_ub_matrix + ub_peaks (list(tuple)): peaks used to determine UB matrix + ub_angles (list(tuple)): goniometer angles for the peaks used to determine UB matrix + ub_matrix (np.adarray): UB matrix + inv_ub_matrix (np.ndarray): inverse of UB matrix in_plane_ref: in plane vector in Qsample frame, goniometers at zero plane_normal: normal vector in Qsample frame, goniometers at zero - u (tuple) - v (tuple) + u (tuple): u vector, (h, k, l) along the beam direction when all goniometer angles are zero + v (tuple): v vector, (h ,k, l) in the scattering plane Methods: """ - def __init__(self, lattice_params=(1, 1, 1, 90, 90, 90)): + def __init__( + self, + lattice_params=(1, 1, 1, 90, 90, 90), + ) -> None: super().__init__(lattice_params) self.type = "xtal" @@ -31,46 +37,13 @@ def __init__(self, lattice_params=(1, 1, 1, 90, 90, 90)): self.plane_normal = None self.in_plane_ref = None - self.i_star, self.j_star, self.k_star = self._reciprocal_basis() - - @classmethod - def from_json(cls, sample_params): - """Alternate constructor from json""" - lattice_params = ( - sample_params["a"], - sample_params["b"], - sample_params["c"], - sample_params["alpha"], - sample_params["beta"], - sample_params["gamma"], - ) - - sample = cls(lattice_params=lattice_params) - - sample.ub_matrix = np.array(sample_params["ub_matrix"]).reshape(3, 3) - sample.plane_normal = np.array(sample_params["plane_normal"]) - - param_dict = ("shape", "width", "height", "depth", "mosaic", "mosaic_v") - - for key, val in sample_params.items(): - match key: - case "height" | "width" | "depth": - setattr(sample, key, val * cm2angstrom) - # case "mosaic" | "mosaic_v": - # setattr(sample, key, val * min2rad) - case _: - if key in param_dict: - setattr(sample, key, val) - sample.update_lattice_parameters(lattice_params) - return sample - @property - def u(self): + def u(self) -> np.ndarray: """u vector, in reciprocal lattice unit, along beam""" return self.ub_matrix_to_uv(self.ub_matrix)[0] @property - def v(self): + def v(self) -> np.ndarray: """ v vector, in reciprocal lattice unit, in the horizaontal scattering plane @@ -78,7 +51,7 @@ def v(self): return self.ub_matrix_to_uv(self.ub_matrix)[1] @staticmethod - def ub_matrix_to_uv(ub_matrix): + def ub_matrix_to_uv(ub_matrix) -> tuple[np.ndarray]: """ "Calculate u and v vector from UB matrix""" inv_ub_matrix = np.linalg.inv(ub_matrix) u = inv_ub_matrix @ np.array([0, 0, 1]) @@ -92,12 +65,16 @@ def ub_matrix_to_lattice_params(ub_matrix): a = np.sqrt(g_mat[0, 0]) b = np.sqrt(g_mat[1, 1]) c = np.sqrt(g_mat[2, 2]) - alpha = np.arccos((g_mat[1, 2] + g_mat[2, 1]) / (2 * b * c)) * rad2deg - beta = np.arccos((g_mat[0, 2] + g_mat[2, 0]) / (2 * a * c)) * rad2deg - gamma = np.arccos((g_mat[0, 1] + g_mat[1, 0]) / (2 * a * b)) * rad2deg + alpha_rad = np.arccos((g_mat[1, 2] + g_mat[2, 1]) / (2 * b * c)) + beta_rad = np.arccos((g_mat[0, 2] + g_mat[2, 0]) / (2 * a * c)) + gamma_rad = np.arccos((g_mat[0, 1] + g_mat[1, 0]) / (2 * a * b)) + alpha = np.rad2deg(alpha_rad) + beta = np.rad2deg(beta_rad) + gamma = np.rad2deg(gamma_rad) + return (a, b, c, alpha, beta, gamma) - def uv_to_ub_matrix(self, u, v): + def uv_to_ub_matrix(self, u, v) -> np.ndarray: """Calculate UB matrix from u and v vector, and lattice parameters""" b_mat = self.b_mat() diff --git a/src/tavi/utilities.py b/src/tavi/utilities.py index e4fcc191..563319b0 100644 --- a/src/tavi/utilities.py +++ b/src/tavi/utilities.py @@ -13,24 +13,6 @@ min2rad = 1.0 / 60.0 / 180.0 * np.pi rad2deg = 180.0 / np.pi -# -------------------------------------------------------------------------- -# d_spacing table from Shirane Appendix 3, in units of Angstrom, from -# -------------------------------------------------------------------------- -mono_ana_xtal = { - "PG002": 3.35416, - "Pg002": 3.35416, - "PG004": 1.67708, - "Cu111": 2.08717, - "Cu220": 1.27813, - "Ge111": 3.26627, - "Ge220": 2.00018, - "Ge311": 1.70576, - "Ge331": 1.29789, - "Be002": 1.79160, - "Be110": 1.14280, - "Heusler": 3.435, # Cu2MnAl(111) -} - # -------------------------------------------------------------------------- # helper functions @@ -85,72 +67,3 @@ def rotation_matrix_2d(phi): ] ) return mat - - -def rot_x(nu): - """rotation matrix about y-axis by angle nu - - Args: - nu (float): angle in degrees - - Note: - Using Mantid convention, beam along z, y is up, x in plane - """ - - angle = nu / rad2deg - c = np.cos(angle) - s = np.sin(angle) - mat = np.array( - [ - [1, 0, 0], - [0, c, -s], - [0, s, c], - ] - ) - return mat - - -def rot_y(omega): - """rotation matrix about y-axis by angle omega - - Args: - omega (float): angle in degrees - - Note: - Using Mantid convention, beam along z, y is up, x in plane - """ - - angle = omega / rad2deg - c = np.cos(angle) - s = np.sin(angle) - mat = np.array( - [ - [c, 0, s], - [0, 1, 0], - [-s, 0, c], - ] - ) - return mat - - -def rot_z(mu): - """rotation matrix about z-axis by angle mu - - Args: - mu (float): angle in degrees - - Note: - Using Mantid convention, beam along z, y is up, x in plane - """ - - angle = mu / rad2deg - c = np.cos(angle) - s = np.sin(angle) - mat = np.array( - [ - [c, -s, 0], - [s, c, 0], - [0, 0, 1], - ] - ) - return mat diff --git a/test_data/nexus_exp1031.h5 b/test_data/nexus_exp1031.h5 index eb356380f29aa5ec7fc5be3b48aaae558afaa3d7..b23e38ec1b906afa18333ff6597dbf293c92732c 100644 GIT binary patch delta 401 zcmXZN*G^Lb0Dxf}1t*9D6c_G=M@|cr)`5GW;;4!f6u00WbyxXBg}TMMwO8I5FHC$4 z6Q6)b(TIQY=j;2`73uqX#sC8iGFU{RB0~%{%y1)&RIJ1(rA8ZLtZ~Ye8*hRN6HPMN z6jN22X1W<>nq{^*=9*`|1*$9*SfpAgDyBxQI&t+9k``NHsbv~0x57%RthUBl>#VoI zMvXSvY>OsaZL{4DDLd`5TeCg(+GoE54ranv!#op91Ar?u;F##!f_cfm!M zTz17(Y1dqL!%er`c1Olt_jKaEE)P8PNVgtYk3I3!Gta&7(ks1Q%X#Cici#Kpqfb8j z@*|%~2hlJU#*?uySrZ1eiFhJfU+}N0ehLcz8w(;*Vu2*J!0R z1^Or&{APah4E*c~4-B4`V~C-0g$*;@2qTR$+8BAp8fUzG1tyrN&?J*hG1WBF6`5hC zS&Ge8VvbT}%FQ*;d=(a`v`}D?D%BQ?h^nzfOk6@zt)-T!vs}FuR$67XH5xQ(vR1Pe zt=3tuO}h;?+GMjWw%TU99d_!l%Wivg+H0Tv4mjwL!;U!WnBz`3>6FvXIP0AAF1YBD z%dSYd>YD3rxapQIX}8^R7u|Z?bKe6G^?Ky7C!Tudxffn~<+V35-g@V~4?g_OV-4L>O?FWkA(hK`HxWezmXsjRE56xr?S6-O!4owaHjifw&M@{ Cxtg5- diff --git a/test_data/test_samples/sample.json b/test_data/test_samples/sample.json index 2b4f9a7b..5a0f0e0d 100644 --- a/test_data/test_samples/sample.json +++ b/test_data/test_samples/sample.json @@ -1,6 +1,6 @@ { - "a": 5.034785, - "b": 5.034785, + "a": 5.1, + "b": 5.1, "c": 13.812004, "alpha": 90, "beta": 90, diff --git a/tests/test_cooper_nathans.py b/tests/test_cooper_nathans.py index 15b4d78e..2734ab4e 100755 --- a/tests/test_cooper_nathans.py +++ b/tests/test_cooper_nathans.py @@ -1,5 +1,6 @@ import numpy as np import pytest + from tavi.instrument.resolution.cooper_nathans import CN np.set_printoptions(floatmode="fixed", precision=4) @@ -53,12 +54,11 @@ def test_copper_nathans_projection(tas_params): @pytest.fixture def tas_params(): # cooper_nathans_CTAX - tas = CN() instrument_config_json_path = "./src/tavi/instrument/instrument_params/cg4c.json" - sample_json_path = "./test_data/test_samples/nitio3.json" + tas = CN(instrument_config_json_path) - tas.load_instrument_from_json(instrument_config_json_path) + sample_json_path = "./test_data/test_samples/nitio3.json" tas.load_sample_from_json(sample_json_path) ei = 4.8 diff --git a/tests/test_sample.py b/tests/test_sample.py index d1391d69..2304f018 100644 --- a/tests/test_sample.py +++ b/tests/test_sample.py @@ -1,7 +1,7 @@ import numpy as np import pytest -from tavi.instrument.tas import TAS +from tavi.sample.sample import Sample from tavi.sample.xtal import Xtal np.set_printoptions(floatmode="fixed", precision=4) @@ -69,7 +69,15 @@ def test_uv_to_ub_matrix(xtal_info): assert np.allclose(ub_matrix_calc, ub_matrix, atol=1e-2) -def test_load_sample_from_json(): +def test_load_generic_sample_from_json(): sample_json_path = "./test_data/test_samples/sample.json" - tas = TAS() - tas.load_sample_from_json(sample_json_path) + sample = Sample.from_json(sample_json_path) + assert np.allclose(sample.a, 5.1) + + +def test_load_xtal_from_json(): + xtal_json_path = "./test_data/test_samples/nitio3.json" + xtal = Xtal.from_json(xtal_json_path) + assert xtal.type == "xtal" + assert np.allclose(xtal.a, 5.034785) + assert np.shape(xtal.ub_matrix) == (3, 3) diff --git a/tests/test_tas.py b/tests/test_tas.py index 855523e4..43e885e2 100644 --- a/tests/test_tas.py +++ b/tests/test_tas.py @@ -1,13 +1,14 @@ import numpy as np -from tavi.instrument.tas import TAS + +from tavi.instrument.tas import TripleAxisSpectrometer from tavi.sample.xtal import Xtal def test_load_from_json(): - tax = TAS() + tax = TripleAxisSpectrometer() cg4c = "./src/tavi/instrument/instrument_params/cg4c.json" nitio3 = "./test_data/test_samples/nitio3.json" - tax.load_instrument_from_json(cg4c) + tax.load_instrument_params_from_json(cg4c) tax.load_sample_from_json(nitio3) assert tax.analyzer.type == "Pg002" @@ -39,9 +40,9 @@ def test_calc_ub_from_2_peaks_hb3(): plane_normal = [0.000009, 0.999047, 0.043637] in_plane_ref = [0.942840, 0.014534, -0.332928] - tas = TAS() + tas = TripleAxisSpectrometer() takin = "./src/tavi/instrument/instrument_params/takin_test.json" - tas.load_instrument_from_json(takin) + tas.load_instrument_params_from_json(takin) tas.load_sample(Xtal(lattice_params)) peak_list = [ @@ -93,10 +94,10 @@ def test_calc_ub_from_2_peaks_ctax(): plane_normal = [-0.04032, 0.998565, -0.035237] in_plane_ref = [-0.993257, -0.043892, -0.107299] - ctax = TAS() + ctax = TripleAxisSpectrometer() cg4c = "./src/tavi/instrument/instrument_params/cg4c.json" nitio3 = "./test_data/test_samples/nitio3.json" - ctax.load_instrument_from_json(cg4c) + ctax.load_instrument_params_from_json(cg4c) ctax.load_sample_from_json(nitio3) peak_list = [ @@ -126,8 +127,8 @@ def test_calc_ub_from_2_peaks_ctax(): def test_calc_ub_from_2_peaks_hb1(): - hb1 = TAS() - hb1.load_instrument_from_json("./src/tavi/instrument/instrument_params/takin_test.json") + hb1 = TripleAxisSpectrometer() + hb1.load_instrument_params_from_json("./src/tavi/instrument/instrument_params/takin_test.json") lattice_params = (3.939520, 3.939520, 3.941957, 90.0, 90.0, 90.0)