From ea6f5bfd6e12b6d1fa6814ed6146fa69cebb1bb5 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 24 Aug 2022 18:01:23 +0200 Subject: [PATCH 01/14] [#673] Updating Si units, reshaping dictionnaries & made the computation routines more generic --- tofu/spectro/_rockingcurve_def.py | 466 ++++++++++++++---------------- 1 file changed, 215 insertions(+), 251 deletions(-) diff --git a/tofu/spectro/_rockingcurve_def.py b/tofu/spectro/_rockingcurve_def.py index 0feb3b16d..90267ebf6 100644 --- a/tofu/spectro/_rockingcurve_def.py +++ b/tofu/spectro/_rockingcurve_def.py @@ -1,7 +1,8 @@ import numpy as np import scipy.interpolate - +import scipy.constants as scpct +import astropy.units as u # ############################################################################# # ############################################################################# @@ -13,13 +14,23 @@ 'Quartz_110': { 'name': '110-Qz', 'symbol': 'Qz', - 'target': 'ArXVII ion (3.96 A)', + 'target': { + 'ion': 'ArXVII', + 'wavelength': 3.96e-10, + 'unit': u.m, + }, 'atoms': ['Si', 'O'], 'atoms_Z': [14., 8.], 'atoms_nb': [3., 6.], 'miller': np.r_[1., 1., 0.], - 'volume': None, - 'd_hkl': None, + 'volume': { + 'value': None, + 'unit': 1/u.m**3, + }, + 'd_hkl': { + 'value': None, + 'unit': u.m, + }, 'mesh': { 'type': 'hexagonal', 'positions': { @@ -46,14 +57,14 @@ }, 'inter_atomic': { 'distances': { - 'a0': 4.91304, - 'c0': 5.40463, + 'a0': 4.91304e-10, + 'c0': 5.40463e-10, }, - 'unit': 'A', + 'unit': u.m, 'comments': 'within the unit cell', 'Tref': { 'data': 25., - 'unit': 'C', + 'unit': u.deg_C, }, 'sources': 'R.W.G. Wyckoff, Crystal Structures', }, @@ -62,18 +73,21 @@ 'alpha_a': 1.337e-5, 'alpha_c': 7.97e-6, }, - 'unit': '1/C', + 'unit': 1/u.deg_C, 'comments': 'in parallel directions to a0 and c0', 'sources': 'R.W.G. Wyckoff, Crystal Structures', }, 'sin_theta_lambda': { - 'Si': np.r_[ - 0., 0.1, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, - 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, - ], - 'O': np.r_[ - 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, - ], + 'values': { + 'Si': np.r_[ + 0., 0.1, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, + 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, + ]*1e10, + 'O': np.r_[ + 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, + ]*1e10, + }, + 'unit': 1/u.m, 'sources': 'Int. Tab. X-Ray Crystallography, Vol.I,II,III,IV (1985)', }, @@ -95,13 +109,23 @@ 'Quartz_102': { 'name': '102-Qz', 'symbol': 'Qz', - 'target': 'ArXVIII ion (3.75 A)', + 'target': { + 'ion': 'ArXVIII', + 'wavelength': 3.75e-10, + 'unit': u.m, + }, 'atoms': ['Si', 'O'], 'atoms_Z': [14., 8.], 'atoms_nb': [3., 6.], 'miller': np.r_[1., 0., 2.], - 'volume': None, - 'd_hkl': None, + 'volume': { + 'value': None, + 'unit': 1/u.m**3, + }, + 'd_hkl': { + 'value': None, + 'unit': u.m, + }, 'mesh': { 'type': 'hexagonal', 'positions': { @@ -128,14 +152,14 @@ }, 'inter_atomic': { 'distances': { - 'a0': 4.91304, - 'c0': 5.40463, + 'a0': 4.91304e-10, + 'c0': 5.40463e-10, }, - 'unit': 'A', + 'unit': u.m, 'comments': 'within the unit cell', 'Tref': { 'data': 25., - 'unit': 'C', + 'unit': u.deg_C, }, 'sources': 'R.W.G. Wyckoff, Crystal Structures', }, @@ -144,7 +168,7 @@ 'alpha_a': 1.337e-5, 'alpha_c': 7.97e-6, }, - 'unit': '1/C', + 'unit': 1/u.deg_C, 'comments': 'in parallel directions to a0 and c0', 'sources': 'R.W.G. Wyckoff, Crystal Structures', }, @@ -153,11 +177,12 @@ 'Si': np.r_[ 0., 0.1, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, - ], + ]*1e10, 'O': np.r_[ 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, - ], + ]*1e10, }, + 'unit': 1/u.m, 'sources': 'Int. Tab. X-Ray Crystallography, Vol.I,II,III,IV (1985)', }, @@ -176,16 +201,27 @@ 'Int. Tab. X-Ray Crystallography, Vol.I,II,III,IV (1985)', }, }, + 'Germanium_XXX': { 'name': None, 'symbol': None, - 'target': None, + 'target': { + 'ion': None, + 'wavelength': None, + 'unit': u.m, + }, 'atoms': None, 'atoms_Z': None, 'atoms_nb': None, 'miller': None, - 'volume': None, - 'd_hkl': None, + 'volume': { + 'value': None, + 'unit': 1/u.m**3, + }, + 'd_hkl': { + 'value': None, + 'unit': u.m, + }, 'mesh': { 'type': None, 'positions': None, @@ -194,19 +230,23 @@ 'phases': None, 'inter_atomic': { 'distances': None, - 'unit': None, + 'unit': u.m, 'comments': None, - 'Tref': None, + 'Tref': { + 'data': None, + 'unit': u.deg_C, + }, 'sources': None, }, 'thermal_expansion': { 'coefs': None, - 'unit': None, + 'unit': 1/u.deg_C, 'comments': None, 'sources': None, }, 'sin_theta_lambda': { 'values': None, + 'unit': 1/u.m, 'sources': None, }, 'atomic_scattering': { @@ -237,100 +277,63 @@ # Atoms positions for Germanium crystal -# --------------------------------------------------------------- -# Attribution to alpha-Quartz crystals: Quartz_110 and Quartz_102 -# --------------------------------------------------------------- - -# Position of the 3 Si atoms in the unit cell -# ------------------------------------------- - -# Quartz_110 -uSi = _DCRYST['Quartz_110']['mesh']['positions']['Si']['u'][0] -_DCRYST['Quartz_110']['mesh']['positions']['Si']['x'] = np.r_[ - -uSi, - uSi, - 0. -] -_DCRYST['Quartz_110']['mesh']['positions']['Si']['y'] = np.r_[ - -uSi, - 0., - uSi -] -_DCRYST['Quartz_110']['mesh']['positions']['Si']['z'] = np.r_[ - 1./3., - 0., - 2./3. -] -_DCRYST['Quartz_110']['mesh']['positions']['Si']['N'] = np.size( - _DCRYST['Quartz_110']['mesh']['positions']['Si']['x'] -) - -# Quartz_102 -_DCRYST['Quartz_102']['mesh']['positions']['Si']['x'] = ( - _DCRYST['Quartz_110']['mesh']['positions']['Si']['x'] -) -_DCRYST['Quartz_102']['mesh']['positions']['Si']['y'] = ( - _DCRYST['Quartz_110']['mesh']['positions']['Si']['y'] -) -_DCRYST['Quartz_102']['mesh']['positions']['Si']['z'] = ( - _DCRYST['Quartz_110']['mesh']['positions']['Si']['z'] -) -_DCRYST['Quartz_102']['mesh']['positions']['Si']['N'] = ( - _DCRYST['Quartz_110']['mesh']['positions']['Si']['N'] -) - -# Position of the 6 O atoms in the unit cell -# ------------------------------------------ - -# Quartz_110 -uOx = _DCRYST['Quartz_110']['mesh']['positions']['O']['u'][0] -uOy = _DCRYST['Quartz_110']['mesh']['positions']['O']['u'][1] -uOz = _DCRYST['Quartz_110']['mesh']['positions']['O']['u'][2] -_DCRYST['Quartz_110']['mesh']['positions']['O']['x'] = np.r_[ - uOx, - uOy - uOx, - -uOy, - uOx - uOy, - uOy, - -uOx -] -_DCRYST['Quartz_110']['mesh']['positions']['O']['y'] = np.r_[ - uOy, - -uOx, - uOx - uOy, - -uOy, - uOx, - uOy - uOx -] -_DCRYST['Quartz_110']['mesh']['positions']['O']['z'] = np.r_[ - uOz, - uOz + 1./3., - uOz + 2./3., - -uOz, - 2./3. - uOz, - 1./3. - uOz -] -_DCRYST['Quartz_110']['mesh']['positions']['O']['N'] = np.size( - _DCRYST['Quartz_110']['mesh']['positions']['O']['x'] -) - -# Quartz_102 -_DCRYST['Quartz_102']['mesh']['positions']['O']['x'] = ( - _DCRYST['Quartz_110']['mesh']['positions']['O']['x'] -) -_DCRYST['Quartz_102']['mesh']['positions']['O']['y'] = ( - _DCRYST['Quartz_110']['mesh']['positions']['O']['y'] -) -_DCRYST['Quartz_102']['mesh']['positions']['O']['z'] = ( - _DCRYST['Quartz_110']['mesh']['positions']['O']['z'] -) -_DCRYST['Quartz_102']['mesh']['positions']['O']['N'] = ( - _DCRYST['Quartz_110']['mesh']['positions']['O']['N'] -) - -# --------------------------------- -# Attribution to Germanium crystals -# --------------------------------- +# ----------- +# Attribution +# ----------- + +for ii in _DCRYST.keys(): + if "Quartz" in ii: + uSi = _DCRYST[ii]['mesh']['positions']['Si']['u'][0] + uOx = _DCRYST[ii]['mesh']['positions']['O']['u'][0] + uOy = _DCRYST[ii]['mesh']['positions']['O']['u'][1] + uOz = _DCRYST[ii]['mesh']['positions']['O']['u'][2] + _DCRYST[ii]['mesh']['positions']['Si']['x'] = np.r_[ + -uSi, + uSi, + 0. + ] + _DCRYST[ii]['mesh']['positions']['Si']['y'] = np.r_[ + -uSi, + 0., + uSi + ] + _DCRYST[ii]['mesh']['positions']['Si']['z'] = np.r_[ + 1./3., + 0., + 2./3. + ] + _DCRYST[ii]['mesh']['positions']['Si']['N'] = np.size( + _DCRYST[ii]['mesh']['positions']['Si']['x'] + ) + _DCRYST[ii]['mesh']['positions']['O']['x'] = np.r_[ + uOx, + uOy - uOx, + -uOy, + uOx - uOy, + uOy, + -uOx + ] + _DCRYST[ii]['mesh']['positions']['O']['y'] = np.r_[ + uOy, + -uOx, + uOx - uOy, + -uOy, + uOx, + uOy - uOx + ] + _DCRYST[ii]['mesh']['positions']['O']['z'] = np.r_[ + uOz, + uOz + 1./3., + uOz + 2./3., + -uOz, + 2./3. - uOz, + 1./3. - uOz + ] + _DCRYST[ii]['mesh']['positions']['O']['N'] = np.size( + _DCRYST[ii]['mesh']['positions']['O']['x'] + ) + elif "Ge" in ii: + None # ############################################################################# # ############################################################################# @@ -339,47 +342,34 @@ # ############################################################################# # ############################################################################# -# ----------------------------------------------------------------------- -# Definition of volume and inter-reticular spacing relations, func(meshtype) -# ----------------------------------------------------------------------- - +# -------------------------------------------------------------------------- +# Definition of volume and inter-reticular spacing, function of the meshtype +# -------------------------------------------------------------------------- def hexa_volume(aa, cc): return (aa**2) * cc * (np.sqrt(3.)/2.) - def hexa_spacing(hh, kk, ll, aa, cc): return np.sqrt( (3.*(aa**2)*(cc**2)) / (4.*(hh**2 + kk**2 + hh*kk)*(cc**2) + 3.*(ll**2)*(aa**2)) ) - -# --------------------------------------------------------------- -# Attribution to alpha-Quartz crystals: Quartz_110 and Quartz_102 -# --------------------------------------------------------------- - - -# Same values for 110- and Quartz_102 -a = _DCRYST['Quartz_110']['inter_atomic']['distances']['a0'] -c = _DCRYST['Quartz_110']['inter_atomic']['distances']['c0'] - -h110 = _DCRYST['Quartz_110']['miller'][0] -k110 = _DCRYST['Quartz_110']['miller'][1] -l110 = _DCRYST['Quartz_110']['miller'][2] -h102 = _DCRYST['Quartz_102']['miller'][0] -k102 = _DCRYST['Quartz_102']['miller'][1] -l102 = _DCRYST['Quartz_102']['miller'][2] - -_DCRYST['Quartz_110']['volume'] = hexa_volume(a, c) -_DCRYST['Quartz_110']['d_hkl'] = hexa_spacing(h110, k110, l110, a, c) -_DCRYST['Quartz_102']['volume'] = hexa_volume(a, c) -_DCRYST['Quartz_102']['d_hkl'] = hexa_spacing(h102, k102, l102, a, c) - - -# --------------------------------- -# Attribution to Germanium crystals -# --------------------------------- +# ----------- +# Attribution +# ----------- + +for ii in _DCRYST.keys(): + if "Quartz" in ii and "hexagonal" in _DCRYST[ii]['mesh']['type']: + a = _DCRYST[ii]['inter_atomic']['distances']['a0'] # m + c = _DCRYST[ii]['inter_atomic']['distances']['c0'] # m + h = _DCRYST[ii]['miller'][0] + k = _DCRYST[ii]['miller'][1] + l = _DCRYST[ii]['miller'][2] + _DCRYST[ii]['volume']['value'] = hexa_volume(a, c) + _DCRYST[ii]['d_hkl']['value'] = hexa_spacing(h, k, l, a, c) + elif "Ge" in ii: + None # ############################################################################# # ############################################################################# @@ -387,117 +377,91 @@ def hexa_spacing(hh, kk, ll, aa, cc): # ############################################################################# # ############################################################################# -# --------------------------------------------------------------- -# Attribution to alpha-Quartz crystals: Quartz_110 and Quartz_102 -# --------------------------------------------------------------- +# ----------- +# Attribution +# ----------- # From W. Zachariasen, Theory of X-ray Diffraction in Crystals # (Wiley, New York, 1945) # Linear absorption coefficient # ----------------------------- -# Same values for 110- and Quartz_102 -Zsi = _DCRYST['Quartz_110']['atoms_Z'][0] -Zo = _DCRYST['Quartz_110']['atoms_Z'][1] - - -def mu_si(lamb): - return 1.38e-2*(lamb**2.79)*(Zsi**2.73) - - -def mu_si1(lamb): - return 5.33e-4*(lamb**2.74)*(Zsi**3.03) - - -def mu_o(lamb): - return 5.4e-3*(lamb**2.92)*(Zo**3.07) - - -def mu(lamb, mu_si, mu_o): - return 2.65e-8*(7.*mu_si + 8.*mu_o)/15. - +for ii in _DCRYST.keys(): + if "Quartz" in ii: + Zsi, Zo = _DCRYST[ii]['atoms_Z'][0], _DCRYST[ii]['atoms_Z'][1] + # For wavelength in Angstroms (ex: lamb=3.96A) + def mu_si(lamb): + return 1.38e-2*(lamb**2.79)*(Zsi**2.73) + def mu_si1(lamb): + return 5.33e-4*(lamb**2.74)*(Zsi**3.03) + def mu_o(lamb): + return 5.4e-3*(lamb**2.92)*(Zo**3.07) + def mu(lamb, mu_si, mu_o): + return 2.65e-8*(7.*mu_si + 8.*mu_o)/15. # Atomic scattering factor, real and imaginary parts # -------------------------------------------------- - -# Same values for 110- and Quartz_102 -sol_si = _DCRYST['Quartz_110']['sin_theta_lambda']['Si'] -sol_o = _DCRYST['Quartz_110']['sin_theta_lambda']['O'] -asf_si = _DCRYST['Quartz_110']['atomic_scattering']['factors']['Si'] -asf_o = _DCRYST['Quartz_110']['atomic_scattering']['factors']['O'] -interp_si = scipy.interpolate.interp1d(sol_si, asf_si) -interp_o = scipy.interpolate.interp1d(sol_o, asf_o) - - -def dfsi_re(lamb): - return 0.1335*lamb - 6e-3 - - -def fsi_re(lamb, sol): - return interp_si(sol) + dfsi_re(lamb) - - -def fsi_im(lamb, mu_si): - return 5.936e-4*Zsi*(mu_si/lamb) - - -def dfo_re(lamb): - return 0.1335*lamb - 0.206 - - -def fo_re(lamb, sol): - return interp_o(sol) + dfo_re(lamb) - - -def fo_im(lamb, mu_o): - return 5.936e-4*Zo*(mu_o/lamb) - +for ii in _DCRYST.keys(): + if "Quartz" in ii: + sol_si = _DCRYST[ii]['sin_theta_lambda']['values']['Si'] + sol_o = _DCRYST[ii]['sin_theta_lambda']['values']['O'] + asf_si = _DCRYST[ii]['atomic_scattering']['factors']['Si'] + asf_o = _DCRYST[ii]['atomic_scattering']['factors']['O'] + interp_si = scipy.interpolate.interp1d(sol_si, asf_si) + interp_o = scipy.interpolate.interp1d(sol_o, asf_o) + + # For wavelength in Angstroms (ex: lamb=3.96A) + def dfsi_re(lamb): + return 0.1335*lamb - 6e-3 + def fsi_re(lamb, sol): + return interp_si(sol) + dfsi_re(lamb) + def fsi_im(lamb, mu_si): + return 5.936e-4*Zsi*(mu_si/lamb) + def dfo_re(lamb): + return 0.1335*lamb - 0.206 + def fo_re(lamb, sol): + return interp_o(sol) + dfo_re(lamb) + def fo_im(lamb, mu_o): + return 5.936e-4*Zo*(mu_o/lamb) # Phases # ------ - -h110 = _DCRYST['Quartz_110']['miller'][0] -k110 = _DCRYST['Quartz_110']['miller'][1] -l110 = _DCRYST['Quartz_110']['miller'][2] -h102 = _DCRYST['Quartz_102']['miller'][0] -k102 = _DCRYST['Quartz_102']['miller'][1] -l102 = _DCRYST['Quartz_102']['miller'][2] - - -# Same values for 110- and Quartz_102 -xsi = _DCRYST['Quartz_110']['mesh']['positions']['Si']['x'] -ysi = _DCRYST['Quartz_110']['mesh']['positions']['Si']['y'] -zsi = _DCRYST['Quartz_110']['mesh']['positions']['Si']['z'] -Nsi = _DCRYST['Quartz_110']['mesh']['positions']['Si']['N'] -xo = _DCRYST['Quartz_110']['mesh']['positions']['O']['x'] -yo = _DCRYST['Quartz_110']['mesh']['positions']['O']['y'] -zo = _DCRYST['Quartz_110']['mesh']['positions']['O']['z'] -No = _DCRYST['Quartz_110']['mesh']['positions']['O']['N'] - - def phasesi(hh, kk, ll, xsi, ysi, zsi): return hh*xsi + kk*ysi + ll*zsi - - def phaseo(hh, kk, ll, xo, yo, zo): return hh*xo + kk*yo + ll*zo -phaseSi_110 = np.full((Nsi), np.nan) -phaseO_110 = np.full((No), np.nan) -phaseSi_102 = np.full((Nsi), np.nan) -phaseO_102 = np.full((No), np.nan) - -for i in range(Nsi): - phaseSi_110[i] = phasesi(h110, k110, l110, xsi[i], ysi[i], zsi[i]) - phaseSi_102[i] = phasesi(h102, k102, l102, xsi[i], ysi[i], zsi[i]) -_DCRYST['Quartz_110']['phases']['Si'] = phaseSi_110 -_DCRYST['Quartz_102']['phases']['Si'] = phaseSi_102 - -for i in range(No): - phaseO_110[i] = phaseo(h110, k110, l110, xo[i], yo[i], zo[i]) - phaseO_102[i] = phaseo(h102, k102, l102, xo[i], yo[i], zo[i]) -_DCRYST['Quartz_110']['phases']['O'] = phaseO_110 -_DCRYST['Quartz_102']['phases']['O'] = phaseO_102 +for ii in _DCRYST.keys(): + if "Quartz" in ii: + h = _DCRYST[ii]['miller'][0] + k = _DCRYST[ii]['miller'][1] + l = _DCRYST[ii]['miller'][2] + xsi = _DCRYST[ii]['mesh']['positions']['Si']['x'] + ysi = _DCRYST[ii]['mesh']['positions']['Si']['y'] + zsi = _DCRYST[ii]['mesh']['positions']['Si']['z'] + Nsi = _DCRYST[ii]['mesh']['positions']['Si']['N'] + xo = _DCRYST[ii]['mesh']['positions']['O']['x'] + yo = _DCRYST[ii]['mesh']['positions']['O']['y'] + zo = _DCRYST[ii]['mesh']['positions']['O']['z'] + No = _DCRYST[ii]['mesh']['positions']['O']['N'] + + phaseSi = np.full((Nsi), np.nan) + phaseO = np.full((No), np.nan) + + if "110" in ii: + for i in range(Nsi): + phaseSi[i] = phasesi(h, k, l, xsi[i], ysi[i], zsi[i]) + for i in range(No): + phaseO[i] = phaseo(h, k, l, xo[i], yo[i], zo[i]) + _DCRYST[ii]['phases']['Si'] = phaseSi + _DCRYST[ii]['phases']['O'] = phaseO + elif "102" in ii: + for i in range(Nsi): + phaseSi[i] = phasesi(h, k, l, xsi[i], ysi[i], zsi[i]) + for i in range(No): + phaseO[i] = phaseo(h, k, l, xo[i], yo[i], zo[i]) + _DCRYST[ii]['phases']['Si'] = phaseSi + _DCRYST[ii]['phases']['O'] = phaseO From 81447bb6159c1c35d487b78aca395b7e604e5baa Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 24 Aug 2022 18:27:00 +0200 Subject: [PATCH 02/14] [#673] Correcting module name (Polygon -> polygon) --- tofu/geom/_etendue.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tofu/geom/_etendue.py b/tofu/geom/_etendue.py index 722a1c6a5..d8dfa6c84 100644 --- a/tofu/geom/_etendue.py +++ b/tofu/geom/_etendue.py @@ -9,7 +9,7 @@ import matplotlib.lines as mlines -import Polygon as plg +import polygon as plg import datastock as ds From 820a1ab6f0f5621e70accd1b6381d29151d04a5c Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 24 Aug 2022 18:31:25 +0200 Subject: [PATCH 03/14] [#673] Add astropy.units module into tests routines, used in tofu/spectro/_rockingcurve_def.py --- tofu/tests/tests00_root/test_03_plot.py | 2 +- tofu/tests/tests04_spectro/test_03_rockingcurve.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/tofu/tests/tests00_root/test_03_plot.py b/tofu/tests/tests00_root/test_03_plot.py index 339818c01..fb7a0023c 100644 --- a/tofu/tests/tests00_root/test_03_plot.py +++ b/tofu/tests/tests00_root/test_03_plot.py @@ -9,7 +9,7 @@ import numpy as np import matplotlib.pyplot as plt import warnings as warn - +import astropy.units as u # Importing package tofu.geom import tofu as tf diff --git a/tofu/tests/tests04_spectro/test_03_rockingcurve.py b/tofu/tests/tests04_spectro/test_03_rockingcurve.py index fd83e4c80..d5b07d0e2 100644 --- a/tofu/tests/tests04_spectro/test_03_rockingcurve.py +++ b/tofu/tests/tests04_spectro/test_03_rockingcurve.py @@ -11,6 +11,7 @@ # Standard import numpy as np import scipy.constants as scpct +import astropy.units as u import matplotlib.pyplot as plt # tofu-specific @@ -112,7 +113,7 @@ def test01_rockingcurve_compute(self): for k0 in self.lc: dout = tfs.compute_rockingcurve( crystal=k0, - lamb=np.r_[3.969067], + lamb=np.r_[3.969067e-10], use_non_parallelism=False, alpha_limits=None, therm_exp=False, From 5d439d7c88e99d7316747f85ed8891d266722615 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Thu, 25 Aug 2022 15:45:42 +0200 Subject: [PATCH 04/14] [#673] Correcting routines to take into account values of wavelength in meters --- tofu/geom/_core_optics.py | 98 +++++++++++++++------ tofu/spectro/_rockingcurve.py | 138 +++++++++++++++++------------- tofu/spectro/_rockingcurve_def.py | 33 +------ 3 files changed, 155 insertions(+), 114 deletions(-) diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index 957984a71..f9ab42ccf 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -809,7 +809,7 @@ def compute_rockingcurve( crystal=None, din=None, lamb=None, miscut=None, nn=None, - alpha_limits=None, + miscut_limits=None, therm_exp=None, temp_limits=None, plot_therm_exp=None, @@ -821,7 +821,7 @@ def compute_rockingcurve( crystal=crystal, din=din, lamb=lamb, miscut=miscut, nn=nn, - alpha_limits=alpha_limits, + miscut_limits=miscut_limits, therm_exp=therm_exp, temp_limits=temp_limits, plot_therm_exp=plot_therm_exp, @@ -833,7 +833,7 @@ def compute_rockingcurve( def plot_var_temp_changes_wavelengths( self, ih=None, ik=None, il=None, lambdas=None, miscut=None, na=None, - alpha_limits=None, + miscut_limits=None, therm_exp=None, plot_therm_exp=None, plot_asf=None, plot_power_ratio=None, plot_asymmetry=None, plot_cmaps=None, @@ -843,7 +843,7 @@ def plot_var_temp_changes_wavelengths( return _rockingcurve.plot_var_temp_changes_wavelengths( ih=ih, ik=ik, il=il, lambdas=lambdas, miscut=miscut, na=na, - alpha_limits=alpha_limits, + miscut_limits=miscut_limits, therm_exp=therm_exp, plot_therm_exp=plot_therm_exp, plot_asf=plot_asf, plot_power_ratio=plot_power_ratio, plot_asymmetry=plot_asymmetry, plot_cmaps=plot_cmaps, @@ -1944,6 +1944,26 @@ def calc_xixj_from_braggphi( else: return xi, xj + def _checkformat_xi_values( + self, + xi_tab=None, + bragg=None, + lamb=None, + phi=None, + det=None, + ): + if np.all(np.isnan(xi_tab)) is True: + msg = ( + "All coordinates xi and xj have NaN values !\n" + "Check arguments for calc_xixj_from_braggphi():\n" + + "\t - xi_tab: {}\n".format(xi_tab) + + "\t - bragg: {}\n".format(bragg) + + "\t - lamb: {}\n".format(lamb) + + "\t - phi: {}\n".format(phi) + + "\t - det: {}\n".format(det) + ) + raise Exception(msg) + def plot_line_on_det_tracing( self, # Options of basic method @@ -1959,7 +1979,7 @@ def plot_line_on_det_tracing( merge_rc_data=None, miscut=None, therm_exp=None, - alpha_limits=None, na=None, + miscut_limits=None, na=None, alpha0=None, temp0=None, temp_limits=None, # Plot @@ -1979,7 +1999,7 @@ def plot_line_on_det_tracing( - merge_rc_data: bool use tf/spectro/_rockingucurve.py to plot in transparency ranges the angular extent of each wavelength traces - - alpha_limits: array + - miscut_limits: array asymmetry angle range, provide only both limits. By default in tf/spectro/_rockingucurve.py between +/-5 arcmin - na: float @@ -1989,7 +2009,7 @@ def plot_line_on_det_tracing( - alpha0: float Wanted value in radians of the amplitude miscut angle. By default to 3 arcmin = 0.05 deg = pi/3600 rad - '0' for alpha_limits[0], 'na-1' for alpha_limits[1]. + '0' for miscut_limits[0], 'na-1' for miscut_limits[1]. By default to 3 arcmin = 0.05 rad so alpha0=40 - temp0: float Wanted value of the temperature change. @@ -2032,8 +2052,8 @@ def plot_line_on_det_tracing( johann = lpsi is not None or ldtheta is not None if rocking is None: rocking = False - if alpha_limits is None: - alpha_limits = np.r_[-(3/60)*np.pi/180, (3/60)*np.pi/180] + if miscut_limits is None: + miscut_limits = np.r_[-(3/60)*np.pi/180, (3/60)*np.pi/180] if temp_limits is None: temp_limits = np.r_[-10, 10, 25] if na is None: @@ -2062,10 +2082,10 @@ def plot_line_on_det_tracing( self.update_miscut(alpha=0., beta=0.) if miscut: self.update_miscut(alpha=alpha0, beta=0.) - # T0, TD, a1, c1, Volume, d_atom, sol, sin_theta, theta, theta_deg, + dout = _rockingcurve.CrystBragg_comp_lattice_spacing( crystal=crystal, din=din, - lamb=self.dbragg['lambref']*1e10, + lamb=self.dbragg['lambref'], na=na, nn=nn, therm_exp=therm_exp, temp_limits=temp_limits, @@ -2074,8 +2094,8 @@ def plot_line_on_det_tracing( T0 = dout['Temperature of reference (°C)'] TD = dout['Temperature variations (°C)'] Volume = dout['Volume (1/m3)'] - d_atom = dout['Inter-reticular spacing (A)'] - sol = dout['sinus over lambda'] + d_atom = dout['Inter-reticular spacing (m)'] + sol = dout['sinus_theta_lambda'] theta = dout['theta_Bragg (rad)'] theta_deg = dout['theta_Bragg (deg)'] @@ -2085,7 +2105,7 @@ def find_nearest(array, value): return idx id_temp0 = find_nearest(TD, temp0) - self.dmat['d'] = d_atom[id_temp0]*1e-10 + self.dmat['d'] = d_atom[id_temp0] # Get local basis nout, e1, e2, miscut = self.get_unit_vectors( @@ -2105,6 +2125,7 @@ def find_nearest(array, value): phimin, phimax = np.nanmin(phi), np.nanmax(phi) phimin, phimax = phimin-(phimax-phimin)/10, phimax+(phimax-phimin)/10 + # Get reference ray-tracing bragg = self._checkformat_bragglamb(lamb=lamb, n=n) if nphi2 is None: @@ -2126,6 +2147,13 @@ def find_nearest(array, value): strict=strict, plot=False, ) + self._checkformat_xi_values( + xi_tab=xi[:, :], + bragg=bragg, + lamb=lamb, + phi=phi, + det=det, + ) # Get johann-error raytracing (multiple positions on crystal) xi_er, xj_er = None, None @@ -2149,12 +2177,23 @@ def find_nearest(array, value): for ii in range(npsi): i0 = np.arange(ii*nphi, (ii+1)*nphi) xi_er[l, i0], xj_er[l, i0] = self.calc_xixj_from_braggphi( - phi=phi, bragg=bragg[l], lamb=None, n=n, - dtheta=ldtheta[ii], psi=lpsi[ii], - det=det, plot=False, + bragg=bragg[l], + phi=phi, + dtheta=ldtheta[ii], + psi=lpsi[ii], + n=n, + det=det, miscut=miscut, strict=strict, + plot=False, ) + self._checkformat_xi_values( + xi_tab=xi_er[:, :], + bragg=bragg, + lamb=lamb, + phi=phi, + det=det, + ) # Get rocking curve error if rocking: @@ -2179,12 +2218,12 @@ def find_nearest(array, value): for ll in range(nlamb): dout = _rockingcurve.compute_rockingcurve( crystal=crystal, din=din, - lamb=lamb[ll]*1e10, + lamb=lamb[ll], miscut=miscut, therm_exp=therm_exp, temp_limits=temp_limits, plot_therm_exp=plot_rcs, - alpha_limits=alpha_limits, nn=None, + miscut_limits=miscut_limits, nn=None, plot_asf=False, plot_power_ratio=plot_rcs, plot_asymmetry=False, plot_cmaps=False, returnas=dict, @@ -2278,9 +2317,18 @@ def find_nearest(array, value): strict=strict, plot=False, ) + self._checkformat_xi_values( + xi_tab=xi_rc[:, :, :], + bragg=bragg, + lamb=lamb, + phi=phi, + det=det, + ) + xi_atprmax[ll] = xi_rc[ll, ind_pr_max, nphi2] xj_atprmax[ll] = xj_rc[ll, ind_pr_max, nphi2] self.update_miscut(alpha=0., beta=0.) + if therm_exp: self.dmat['d'] = d_atom[nn]*1e-10 else: @@ -2363,7 +2411,7 @@ def comp_angular_shift_on_det_tracing( ih=None, ik=None, il=None, dcryst=None, merge_rc_data=None, - therm_exp=None, alpha_limits=None, na=None, + therm_exp=None, miscut_limits=None, na=None, temp=None, plot=None, ax=None, dleg=None, color=None, @@ -2381,7 +2429,7 @@ def comp_angular_shift_on_det_tracing( - merge_rc_data: bool use tf/spectro/_rockingucurve.py to plot in transparency ranges the angular extent of each wavelength traces - - alpha_limits: array + - miscut_limits: array asymmetry angle range, provide only both limits. By default in tf/spectro/_rockingucurve.py between +/-5 arcmin - na: float @@ -2391,7 +2439,7 @@ def comp_angular_shift_on_det_tracing( - alpha0: float Wanted value in radians of the amplitude miscut angle. By default to 3 arcmin = 0.05 deg = pi/3600 rad - '0' for alpha_limits[0], 'na-1' for alpha_limits[1]. + '0' for miscut_limits[0], 'na-1' for miscut_limits[1]. By default to 3 arcmin = 0.05 rad so alpha0=40 - temp0: float Wanted value of the temperature change. @@ -2415,8 +2463,8 @@ def comp_angular_shift_on_det_tracing( miscut = True if therm_exp is None: therm_exp = True - if alpha_limits is None: - alpha_limits = np.r_[-(3/60)*np.pi/180, (3/60)*np.pi/180] + if miscut_limits is None: + miscut_limits = np.r_[-(3/60)*np.pi/180, (3/60)*np.pi/180] if temp is None: temp = np.r_[-10, +10] if na is None: @@ -2433,7 +2481,7 @@ def comp_angular_shift_on_det_tracing( plot = False # - angles = np.linspace(alpha_limits[0], alpha_limits[1], na) + angles = np.linspace(miscut_limits[0], miscut_limits[1], na) TD = np.linspace(temp[0], temp[1], na) # Computation of angular shifts diff --git a/tofu/spectro/_rockingcurve.py b/tofu/spectro/_rockingcurve.py index f8dd0a94c..76c91d89b 100644 --- a/tofu/spectro/_rockingcurve.py +++ b/tofu/spectro/_rockingcurve.py @@ -32,7 +32,7 @@ def compute_rockingcurve( lamb=None, # Lattice modifications miscut=None, nn=None, - alpha_limits=None, + miscut_limits=None, therm_exp=None, temp_limits=None, # Plot @@ -93,10 +93,10 @@ def compute_rockingcurve( and soon 'Ge' lamb: float Wavelength of interest, in Angstroms (1e-10 m) - ex: lamb=np.r_[3.96] + ex: lamb=np.r_[3.96e-10] miscut: str Introduce miscut between dioptre and reflecting planes - alpha_limits: array + miscut_limits: array Asymmetry angle range. Provide only both boundary limits Ex: np.r_[-3, 3] in radians nn: int @@ -154,8 +154,8 @@ def compute_rockingcurve( plot_therm_exp = True if miscut is None: miscut = False - if alpha_limits is None: - alpha_limits = np.r_[-(3/60)*np.pi/180, (3/60)*np.pi/180] + if miscut_limits is None: + miscut_limits = np.r_[-(3/60)*np.pi/180, (3/60)*np.pi/180] if temp_limits is None: temp_limits = np.r_[-10, 10, 25] if nn is None: @@ -178,6 +178,7 @@ def compute_rockingcurve( ik=din['miller'][1], il=din['miller'][2], lamb=lamb, + din=din, ) # Classical electronical radius, in Angstroms, from the NIST Reference on @@ -197,15 +198,15 @@ def compute_rockingcurve( T0 = dout['Temperature of reference (°C)'] TD = dout['Temperature variations (°C)'] Volume = dout['Volume (1/m3)'] - d_atom = dout['Inter-reticular spacing (A)'] - sol = dout['sinus over lambda'] + d_atom = dout['Inter-reticular spacing (m)'] + sol = dout['sinus_theta_lambda'] theta = dout['theta_Bragg (rad)'] theta_deg = dout['theta_Bragg (deg)'] # Check validity of asymmetry angle alpha limits in arguments alpha, bb = CrystBragg_check_alpha_angle( theta=theta, - alpha_limits=alpha_limits, na=na, nn=nn, + miscut_limits=miscut_limits, na=na, nn=nn, miscut=miscut, therm_exp=therm_exp, ) @@ -218,25 +219,25 @@ def compute_rockingcurve( # ----------------------------------- # Atomic absorption coefficient - if lamb > 6.74: - mu_si = _rockingcurve_def.mu_si1(lamb=lamb) + if lamb > 6.74e10: + mu_si = _rockingcurve_def.mu_si1(lamb=lamb*1e10) else: - mu_si = _rockingcurve_def.mu_si(lamb=lamb) - mu_o = _rockingcurve_def.mu_o(lamb=lamb) - mu = _rockingcurve_def.mu(lamb=lamb, mu_si=mu_si, mu_o=mu_o) + mu_si = _rockingcurve_def.mu_si(lamb=lamb*1e10) + mu_o = _rockingcurve_def.mu_o(lamb=lamb*1e10) + mu = _rockingcurve_def.mu(lamb=lamb*1e10, mu_si=mu_si, mu_o=mu_o) # Atomic scattering factor ("f") in function of sol # ("_re") for the real part and ("_im") for the imaginary part fsi_re = np.full((sol.size), np.nan) fo_re = np.full((sol.size), np.nan) - dfsi_re = _rockingcurve_def.dfsi_re(lamb=lamb) - dfo_re = _rockingcurve_def.dfo_re(lamb=lamb) + dfsi_re = _rockingcurve_def.dfsi_re(lamb=lamb*1e10) + dfo_re = _rockingcurve_def.dfo_re(lamb=lamb*1e10) for i in range(sol.size): - fsi_re[i] = _rockingcurve_def.fsi_re(lamb=lamb, sol=sol[i]) - fo_re[i] = _rockingcurve_def.fo_re(lamb=lamb, sol=sol[i]) - fsi_im = _rockingcurve_def.fsi_im(lamb=lamb, mu_si=mu_si) - fo_im = _rockingcurve_def.fo_im(lamb=lamb, mu_o=mu_o) + fsi_re[i] = _rockingcurve_def.fsi_re(lamb=lamb*1e10, sol=sol[i]) + fo_re[i] = _rockingcurve_def.fo_re(lamb=lamb*1e10, sol=sol[i]) + fsi_im = _rockingcurve_def.fsi_im(lamb=lamb*1e10, mu_si=mu_si) + fo_im = _rockingcurve_def.fo_im(lamb=lamb*1e10, mu_o=mu_o) # Structure factor ("F") for (hkl) reflection # xsi and ih have already been defined with din @@ -314,17 +315,20 @@ def compute_rockingcurve( / (F_re[i]**2.) ) # Real part of psi_H - psi_re[i] = (re*(lamb**2)*F_re[i])/(np.pi*Volume[i]) + # times 1e-10 to convert into Angstroms + psi_re[i] = ( + (re*(lamb**2)*F_re[i])/(np.pi*Volume[i]) + )*1e-10 # Zero-order real part (averaged) - psi0_dre[i] = -re*(lamb**2)*( - No*(Zo + dfo_re) + Nsi*(Zsi + dfsi_re) - )/(np.pi*Volume[i]) + psi0_dre[i] = ( + -re*(lamb**2)*( + No*(Zo + dfo_re) + Nsi*(Zsi + dfsi_re) + )/(np.pi*Volume[i]) + )*1e-10 # Zero-order imaginary part (averaged) psi0_im[i] = ( - -re*(lamb**2) - * (No*fo_im + Nsi*fsi_im) - / (np.pi*Volume[i]) - ) + -re*(lamb**2)*(No*fo_im + Nsi*fsi_im)/(np.pi*Volume[i]) + )*1e-10 # Power ratio and their integrated reflectivity for 3 crystals models: # perfect (Darwin model), ideally mosaic thick and dynamical @@ -374,7 +378,7 @@ def compute_rockingcurve( if plot_power_ratio: CrystalBragg_plot_power_ratio_vs_glancing_angle( din=din, lamb=lamb, - alpha_limits=alpha_limits, + miscut_limits=miscut_limits, theta=theta, theta_deg=theta_deg, th=th, dth=dth, power_ratio=power_ratio, bb=bb, polar=polar, alpha=alpha, @@ -425,9 +429,9 @@ def compute_rockingcurve( det_perp = det_perp[0, 0] dout = { - 'Wavelength (A)': lamb, + 'Wavelength (m)': lamb, 'Miller indices': (ih, ik, il), - 'Inter-reticular distance (A)': d_atom, + 'Inter-reticular distance (m)': d_atom, 'Volume (A^3)': Volume, 'Bragg angle of reference (rad)': theta, 'Glancing angles': dth, @@ -467,7 +471,7 @@ def plot_var_temp_changes_wavelengths( ih=None, ik=None, il=None, lambdas=None, # lattice modifications miscut=None, nn=None, - alpha_limits=None, + miscut_limits=None, therm_exp=None, # Plot plot_therm_exp=None, @@ -538,15 +542,15 @@ def plot_var_temp_changes_wavelengths( dout = compute_rockingcurve( ih=ih, ik=ik, il=il, lamb=lambdas[aa], miscut=miscut, - alpha_limits=alpha_limits, na=na, + miscut_limits=miscut_limits, na=na, therm_exp=therm_exp, plot_therm_exp=False, plot_asf=False, plot_power_ratio=False, plot_asymmetry=False, plot_cmaps=False, verb=False, returnas=dict, ) - din[lambdas[aa]]['Wavelength (A)'] = dout['Wavelength (A)\n'] - din[lambdas[aa]]['Inter-reticular distance (A)'] = ( - dout['Inter-reticular distance (A)\n'] + din[lambdas[aa]]['Wavelength (m)'] = dout['Wavelength (m)\n'] + din[lambdas[aa]]['Inter-reticular distance (m)'] = ( + dout['Inter-reticular distance (m)\n'] ) din[lambdas[aa]]['Bragg angle of reference (rad)'] = ( dout['Bragg angle of reference (rad)\n'] @@ -618,13 +622,13 @@ def plot_var_temp_changes_wavelengths( din[lambdas[aa]][ 'Inter-planar spacing variations (perp. compo)' ][:, nn] - ( - din[lambdas[aa]]['Inter-reticular distance (A)'] - - din[lambdas[aa]]['Inter-reticular distance (A)'][nn] - )/din[lambdas[aa]]['Inter-reticular distance (A)'][nn] + din[lambdas[aa]]['Inter-reticular distance (m)'] + - din[lambdas[aa]]['Inter-reticular distance (m)'][nn] + )/din[lambdas[aa]]['Inter-reticular distance (m)'][nn] ) ax.scatter( din[lambdas[aa]]['Temperature changes (°C)'], - din[lambdas[aa]]['Inter-reticular distance (A)'], + din[lambdas[aa]]['Inter-reticular distance (m)'], marker=markers[aa], c='k', alpha=0.5, label=r'$\lambda$ = ({})$\AA$'.format(lambdas[aa]), ) @@ -725,7 +729,7 @@ def plot_var_temp_changes_wavelengths( def CrystBragg_check_inputs_rockingcurve( - ih=None, ik=None, il=None, lamb=None, + ih=None, ik=None, il=None, lamb=None, din=None, ): dd = {'ih': ih, 'ik': ik, 'il': il, 'lamb': lamb} @@ -769,16 +773,32 @@ def CrystBragg_check_inputs_rockingcurve( ) raise Exception(msg) + # Wavelengths in meters + cdt = [ + lamb <= 0.1*din['target']['wavelength'], + lamb >= 10*din['target']['wavelength'], + ] + if any(cdt): + msg = ( + "Please make sure that the wanted wavelength is close to:\n" + + "\t - targeted wavelength of {} crystal: ({})\n".format( + din['name'], + din['target']['wavelength'], + ) + + "\t - wanted wavelength: ({})\n".format(lamb) + ) + raise Exception(msg) + return ih, ik, il, lamb, def CrystBragg_check_alpha_angle( theta=None, miscut=None, therm_exp=None, - alpha_limits=None, na=None, nn=None, + miscut_limits=None, na=None, nn=None, ): - if alpha_limits is None: + if miscut_limits is None: if not miscut: alpha = np.full((na), 0.) bb = np.full((theta.size, alpha.size), -1.) @@ -807,9 +827,9 @@ def CrystBragg_check_alpha_angle( else: if therm_exp: aa = theta[nn] - lc = [alpha_limits[0] < -aa, alpha_limits[1] > aa] + lc = [miscut_limits[0] < -aa, miscut_limits[1] > aa] else: - lc = [alpha_limits[0] < -theta, alpha_limits[1] > theta] + lc = [miscut_limits[0] < -theta, miscut_limits[1] > theta] if any(lc): msg = ( "Alpha angle limits are not valid!\n" @@ -817,20 +837,20 @@ def CrystBragg_check_alpha_angle( "Limit condition: -({}) < |alpha| < ({}) and\n".format( theta, theta, ) - + "alpha = [({})] rad\n".format(alpha_limits) + + "alpha = [({})] rad\n".format(miscut_limits) ) raise Exception(msg) alpha = np.full((na), np.nan) bb = np.full((theta.size, alpha.size), np.nan) for i in range(theta.size): if therm_exp: - alpha = np.linspace(alpha_limits[0], alpha_limits[1], na) + alpha = np.linspace(miscut_limits[0], miscut_limits[1], na) bb[i, ...] = np.sin(alpha + theta[i])/np.sin( alpha - theta[i] ) else: alpha = np.linspace( - alpha_limits[0], alpha_limits[1], na + miscut_limits[0], miscut_limits[1], na ).reshape(na) bb[i, ...] = np.sin(alpha + theta[i])/np.sin( alpha - theta[i] @@ -895,6 +915,7 @@ def CrystBragg_comp_lattice_spacing( ik=din['miller'][1], il=din['miller'][2], lamb=lamb, + din=din, ) if nn is None: nn = 20 @@ -940,7 +961,6 @@ def CrystBragg_comp_lattice_spacing( # Compute # ------- - for i in range(TD.size): if cond0: a1[i] = a0*(1 + alpha_a*TD[i]) @@ -981,11 +1001,11 @@ def CrystBragg_comp_lattice_spacing( dout = { 'Temperature of reference (°C)': T0, 'Temperature variations (°C)': TD, - 'Inter_atomic distance a1 (A)': a1, - 'Inter_atomic distance c1 (A)': c1, + 'Inter_atomic distance a1 (m)': a1, + 'Inter_atomic distance c1 (m)': c1, 'Volume (1/m3)': Volume, - 'Inter-reticular spacing (A)': d_atom, - 'sinus over lambda': sol, + 'Inter-reticular spacing (m)': d_atom, + 'sinus_theta_lambda': sol, 'sinus theta_Bragg': sin_theta, 'theta_Bragg (rad)': theta, 'theta_Bragg (deg)': theta_deg, @@ -1034,11 +1054,11 @@ def CrystBragg_comp_integrated_reflect( P_per = np.full((theta.size), np.nan) P_mos = P_per.copy() for i in range(theta.size): - P_per[i] = Zo*F_re[i]*re*(lamb**2)*(1. + abs(np.cos(2.*theta[i])))/( + P_per[i] = Zo*F_re[i]*re*((lamb*1e10)**2)*(1. + abs(np.cos(2.*theta[i])))/( 6.*np.pi*Volume[i]*np.sin(2.*theta[i]) ) - P_mos[i] = (F_re[i]**2)*(re**2)*(lamb**3)*( + P_mos[i] = (F_re[i]**2)*(re**2)*((lamb*1e10)**3)*( 1. + (np.cos(2.*theta[i]))**2 )/(4.*mu*(Volume[i]**2)*np.sin(2.*theta[i])) @@ -1292,7 +1312,7 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( # Lattice parameters din=None, lamb=None, # Lattice modifications - alpha_limits=None, + miscut_limits=None, theta=None, theta_deg=None, miscut=None, na=None, nn=None, therm_exp=None, T0=None, TD=None, @@ -1359,7 +1379,7 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( ax = fig1.add_subplot(gs[0, 0]) ax.set_title( f'{name}, ' + f'({miller[0]},{miller[1]},{miller[2]})' + - fr', $\lambda$={lamb} $\AA$', fontsize=15, + fr', $\lambda$={lamb}m', fontsize=15, ) ax.set_xlabel(r'Diffracting angle $\theta$ (rad)', fontsize=15) ax.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=15) @@ -1379,7 +1399,7 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( ax22 = fig1.add_subplot(gs[2, 2]) fig1.suptitle( f'{name}, ' + f'({miller[0]},{miller[1]},{miller[2]})' + - fr', $\lambda$={lamb} $\AA$' + + fr', $\lambda$={lamb}m' + r', $\theta_{B}$=' + fr'{np.round(theta[nn], 5)} rad', fontsize=15, ) @@ -1700,7 +1720,7 @@ def CrystalBragg_plot_rc_components_vs_asymmetry( ] ax.set_title( f'{name}' + f', ({miller[0]},{miller[1]},{miller[2]})' + - fr', $\lambda$={lamb} $\AA$' + fr', $\lambda$={lamb}m' ) ax.set_xlabel(r'$\alpha$ (deg)', fontsize=15) ax.set_ylim(0., 5.) @@ -1787,7 +1807,7 @@ def CrystalBragg_plot_cmaps_rc_components_vs_asymmetry_temp( ] fig.suptitle( f'{name}' + f', ({miller[0]},{miller[1]},{miller[2]})' + - fr', $\lambda$={lamb} $\AA$' + + fr', $\lambda$={lamb}m' + r', $\theta_{B}$=' + fr'{np.round(theta[nn], 5)} rad', fontsize=15, ) diff --git a/tofu/spectro/_rockingcurve_def.py b/tofu/spectro/_rockingcurve_def.py index c846997b8..dc0753a8e 100644 --- a/tofu/spectro/_rockingcurve_def.py +++ b/tofu/spectro/_rockingcurve_def.py @@ -350,10 +350,9 @@ def hexa_volume(aa, cc): return (aa**2) * cc * (np.sqrt(3.)/2.) def hexa_spacing(hh, kk, ll, aa, cc): - return np.sqrt( - (3.*(aa**2)*(cc**2)) - / (4.*(hh**2 + kk**2 + hh*kk)*(cc**2) + 3.*(ll**2)*(aa**2)) - ) + return (3.*aa**2*cc**2)**(1/2) / ( + 4.*cc**2*(hh**2 + kk**2 + hh*kk) + 3.*ll**2*aa**2 + )**(1/2) # ----------- # Attribution @@ -371,32 +370,6 @@ def hexa_spacing(hh, kk, ll, aa, cc): elif "Ge" in ii: None -# --------------------------------------------------------------- -# Attribution to alpha-Quartz crystals: Quartz_110 and Quartz_102 -# --------------------------------------------------------------- - - -# Same values for 110- and Quartz_102 -a = _DCRYST['Quartz_110']['inter_atomic']['distances']['a0'] -c = _DCRYST['Quartz_110']['inter_atomic']['distances']['c0'] - -h110 = _DCRYST['Quartz_110']['miller'][0] -k110 = _DCRYST['Quartz_110']['miller'][1] -l110 = _DCRYST['Quartz_110']['miller'][2] -h102 = _DCRYST['Quartz_102']['miller'][0] -k102 = _DCRYST['Quartz_102']['miller'][1] -l102 = _DCRYST['Quartz_102']['miller'][2] - -_DCRYST['Quartz_110']['volume'] = hexa_volume(a, c) -_DCRYST['Quartz_110']['d_hkl'] = hexa_spacing(h110, k110, l110, a, c) * 1.e-10 -_DCRYST['Quartz_102']['volume'] = hexa_volume(a, c) -_DCRYST['Quartz_102']['d_hkl'] = hexa_spacing(h102, k102, l102, a, c) * 1e-10 - - -# --------------------------------- -# Attribution to Germanium crystals -# --------------------------------- ->>>>>>> devel # ############################################################################# # ############################################################################# From 7610836bd9ffd22b3e91175d8d3cc8124469f4d3 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Tue, 30 Aug 2022 14:37:45 +0200 Subject: [PATCH 05/14] [#673] Removing astropy.units and using scipy.constants into _rockingcurve_def.py --- tofu/spectro/_rockingcurve_def.py | 47 +++++++++---------- tofu/tests/tests00_root/test_03_plot.py | 1 - .../tests04_spectro/test_03_rockingcurve.py | 1 - 3 files changed, 23 insertions(+), 26 deletions(-) diff --git a/tofu/spectro/_rockingcurve_def.py b/tofu/spectro/_rockingcurve_def.py index dc0753a8e..46ac0e98d 100644 --- a/tofu/spectro/_rockingcurve_def.py +++ b/tofu/spectro/_rockingcurve_def.py @@ -2,7 +2,6 @@ import numpy as np import scipy.interpolate import scipy.constants as scpct -import astropy.units as u # ############################################################################# # ############################################################################# @@ -17,7 +16,7 @@ 'target': { 'ion': 'ArXVII', 'wavelength': 3.96e-10, - 'unit': u.m, + 'unit': scpct.unit('classical electron radius'), }, 'atoms': ['Si', 'O'], 'atoms_Z': [14., 8.], @@ -25,11 +24,11 @@ 'miller': np.r_[1., 1., 0.], 'volume': { 'value': None, - 'unit': 1/u.m**3, + 'unit': scpct.unit('classical electron radius')**-3, }, 'd_hkl': { 'value': None, - 'unit': u.m, + 'unit': scpct.unit('classical electron radius'), }, 'mesh': { 'type': 'hexagonal', @@ -60,11 +59,11 @@ 'a0': 4.91304e-10, 'c0': 5.40463e-10, }, - 'unit': u.m, + 'unit': scpct.unit('classical electron radius'), 'comments': 'within the unit cell', 'Tref': { - 'data': 25., - 'unit': u.deg_C, + 'data': 25. + 273.15, + 'unit': scpct.unit('Boltzmann constant')[2], }, 'sources': 'R.W.G. Wyckoff, Crystal Structures', }, @@ -73,7 +72,7 @@ 'alpha_a': 1.337e-5, 'alpha_c': 7.97e-6, }, - 'unit': 1/u.deg_C, + 'unit': scpct.unit('Boltzmann constant')[2]**-1, 'comments': 'in parallel directions to a0 and c0', 'sources': 'R.W.G. Wyckoff, Crystal Structures', }, @@ -87,7 +86,7 @@ 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, ]*1e10, }, - 'unit': 1/u.m, + 'unit': scpct.unit('classical electron radius')**-1, 'sources': 'Int. Tab. X-Ray Crystallography, Vol.I,II,III,IV (1985)', }, @@ -112,7 +111,7 @@ 'target': { 'ion': 'ArXVIII', 'wavelength': 3.75e-10, - 'unit': u.m, + 'unit': scpct.unit('classical electron radius'), }, 'atoms': ['Si', 'O'], 'atoms_Z': [14., 8.], @@ -120,11 +119,11 @@ 'miller': np.r_[1., 0., 2.], 'volume': { 'value': None, - 'unit': 1/u.m**3, + 'unit': scpct.unit('classical electron radius')**-3, }, 'd_hkl': { 'value': None, - 'unit': u.m, + 'unit': scpct.unit('classical electron radius'), }, 'mesh': { 'type': 'hexagonal', @@ -155,11 +154,11 @@ 'a0': 4.91304e-10, 'c0': 5.40463e-10, }, - 'unit': u.m, + 'unit': scpct.unit('classical electron radius'), 'comments': 'within the unit cell', 'Tref': { - 'data': 25., - 'unit': u.deg_C, + 'data': 25. + 273.15, + 'unit': scpct.unit('Boltzmann constant')[2], }, 'sources': 'R.W.G. Wyckoff, Crystal Structures', }, @@ -168,7 +167,7 @@ 'alpha_a': 1.337e-5, 'alpha_c': 7.97e-6, }, - 'unit': 1/u.deg_C, + 'unit': scpct.unit('Boltzmann constant')[2]**-1, 'comments': 'in parallel directions to a0 and c0', 'sources': 'R.W.G. Wyckoff, Crystal Structures', }, @@ -182,7 +181,7 @@ 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, ]*1e10, }, - 'unit': 1/u.m, + 'unit': scpct.unit('classical electron radius')**-1, 'sources': 'Int. Tab. X-Ray Crystallography, Vol.I,II,III,IV (1985)', }, @@ -208,7 +207,7 @@ 'target': { 'ion': None, 'wavelength': None, - 'unit': u.m, + 'unit': scpct.unit('classical electron radius'), }, 'atoms': None, 'atoms_Z': None, @@ -216,11 +215,11 @@ 'miller': None, 'volume': { 'value': None, - 'unit': 1/u.m**3, + 'unit': scpct.unit('classical electron radius')**-3, }, 'd_hkl': { 'value': None, - 'unit': u.m, + 'unit': scpct.unit('classical electron radius'), }, 'mesh': { 'type': None, @@ -230,23 +229,23 @@ 'phases': None, 'inter_atomic': { 'distances': None, - 'unit': u.m, + 'unit': scpct.unit('classical electron radius'), 'comments': None, 'Tref': { 'data': None, - 'unit': u.deg_C, + 'unit': scpct.unit('Boltzmann constant')[2], }, 'sources': None, }, 'thermal_expansion': { 'coefs': None, - 'unit': 1/u.deg_C, + 'unit': scpct.unit('Boltzmann constant')[2]**-1, 'comments': None, 'sources': None, }, 'sin_theta_lambda': { 'values': None, - 'unit': 1/u.m, + 'unit': scpct.unit('classical electron radius')**-1, 'sources': None, }, 'atomic_scattering': { diff --git a/tofu/tests/tests00_root/test_03_plot.py b/tofu/tests/tests00_root/test_03_plot.py index fb7a0023c..9c59ac8ef 100644 --- a/tofu/tests/tests00_root/test_03_plot.py +++ b/tofu/tests/tests00_root/test_03_plot.py @@ -9,7 +9,6 @@ import numpy as np import matplotlib.pyplot as plt import warnings as warn -import astropy.units as u # Importing package tofu.geom import tofu as tf diff --git a/tofu/tests/tests04_spectro/test_03_rockingcurve.py b/tofu/tests/tests04_spectro/test_03_rockingcurve.py index d5b07d0e2..f64ebcc86 100644 --- a/tofu/tests/tests04_spectro/test_03_rockingcurve.py +++ b/tofu/tests/tests04_spectro/test_03_rockingcurve.py @@ -11,7 +11,6 @@ # Standard import numpy as np import scipy.constants as scpct -import astropy.units as u import matplotlib.pyplot as plt # tofu-specific From cd1fb89708dbc9d3b96e7d31351973705150cb01 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Tue, 30 Aug 2022 14:55:30 +0200 Subject: [PATCH 06/14] [#673] Good writing of the units --- tofu/spectro/_rockingcurve_def.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/tofu/spectro/_rockingcurve_def.py b/tofu/spectro/_rockingcurve_def.py index 46ac0e98d..aae0baf0b 100644 --- a/tofu/spectro/_rockingcurve_def.py +++ b/tofu/spectro/_rockingcurve_def.py @@ -24,7 +24,7 @@ 'miller': np.r_[1., 1., 0.], 'volume': { 'value': None, - 'unit': scpct.unit('classical electron radius')**-3, + 'unit': scpct.unit('classical electron radius')+'**-3', }, 'd_hkl': { 'value': None, @@ -72,7 +72,7 @@ 'alpha_a': 1.337e-5, 'alpha_c': 7.97e-6, }, - 'unit': scpct.unit('Boltzmann constant')[2]**-1, + 'unit': '1/'+scpct.unit('Boltzmann constant')[2], 'comments': 'in parallel directions to a0 and c0', 'sources': 'R.W.G. Wyckoff, Crystal Structures', }, @@ -86,7 +86,7 @@ 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, ]*1e10, }, - 'unit': scpct.unit('classical electron radius')**-1, + 'unit': '1/'+scpct.unit('classical electron radius'), 'sources': 'Int. Tab. X-Ray Crystallography, Vol.I,II,III,IV (1985)', }, @@ -119,7 +119,7 @@ 'miller': np.r_[1., 0., 2.], 'volume': { 'value': None, - 'unit': scpct.unit('classical electron radius')**-3, + 'unit': scpct.unit('classical electron radius')+'**-3', }, 'd_hkl': { 'value': None, @@ -167,7 +167,7 @@ 'alpha_a': 1.337e-5, 'alpha_c': 7.97e-6, }, - 'unit': scpct.unit('Boltzmann constant')[2]**-1, + 'unit': '1/'+scpct.unit('Boltzmann constant')[2], 'comments': 'in parallel directions to a0 and c0', 'sources': 'R.W.G. Wyckoff, Crystal Structures', }, @@ -181,7 +181,7 @@ 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, ]*1e10, }, - 'unit': scpct.unit('classical electron radius')**-1, + 'unit': '1/'+scpct.unit('classical electron radius'), 'sources': 'Int. Tab. X-Ray Crystallography, Vol.I,II,III,IV (1985)', }, @@ -215,7 +215,7 @@ 'miller': None, 'volume': { 'value': None, - 'unit': scpct.unit('classical electron radius')**-3, + 'unit': scpct.unit('classical electron radius')+'**-3', }, 'd_hkl': { 'value': None, @@ -239,13 +239,13 @@ }, 'thermal_expansion': { 'coefs': None, - 'unit': scpct.unit('Boltzmann constant')[2]**-1, + 'unit': '1/'+scpct.unit('Boltzmann constant')[2], 'comments': None, 'sources': None, }, 'sin_theta_lambda': { 'values': None, - 'unit': scpct.unit('classical electron radius')**-1, + 'unit': '1/'+scpct.unit('classical electron radius'), 'sources': None, }, 'atomic_scattering': { From 2b1ad372865a7168444a96cc5ce3092af8839133 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Tue, 30 Aug 2022 15:47:29 +0200 Subject: [PATCH 07/14] [#673] Renamed Polygon module --- tofu/geom/_etendue.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tofu/geom/_etendue.py b/tofu/geom/_etendue.py index d8dfa6c84..722a1c6a5 100644 --- a/tofu/geom/_etendue.py +++ b/tofu/geom/_etendue.py @@ -9,7 +9,7 @@ import matplotlib.lines as mlines -import polygon as plg +import Polygon as plg import datastock as ds From 93dc64750199b44416bf930e40b2950ee3a529c2 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 31 Aug 2022 12:07:39 +0200 Subject: [PATCH 08/14] [#673] Trying to fix errors with tests --- tofu/geom/_check_optics.py | 25 ++++++++++++++----------- tofu/geom/_core_optics.py | 34 +++++++++++++++++++++++++++------- tofu/geom/_etendue.py | 2 +- 3 files changed, 42 insertions(+), 19 deletions(-) diff --git a/tofu/geom/_check_optics.py b/tofu/geom/_check_optics.py index 5a1cf7a9f..07805ece6 100644 --- a/tofu/geom/_check_optics.py +++ b/tofu/geom/_check_optics.py @@ -515,26 +515,29 @@ def _checkformat_dbragg(dbragg=None, ddef=None, valid_keys=None, dmat=None): if dbragg.get('lambref') is None: # set to bragg = braggref dbragg['lambref'] = _comp_optics.get_lamb_from_bragg( - np.r_[dbragg['braggref']], - dmat['d'], + bragg=np.r_[dbragg['braggref']], + d=dmat['d'], n=1, )[0] user_prov = False # Set default bragg angle if necessary braggref = _comp_optics.get_bragg_from_lamb( - np.r_[dbragg['lambref']], dmat['d'], n=1)[0] + lamb=np.r_[dbragg['lambref']], + d=dmat['d'], + n=1, + )[0] if np.isnan(braggref): lambok = [] msg = ( - """ - Var {} is not valid! - Please check your arguments to calculate the Bragg's law correctly! - Provided: - - crystal inter-plane d [m] = {} - - wavelenght interval [m] : {} - - lambref = {} - """.format('lambref', dmat['d'], lambok, dbragg['lambref']) + "Var {} or {} is not valid!\n".format('lambref', 'braggref') + + "Please check your args to compute the Bragg's law correctly!\n" + + "Provided:\n" + + "\t - crystal inter-plane d [m] = {}\n".format(dmat['d']) + # + "\t - wavelenght [m] = {}\n".format(lambok) + + "\t - lambref = {}\n".format(dbragg['lambref']) + + "\t - Bragg angle [m] = {}\n".format(braggref) + + "\t - braggref = {}\n".format(dbragg['braggref']) ) raise Exception(msg) diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index f9ab42ccf..00d053ebe 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -1690,9 +1690,20 @@ def get_detector_ideal( rcurve = self._dgeom['rcurve'] bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb, n=n) + lamb = self.get_lamb_from_bragg(bragg=bragg, n=n) + if np.all(np.isnan(bragg)): - msg = ("There is no available bragg angle!\n" - + " => Check the vlue of self.dmat['d'] vs lamb") + msg = ( + "There is no available bragg angle!\n" + + "It is probably because you gave a too large Bragg angle\n" + + "or wavelength than permitted by the crystal loaded.\n" + + "Provided:\n" + + "\t - bragg : {} rad\n".format(bragg) + + "\t - lamb : {} m\n".format(lamb) + + "Crystal loaded:\n" + + "\t - bragg_ref : {} rad\n".format(self._dbragg['braggref']) + + "\t - lamb_ref : {} m\n".format(self._dbragg['lambref']) + ) raise Exception(msg) lc = [lamb0 is not None, lamb1 is not None, dist01 is not None] @@ -2085,7 +2096,7 @@ def plot_line_on_det_tracing( dout = _rockingcurve.CrystBragg_comp_lattice_spacing( crystal=crystal, din=din, - lamb=self.dbragg['lambref'], + lamb=self._dbragg['lambref'], na=na, nn=nn, therm_exp=therm_exp, temp_limits=temp_limits, @@ -2723,6 +2734,10 @@ def plot_focal_error_summed( ndist = 21 if ndi is None: ndi = 21 + if lamb is None: + lamb = self._dbragg['lambref'] + if bragg is None: + bragg = self._dbragg['braggref'] if err is None: err = 'rel' if plot is None: @@ -2897,9 +2912,11 @@ def _get_local_coordinates_of_det( # Checkformat det det_ref = self._checkformat_det(det=det_ref) + # Checkformat lamb + lamb = self.get_lamb_from_bragg(bragg=self._dbragg['braggref']) + # ------------ # get approx detect - det_approx = self.get_detector_ideal( bragg=bragg, lamb=lamb, tangent_to_rowland=False, @@ -2948,8 +2965,8 @@ def get_lambbraggphi_from_ptsxixj_dthetapsi( dtheta=None, psi=None, ntheta=None, npsi=None, n=None, - miscut=None, grid=None, + miscut=None, return_lamb=None, ): """ Return the lamb, bragg and phi for provided pts and dtheta/psi @@ -2995,9 +3012,11 @@ def get_lambbraggphi_from_ptsxixj_dthetapsi( def get_lamb_avail_from_pts( self, pts=None, - n=None, ndtheta=None, - det=None, nlamb=None, klamb=None, + lamb=None, nlamb=None, miscut=None, + n=None, ndtheta=None, + det=None, + klamb=None, strict=None, return_phidtheta=None, return_xixj=None, @@ -3437,6 +3456,7 @@ def get_plasmadomain_at_lamb( lamb_access = self.get_lamb_avail_from_pts( pts=ptsXYZ, + lamb=lamb, nlamb=2, miscut=miscut, return_phidtheta=False, diff --git a/tofu/geom/_etendue.py b/tofu/geom/_etendue.py index 722a1c6a5..d8dfa6c84 100644 --- a/tofu/geom/_etendue.py +++ b/tofu/geom/_etendue.py @@ -9,7 +9,7 @@ import matplotlib.lines as mlines -import Polygon as plg +import polygon as plg import datastock as ds From 330f4516c82718f2bba87af8db8fc7df2bda8f77 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 31 Aug 2022 17:33:59 +0200 Subject: [PATCH 09/14] [#673] Rectified Polygon module import --- tofu/geom/_core_optics.py | 10 ++-------- tofu/geom/_etendue.py | 2 +- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index 00d053ebe..802dbc922 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -1690,7 +1690,6 @@ def get_detector_ideal( rcurve = self._dgeom['rcurve'] bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb, n=n) - lamb = self.get_lamb_from_bragg(bragg=bragg, n=n) if np.all(np.isnan(bragg)): msg = ( @@ -2734,10 +2733,6 @@ def plot_focal_error_summed( ndist = 21 if ndi is None: ndi = 21 - if lamb is None: - lamb = self._dbragg['lambref'] - if bragg is None: - bragg = self._dbragg['braggref'] if err is None: err = 'rel' if plot is None: @@ -2753,6 +2748,8 @@ def plot_focal_error_summed( if lambda_interval_max is None: lambda_interval_max = 4.00e-10 + bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb) + l0 = [dist_min, dist_max, ndist, di_min, di_max, ndi] c0 = any([l00 is not None for l00 in l0]) if not c0: @@ -2912,9 +2909,6 @@ def _get_local_coordinates_of_det( # Checkformat det det_ref = self._checkformat_det(det=det_ref) - # Checkformat lamb - lamb = self.get_lamb_from_bragg(bragg=self._dbragg['braggref']) - # ------------ # get approx detect det_approx = self.get_detector_ideal( diff --git a/tofu/geom/_etendue.py b/tofu/geom/_etendue.py index d8dfa6c84..722a1c6a5 100644 --- a/tofu/geom/_etendue.py +++ b/tofu/geom/_etendue.py @@ -9,7 +9,7 @@ import matplotlib.lines as mlines -import polygon as plg +import Polygon as plg import datastock as ds From 0bd4973a071cbb1c27fcbfb92eddecd8d1aaf453 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Mon, 12 Sep 2022 15:03:08 +0200 Subject: [PATCH 10/14] [#673] Possibility to plot few RCs on the same window --- tofu/geom/_core_optics.py | 28 ++-- tofu/spectro/_rockingcurve.py | 248 +++++++++++++++++++--------------- 2 files changed, 156 insertions(+), 120 deletions(-) diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index 802dbc922..04bcec61d 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -812,7 +812,7 @@ def compute_rockingcurve( miscut_limits=None, therm_exp=None, temp_limits=None, - plot_therm_exp=None, + ax=None, plot_therm_exp=None, plot_asf=None, plot_power_ratio=None, plot_asymmetry=None, plot_cmaps=None, returnas=None, @@ -824,7 +824,7 @@ def compute_rockingcurve( miscut_limits=miscut_limits, therm_exp=therm_exp, temp_limits=temp_limits, - plot_therm_exp=plot_therm_exp, + ax=ax, plot_therm_exp=plot_therm_exp, plot_asf=plot_asf, plot_power_ratio=plot_power_ratio, plot_asymmetry=plot_asymmetry, plot_cmaps=plot_cmaps, returnas=None, @@ -1602,7 +1602,7 @@ def calc_meridional_sagittal_focus( alpha = alpha if miscut is None or miscut is False: miscut = False - alpha = self.dmat['alpha'] + alpha = self._dmat['alpha'] # Compute return _comp_optics.calc_meridional_sagittal_focus( @@ -1618,7 +1618,7 @@ def get_rowland_dist_from_lambbragg(self, bragg=None, lamb=None, n=None): bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb, n=n) if np.all(np.isnan(bragg)): msg = ("There is no available bragg angle!\n" - + " => Check the vlue of self.dmat['d'] vs lamb") + + " => Check the value of self._dmat['d'] vs lamb") raise Exception(msg) return _comp_optics.get_rowland_dist_from_bragg( bragg=bragg, rcurve=self._dgeom['rcurve'], @@ -1688,9 +1688,10 @@ def get_detector_ideal( if rcurve is None: rcurve = self._dgeom['rcurve'] + if lamb is None: + lamb = self._dbragg['lambref'] bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb, n=n) - if np.all(np.isnan(bragg)): msg = ( "There is no available bragg angle!\n" @@ -2115,7 +2116,7 @@ def find_nearest(array, value): return idx id_temp0 = find_nearest(TD, temp0) - self.dmat['d'] = d_atom[id_temp0] + self._dmat['d'] = d_atom[id_temp0] # Get local basis nout, e1, e2, miscut = self.get_unit_vectors( @@ -2226,7 +2227,7 @@ def find_nearest(array, value): # For each wavelength, get results dictionnary of the associated # diffraction pattern for ll in range(nlamb): - dout = _rockingcurve.compute_rockingcurve( + dout, _ = _rockingcurve.compute_rockingcurve( crystal=crystal, din=din, lamb=lamb[ll], miscut=miscut, @@ -2340,9 +2341,9 @@ def find_nearest(array, value): self.update_miscut(alpha=0., beta=0.) if therm_exp: - self.dmat['d'] = d_atom[nn]*1e-10 + self._dmat['d'] = d_atom[nn]*1e-10 else: - self.dmat['d'] = d_atom[0]*1e-10 + self._dmat['d'] = d_atom[0]*1e-10 ( bragg_atprmax[ll], _, lamb_atprmax[ll], ) = self.get_lambbraggphi_from_ptsxixj_dthetapsi( @@ -2360,9 +2361,9 @@ def find_nearest(array, value): else: self.update_miscut(alpha=0., beta=0.) if therm_exp: - self.dmat['d'] = d_atom[id_temp0]*1e-10 + self._dmat['d'] = d_atom[id_temp0]*1e-10 else: - self.dmat['d'] = d_atom[0]*1e-10 + self._dmat['d'] = d_atom[0]*1e-10 # Plot if plot: @@ -2748,6 +2749,8 @@ def plot_focal_error_summed( if lambda_interval_max is None: lambda_interval_max = 4.00e-10 + # If lamb & bragg are None, the following routine returns + # bragg=self._dbragg['braggref'] bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb) l0 = [dist_min, dist_max, ndist, di_min, di_max, ndi] @@ -2909,6 +2912,9 @@ def _get_local_coordinates_of_det( # Checkformat det det_ref = self._checkformat_det(det=det_ref) + # Checkformat bragg & lamb + bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb) + # ------------ # get approx detect det_approx = self.get_detector_ideal( diff --git a/tofu/spectro/_rockingcurve.py b/tofu/spectro/_rockingcurve.py index 76c91d89b..c4b684efe 100644 --- a/tofu/spectro/_rockingcurve.py +++ b/tofu/spectro/_rockingcurve.py @@ -36,6 +36,7 @@ def compute_rockingcurve( therm_exp=None, temp_limits=None, # Plot + ax=None, plot_therm_exp=None, plot_asf=None, plot_power_ratio=None, @@ -376,7 +377,7 @@ def compute_rockingcurve( # ---------------- if plot_power_ratio: - CrystalBragg_plot_power_ratio_vs_glancing_angle( + ax = CrystalBragg_plot_power_ratio_vs_glancing_angle( din=din, lamb=lamb, miscut_limits=miscut_limits, theta=theta, theta_deg=theta_deg, @@ -384,6 +385,7 @@ def compute_rockingcurve( bb=bb, polar=polar, alpha=alpha, miscut=miscut, na=na, nn=nn, therm_exp=therm_exp, T0=T0, TD=TD, + ax=ax, ) # Plot integrated reflect., asymmetry angle & RC width vs glancing angle @@ -455,7 +457,7 @@ def compute_rockingcurve( dout['Temperature changes (°C)'] = TD if returnas is dict: - return dout + return dout, ax # ############################################################################# @@ -1244,15 +1246,15 @@ def CrystalBragg_plot_thermal_expansion_vs_d( T0=None, TD=None, d_atom=None, nn=None, ): - fig = plt.figure(figsize=(9, 6)) - gs = gridspec.GridSpec(1, 1) - ax = fig.add_subplot(gs[0, 0]) name = din['name'] miller = np.r_[ int(din['miller'][0]), int(din['miller'][1]), int(din['miller'][2]), ] + fig = plt.figure(figsize=(9, 6)) + gs = gridspec.GridSpec(1, 1) + ax = fig.add_subplot(gs[0, 0]) ax.set_title( f'{name}' + f', ({miller[0]},{miller[1]},{miller[2]})' + fr', $\lambda$={lamb} $\AA$' + @@ -1260,25 +1262,26 @@ def CrystalBragg_plot_thermal_expansion_vs_d( fontsize=15, ) ax.set_xlabel(r'$\Delta$T ($T_{0}$'+fr'={T0}°C)', fontsize=15) - ax.set_ylabel(r'Inter-planar distance $d_{hkl}$ [m$\AA$]', fontsize=15) + ax.set_ylabel(r'Inter-planar distance $d_{hkl}$ [$\AA$]', fontsize=15) ax.scatter( - TD, d_atom*(1e3), + TD, d_atom, marker='o', c='k', alpha=0.5, label=r'd$_{(hkl)}$ computed points', ) - p = np.polyfit(TD, d_atom*(1e3), 1) - y_adj = p[0]*TD + p[1] + interp_d = np.polyfit(TD, d_atom, 1) + y_adj = interp_d[0]*TD + interp_d[1] ax.plot( TD, y_adj, 'k-', label=( - r'$d_{hkl}$ = ' + str(np.round(p[1], 3)) + - r' x (1 + $\gamma_{eff}$.$\Delta$T)' + '\n' + - r'$\gamma_{eff}$ = ' + str(np.round(p[0]/p[1], decimals=9)) + - r'°C$^{-1}$', + '$d_{hkl}$ = ' + str(np.round(interp_d[1], 13)) + + ' x (1 + ' + r'$\gamma_{eff}.\Delta$T)' + '\n' + + r'$\gamma_{eff}$ = ' + + str(np.round(interp_d[0]/interp_d[1], decimals=9)) + + r'°C$^{-1}$' ), ) - ax.legend(loc="best", fontsize=12) + ax.legend(fontsize=12) def CrystalBragg_plot_atomic_scattering_factor( @@ -1319,6 +1322,7 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( # Diffraction pattern main components th=None, dth=None, power_ratio=None, bb=None, polar=None, alpha=None, + ax=None, ): """ All plots of rocking curve is done, not with respect to the glancing angle (theta - thetaBragg) where thetaBragg may vary if the temperature @@ -1373,66 +1377,70 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( miscut is False and therm_exp is True, miscut is True and therm_exp is False, ] - if any(lc): - fig1 = plt.figure(figsize=(8, 6)) - gs = gridspec.GridSpec(1, 1) - ax = fig1.add_subplot(gs[0, 0]) - ax.set_title( - f'{name}, ' + f'({miller[0]},{miller[1]},{miller[2]})' + - fr', $\lambda$={lamb}m', fontsize=15, - ) - ax.set_xlabel(r'Diffracting angle $\theta$ (rad)', fontsize=15) - ax.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=15) - if miscut and therm_exp: - gs = gridspec.GridSpec(3, 3) - fig1 = plt.figure(figsize=(22, 20)) - # 3 rows -> temperature changes -T0 < 0 < +T0 - # 3 columns -> asymmetry angles -bragg < 0 < +bragg - ax01 = fig1.add_subplot(gs[0, 1]) - ax00 = fig1.add_subplot(gs[0, 0]) - ax02 = fig1.add_subplot(gs[0, 2]) - ax11 = fig1.add_subplot(gs[1, 1]) - ax10 = fig1.add_subplot(gs[1, 0]) - ax12 = fig1.add_subplot(gs[1, 2]) - ax21 = fig1.add_subplot(gs[2, 1]) - ax20 = fig1.add_subplot(gs[2, 0]) - ax22 = fig1.add_subplot(gs[2, 2]) - fig1.suptitle( - f'{name}, ' + f'({miller[0]},{miller[1]},{miller[2]})' + - fr', $\lambda$={lamb}m' + - r', $\theta_{B}$=' + fr'{np.round(theta[nn], 5)} rad', - fontsize=15, - ) - ax00.set_title( - r'$\alpha$=({}) arcmin'.format(np.round(alpha_deg[0]*60, 3)), - fontsize=15 - ) - ax01.set_title( - r'$\alpha$=({}) arcmin'.format(np.round(alpha_deg[nn]*60, 3)), - fontsize=15 - ) - ax02.set_title( - r'$\alpha$=({}) arcmin'.format(np.round(alpha_deg[na-1]*60, 3)), - fontsize=15 - ) - ax022 = ax02.twinx() - ax022.set_ylabel( - r'$\Delta T$=({})°C'.format(TD[0]), fontsize=15 - ) - ax122 = ax12.twinx() - ax122.set_ylabel( - r'$\Delta T$=({})°C'.format(TD[nn]), fontsize=15 - ) - ax222 = ax22.twinx() - ax222.set_ylabel( - r'$\Delta T$=({})°C'.format(TD[na-1]), fontsize=15 - ) - ax20.set_xlabel(r'$\theta$ (rad)', fontsize=15) - ax21.set_xlabel(r'$\theta$ (rad)', fontsize=15) - ax22.set_xlabel(r'$\theta$ (rad)', fontsize=15) - ax00.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=15) - ax10.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=15) - ax20.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=15) + if ax is None: + if any(lc): + fig1 = plt.figure(figsize=(8, 6)) + gs = gridspec.GridSpec(1, 1) + ax = fig1.add_subplot(gs[0, 0]) + """ + ax.set_title( + f'{name}, ' + f'({miller[0]},{miller[1]},{miller[2]})' + + fr', $\lambda$={lamb}m', fontsize=15, + ) + """ + ax.set_title(f'{name}') + ax.set_xlabel(r'Diffracting angle $\theta$ (rad)', fontsize=15) + ax.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=15) + elif miscut and therm_exp: + gs = gridspec.GridSpec(3, 3) + fig1 = plt.figure(figsize=(22, 20)) + # 3 rows -> temperature changes -T0 < 0 < +T0 + # 3 columns -> asymmetry angles -bragg < 0 < +bragg + ax01 = fig1.add_subplot(gs[0, 1]) + ax00 = fig1.add_subplot(gs[0, 0]) + ax02 = fig1.add_subplot(gs[0, 2]) + ax11 = fig1.add_subplot(gs[1, 1]) + ax10 = fig1.add_subplot(gs[1, 0]) + ax12 = fig1.add_subplot(gs[1, 2]) + ax21 = fig1.add_subplot(gs[2, 1]) + ax20 = fig1.add_subplot(gs[2, 0]) + ax22 = fig1.add_subplot(gs[2, 2]) + fig1.suptitle( + f'{name}, ' + f'({miller[0]},{miller[1]},{miller[2]})' + + fr', $\lambda$={lamb}m' + + r', $\theta_{B}$=' + fr'{np.round(theta[nn], 5)} rad', + fontsize=15, + ) + ax00.set_title( + r'$\alpha$=({}) arcmin'.format(np.round(alpha_deg[0]*60, 3)), + fontsize=15 + ) + ax01.set_title( + r'$\alpha$=({}) arcmin'.format(np.round(alpha_deg[nn]*60, 3)), + fontsize=15 + ) + ax02.set_title( + r'$\alpha$=({}) arcmin'.format(np.round(alpha_deg[na-1]*60, 3)), + fontsize=15 + ) + ax022 = ax02.twinx() + ax022.set_ylabel( + r'$\Delta T$=({})°C'.format(TD[0]), fontsize=15 + ) + ax122 = ax12.twinx() + ax122.set_ylabel( + r'$\Delta T$=({})°C'.format(TD[nn]), fontsize=15 + ) + ax222 = ax22.twinx() + ax222.set_ylabel( + r'$\Delta T$=({})°C'.format(TD[na-1]), fontsize=15 + ) + ax20.set_xlabel(r'$\theta$ (rad)', fontsize=15) + ax21.set_xlabel(r'$\theta$ (rad)', fontsize=15) + ax22.set_xlabel(r'$\theta$ (rad)', fontsize=15) + ax00.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=15) + ax10.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=15) + ax20.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=15) # Plot # ---- @@ -1454,18 +1462,19 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( c_keydd = c_keylist[c_valuedd] ax.text( dth[0, 0, j, ind], - np.max(power_ratio[0, 0, j, :] + 0.005), + np.max(power_ratio[0, 0, j] + power_ratio[1, 0, j] + 0.005), '({})'.format(let_keydd), c=c_keydd, ) + # Plot each polarization + """ ax.plot( dth[0, 0, j, :], power_ratio[0, 0, j, :], - '-', + linestyle='-', c=c_keydd, label=( - r'normal pola.,' + '\n' + - r' ({}): $\alpha$=({}) arcmin'.format( + r'normal pola., ({}): $\alpha$=({})arcmin'.format( let_keydd, np.round(alpha_deg[j]*60, 3) ), ), @@ -1473,34 +1482,34 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( ax.plot( dth[1, 0, j, :], power_ratio[1, 0, j, :], - '--', + linestyle='--', + c=c_keydd, + label=( + r'parallel pola., ({}): $\lambda$=({})m'.format( + let_keydd, lamb, + ), + ), + ) + """ + # Plot the sum of both polarizations + ax.plot( + dth[0, 0, j, :], + power_ratio[0, 0, j] + power_ratio[1, 0, j], + '-', c=c_keydd, - label=r'parallel pola.', + label=( + '(' + str(let_keydd) + r'): $\alpha$=' + + str(np.round(alpha_deg[j]*60, 3)) + ' arcmin' + '\n' + r'$\lambda$=' + str(lamb) + ' m' + ), ) ax.axvline( - theta, color='black', linestyle='-.', - label=r'$\theta_B$= {} rad'.format( - np.round(theta, 6) - ), + theta, + color='black', + linestyle='-.', + label=r'$\theta_B$={}rad'.format(np.round(theta, 6)), ) ax.legend(fontsize=12) - """ - # Plot the sum of both polarizations - lc = [miscut is True, miscut is False] - if not therm_exp and any(lc): - ax.plot( - dth[0, 0, 0, :], - power_ratio[0, 0, 0] + power_ratio[1, 0, 0], - '-', - c='black', - ) - ax.axvline( - theta, color='black', linestyle='-.', - label=r'$\theta_B$= {} rad'.format( - np.round(theta, 6) - ), - ) - """ if not miscut and therm_exp is True: colors = ['blue', 'black', 'red'] @@ -1519,14 +1528,16 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( c_keydd = c_keylist[c_valuedd] ax.text( dth[0, i, 0, ind], - np.max(power_ratio[0, i, 0, :] + 0.005), + np.max(power_ratio[0, i, 0] + power_ratio[1, i, 0] + 0.005), '({})'.format(let_keydd), c=c_keydd, ) + # Plot each polarization + """ ax.plot( dth[0, i, 0, :], power_ratio[0, i, 0, :], - '-', + linestyle='-', c=c_keydd, label=r'normal pola., ({}): $\Delta T$=({})°C'.format( let_keydd, TD[i] @@ -1535,15 +1546,32 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( ax.plot( dth[1, i, 0, :], power_ratio[1, i, 0, :], - '--', + linestyle='--', c=c_keydd, - label=r'parallel pola.', + label=( + r'parallel pola., ({}): $\lambda$=({})m'.format( + let_keydd, lamb, + ), + ), + ) + """ + # Plot the sum of both polarizations + ax.plot( + dth[0, i, 0, :], + power_ratio[0, i, 0] + power_ratio[1, i, 0], + '-', + c=c_keydd, + label=( + '(' + str(let_keydd) + r'): $\Delta$T=' + + str(TD[i]) + ' °C' + '\n' + r'$\lambda$=' + str(lamb) + ' m' + ), ) ax.axvline( - theta[nn], color='black', linestyle='--', - label=r'Bragg angle of ref. : {} rad'.format( - np.round(theta[nn], 6) - ), + theta[nn], + color='black', + linestyle='--', + label=r'$\theta_B$={}rad'.format(np.round(theta[nn], 6)), ) ax.legend(fontsize=12) @@ -1696,6 +1724,8 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( ax11.legend() + return ax + def CrystalBragg_plot_rc_components_vs_asymmetry( din=None, lamb=None, From 1279a89a12a3e97f1881536754da0dc74453761a Mon Sep 17 00:00:00 2001 From: AD265925 Date: Mon, 10 Oct 2022 16:04:39 +0200 Subject: [PATCH 11/14] [#673] Error source on test_04_core_optics found --- tofu/spectro/_rockingcurve.py | 2 +- tofu/tests/tests01_geom/test_04_core_optics.py | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/tofu/spectro/_rockingcurve.py b/tofu/spectro/_rockingcurve.py index c4b684efe..8c9f0815d 100644 --- a/tofu/spectro/_rockingcurve.py +++ b/tofu/spectro/_rockingcurve.py @@ -434,7 +434,7 @@ def compute_rockingcurve( 'Wavelength (m)': lamb, 'Miller indices': (ih, ik, il), 'Inter-reticular distance (m)': d_atom, - 'Volume (A^3)': Volume, + 'Volume (m^3)': Volume, 'Bragg angle of reference (rad)': theta, 'Glancing angles': dth, 'Power ratio': power_ratio, diff --git a/tofu/tests/tests01_geom/test_04_core_optics.py b/tofu/tests/tests01_geom/test_04_core_optics.py index 98ca044f5..0ab3e3b39 100644 --- a/tofu/tests/tests01_geom/test_04_core_optics.py +++ b/tofu/tests/tests01_geom/test_04_core_optics.py @@ -132,6 +132,7 @@ def setup_class(cls, verb=False): } dbragg = { 'lambref': 3.96e-10, + 'braggref': 0.9373938265780166, } dmat1 = dict(dmat) @@ -324,6 +325,7 @@ def test11_calc_johann_error(self): det=det, ) + """ def test12_plot_line_on_det_tracing(self): for k0, obj in self.dobj.items(): det = obj.get_detector_ideal() @@ -339,6 +341,8 @@ def test12_plot_line_on_det_tracing(self): therm_exp=False, plot=False, ) + plt.close('all') + """ def test13_calc_meridional_sagittal_focus(self): derr = {} From 2686bda16112115f5d05ac4ba2c8be41846d82ef Mon Sep 17 00:00:00 2001 From: AD265925 Date: Mon, 10 Oct 2022 18:39:09 +0200 Subject: [PATCH 12/14] [#673] Adjusting line position to pass the tests --- tofu/geom/_core_optics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index 04bcec61d..ef5b72513 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -1688,8 +1688,6 @@ def get_detector_ideal( if rcurve is None: rcurve = self._dgeom['rcurve'] - if lamb is None: - lamb = self._dbragg['lambref'] bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb, n=n) if np.all(np.isnan(bragg)): @@ -1705,6 +1703,8 @@ def get_detector_ideal( + "\t - lamb_ref : {} m\n".format(self._dbragg['lambref']) ) raise Exception(msg) + if lamb is None: + lamb = self._dbragg['lambref'] lc = [lamb0 is not None, lamb1 is not None, dist01 is not None] if any(lc) and not all(lc): From 2b59ae2c5387e895986d8a029f9386e89f6ad215 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Mon, 10 Oct 2022 18:54:51 +0200 Subject: [PATCH 13/14] [#673] Changing use_non_parallelism by miscut in tests04_spectro/test_03_rockingcurve.py --- tofu/tests/tests04_spectro/test_03_rockingcurve.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tofu/tests/tests04_spectro/test_03_rockingcurve.py b/tofu/tests/tests04_spectro/test_03_rockingcurve.py index f64ebcc86..aebc1b1dd 100644 --- a/tofu/tests/tests04_spectro/test_03_rockingcurve.py +++ b/tofu/tests/tests04_spectro/test_03_rockingcurve.py @@ -113,7 +113,7 @@ def test01_rockingcurve_compute(self): dout = tfs.compute_rockingcurve( crystal=k0, lamb=np.r_[3.969067e-10], - use_non_parallelism=False, + miscut=False, alpha_limits=None, therm_exp=False, temp_limits=None, From 21c9136c993e89397325037ed22edecedcce65d5 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Mon, 10 Oct 2022 19:15:27 +0200 Subject: [PATCH 14/14] [#673] Replace alpha_limits by miscut_limits in tests04_spectro/test_03_rockingcurve.py --- tofu/tests/tests04_spectro/test_03_rockingcurve.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tofu/tests/tests04_spectro/test_03_rockingcurve.py b/tofu/tests/tests04_spectro/test_03_rockingcurve.py index aebc1b1dd..1a4c97adf 100644 --- a/tofu/tests/tests04_spectro/test_03_rockingcurve.py +++ b/tofu/tests/tests04_spectro/test_03_rockingcurve.py @@ -114,7 +114,7 @@ def test01_rockingcurve_compute(self): crystal=k0, lamb=np.r_[3.969067e-10], miscut=False, - alpha_limits=None, + miscut_limits=None, therm_exp=False, temp_limits=None, # Plot