From dc3387f4f661e180481f5d8c61e3d578395659ae Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 8 Nov 2023 11:05:20 +0100 Subject: [PATCH 1/9] [#870] plot_line_on_det_tracing: can plot mutiple wavelegths, can change pixels coords to lower origin, can convert line+Rcs traces on detector to 2D image --- tofu/geom/_core_optics.py | 423 ++++++++++++++++++++++++++++---------- 1 file changed, 309 insertions(+), 114 deletions(-) diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index 957984a71..1b3c3d7e7 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -1946,29 +1946,50 @@ def calc_xixj_from_braggphi( def plot_line_on_det_tracing( self, + # ----------------------- # Options of basic method dcryst=None, - n=None, nphi2=None, - det=None, johann=None, - lpsi=None, ldtheta=None, + n=None, + nphi=None, + det=None, + johann=None, + lpsi=None, + ldtheta=None, + # ----------------------- # Type of crystal - crystal=None, din=None, + crystal=None, + din=None, + # ----------------------- # Wavelength lamb=None, + # ----------------------- # Options of crystal modifications merge_rc_data=None, miscut=None, therm_exp=None, - alpha_limits=None, na=None, - alpha0=None, temp0=None, + alpha_limits=None, + na=None, + alpha0=None, + temp0=None, temp_limits=None, + # ----------------------- # Plot plot_rcs=None, + plot_line_tracing=None, + plot_perfect=None, + plot_simu_image=None, + mode=None, + nxi=None, + nxj=None, strict=None, - plot=None, ax=None, - dleg=None, color=None, - rocking=None, fs=None, dmargin=None, - wintit=None, tit=None, + ax=None, + dleg=None, + color=None, + rocking=None, + fs=None, + dmargin=None, + wintit=None, + tit=None, ): """ Visualize the de-focusing by ray-tracing of chosen lamb Possibility to plot few wavelength' arcs on the same plot. @@ -1998,7 +2019,9 @@ def plot_line_on_det_tracing( By default to 10°C so temp0=35 """ + # ----------------------- # Check / format inputs + lok = [ k0 for k0 in _rockingcurve_def._DCRYST.keys() if 'xxx' not in k0.lower() @@ -2010,20 +2033,54 @@ def plot_line_on_det_tracing( ) din = _rockingcurve_def._DCRYST[crystal] - if merge_rc_data is None: - merge_rc_data = False - if lamb is None and merge_rc_data is False: + dlamb = {} + if lamb is None: lamb = self._dbragg['lambref'] - elif lamb is None and merge_rc_data is True: - # He-like resonance line w at 3.969067 A, intercombination lines - # line x at 3.965858A and line y at 3.969356A, forbidden line z at - # 3.994145A; Li-like dielectronic satellite line k at 3.98981A - lamb = np.r_[ - 3.949067e-10, 3.965858e-10, 3.969356e-10, - 3.994145e-10, 3.989810e-10, - ] + elif lamb == 'ArXVII': + dlamb = { + 'ArXVII_w_Bruhns': 3.949065e-10, + 'ArXV_n3_Adhoc200408': 3.956000e-10, + 'ArXVII_x_Adhoc200408':3.965857e-10, + 'ArXVII_y_Adhoc200408': 3.969356e-10, + 'ArXVI_q_Adhoc200408': 3.981300e-10, + 'ArXVI_r_Adhoc200408': 3.983400e-10, + 'ArXVI_a_Adhoc200408': 3.984570e-10, + 'ArXVI_k_Adhoc200408': 3.989800e-10, + 'ArXVI_j_Adhoc200408': 3.993800e-10, + 'ArXVII_z_Amaro': 3.994130e-10, + 'WXLIV_0_Adhoc20210119': 3.963000e-10, + 'WXLIV_1_Adhoc20210119': 3.975000e-10, + 'WXLIV_2_Adhoc20210119': 3.988670e-10, + } + values = list(dlamb.values()) + lamb = np.asarray(values) + elif lamb == 'ArXVII-woW': + dlamb = { + 'ArXVII_w_Bruhns': 3.949065e-10, + 'ArXV_n3_Adhoc200408': 3.956000e-10, + 'ArXVII_x_Adhoc200408':3.965857e-10, + 'ArXVII_y_Adhoc200408': 3.969356e-10, + 'ArXVI_q_Adhoc200408': 3.981300e-10, + 'ArXVI_r_Adhoc200408': 3.983400e-10, + 'ArXVI_a_Adhoc200408': 3.984570e-10, + 'ArXVI_k_Adhoc200408': 3.989800e-10, + 'ArXVI_j_Adhoc200408': 3.993800e-10, + 'ArXVII_z_Amaro': 3.994130e-10, + } + values = list(dlamb.values()) + lamb = np.asarray(values) + elif lamb == 'ArXVII-wxyz': + dlamb = { + 'ArXVII_w_Bruhns': 3.949065e-10, + 'ArXVII_x_Adhoc200408':3.965857e-10, + 'ArXVII_y_Adhoc200408': 3.969356e-10, + 'ArXVII_z_Amaro': 3.994130e-10, + } + values = list(dlamb.values()) + lamb = np.asarray(values) lamb = np.atleast_1d(lamb).ravel() nlamb = lamb.size + if miscut is None: miscut = False if therm_exp is None: @@ -2033,9 +2090,12 @@ def plot_line_on_det_tracing( 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] + alpha_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: na = 41 nn = (na/2.) @@ -2043,30 +2103,57 @@ def plot_line_on_det_tracing( nn = int(nn - 1) else: nn = int(nn - 0.5) + if alpha0 is None: alpha0 = (3/60)*np.pi/180. if temp0 is None: - temp0 = 10. + temp0 = 15. + if det is None or det.get('outline') is None: msg = ("Please provide det as a dict with 'outline'!") raise Exception(msg) + if plot_rcs is None: plot_rcs = False - if plot is None: - plot = True + if plot_line_tracing is None: + plot_line_tracing = True + if plot_perfect is None: + plot_perfect = True + + if merge_rc_data is None: + merge_rc_data = False + + if plot_simu_image is None and merge_rc_data: + plot_simu_image = True + elif plot_simu_image is None and not merge_rc_data: + plot_simu_image = False + elif plot_simu_image and not merge_rc_data: + plot_simu_image = False + if strict is None: strict = True + if mode is None: + mode = None + if nxi is None: + nxi = 487 + if nxj is None: + nxj = 1467 + # ----------------------- # Check from args inputs the values of amplitude miscut angle alpha and # inter-reticular spacing + 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, + crystal=crystal, + din=din, lamb=self.dbragg['lambref']*1e10, - na=na, nn=nn, + na=na, + nn=nn, therm_exp=therm_exp, temp_limits=temp_limits, plot_therm_exp=False, @@ -2087,16 +2174,23 @@ def find_nearest(array, value): id_temp0 = find_nearest(TD, temp0) self.dmat['d'] = d_atom[id_temp0]*1e-10 + # ----------------------- # Get local basis + nout, e1, e2, miscut = self.get_unit_vectors( miscut=miscut, ) nin = -nout + # ----------------------- # Compute lamb / phi + _, phi = self.get_lambbraggphi_from_ptsxixj_dthetapsi( - xi=det['outline'][0, :], xj=det['outline'][1, :], det=det, - dtheta=0, psi=0, + xi=det['outline'][0, :], + xj=det['outline'][1, :], + det=det, + dtheta=0, + psi=0, miscut=miscut, n=n, grid=True, @@ -2105,11 +2199,12 @@ 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: - nphi2 = 50 - nphi = 2*nphi2 + if nphi is None: + nphi = 100 phi = np.linspace(phimin, phimax, nphi) xi = np.full((nlamb, nphi), np.nan) @@ -2127,7 +2222,9 @@ def find_nearest(array, value): plot=False, ) + # ----------------------- # Get johann-error raytracing (multiple positions on crystal) + xi_er, xj_er = None, None if johann and not rocking: if lpsi is None: @@ -2156,37 +2253,93 @@ def find_nearest(array, value): strict=strict, ) + # ----------------------- # Get rocking curve error if rocking: pass + # ----------------------- # Picking the number of points used to compute a rocking curve & their # glancing angles associated, computing the coordinates (xi_rc, xj_rc) # related to plot the wavelength arc with a transparency parameter # 'alpha' (cf.plt.plot()) corresponding to the diffracted intensity # value at this glancing angle. + if merge_rc_data: - xi_rc = np.full((1), np.nan) + + # First compute_rockingcurve() for output arrays sizing + dout = _rockingcurve.compute_rockingcurve( + crystal=crystal, + din=din, + lamb=lamb[0]*1e10, + miscut=miscut, + therm_exp=therm_exp, + temp_limits=temp_limits, + plot_therm_exp=plot_rcs, + alpha_limits=alpha_limits, + nn=None, + plot_asf=False, + plot_power_ratio=plot_rcs, + plot_asymmetry=False, + plot_cmaps=False, + returnas=dict, + ) + if miscut and therm_exp: + dth = dout['Glancing angles'][0, id_temp0, id_alpha0, :] + ndth = dth.size + elif not miscut and not therm_exp: + dth = dout['Glancing angles'][0, 0, 0, :] + ndth = dth.size + elif miscut and not therm_exp: + dth = dout['Glancing angles'][0, 0, id_alpha0, :] + ndth = dth.size + elif not miscut and therm_exp: + dth = dout['Glancing angles'][0, id_temp0, 0, :] + ndth = dth.size + + # For each wavelength, get results dictionnary of the associated + # diffraction pattern + power_ratio = np.full(( + nlamb, + dout['Power ratio'].shape[0], + dout['Power ratio'].shape[1], + dout['Power ratio'].shape[2], + dout['Power ratio'].shape[3], + ), np.nan) + xi_rc = np.full((nlamb, ndth, nphi), np.nan) xj_rc = xi_rc.copy() - power_ratio = xi_rc.copy() - xi_atprmax = xi_rc.copy() - xj_atprmax = xi_rc.copy() + xi_atprmax = np.full((nlamb, 1), np.nan) + xj_atprmax = xi_atprmax.copy() bragg_atprmax = xi_atprmax.copy() lamb_atprmax = xi_atprmax.copy() - # For each wavelength, get results dictionnary of the associated - # diffraction pattern + pix_horiz = np.linspace( + det['outline'][0, 0], + det['outline'][0, 1], + nxi + ) + pix_verti = np.linspace( + det['outline'][1, 1], + det['outline'][1, 2], + nxj + ) + data = np.full((1, pix_verti.size, pix_horiz.size), 0) + for ll in range(nlamb): dout = _rockingcurve.compute_rockingcurve( - crystal=crystal, din=din, + crystal=crystal, + din=din, lamb=lamb[ll]*1e10, miscut=miscut, therm_exp=therm_exp, temp_limits=temp_limits, plot_therm_exp=plot_rcs, - alpha_limits=alpha_limits, nn=None, - plot_asf=False, plot_power_ratio=plot_rcs, - plot_asymmetry=False, plot_cmaps=False, + alpha_limits=alpha_limits, + nn=None, + plot_asf=False, + plot_power_ratio=plot_rcs, + plot_asymmetry=False, + plot_cmaps=False, returnas=dict, ) TD = np.zeros((na,), dtype=float) @@ -2197,26 +2350,10 @@ def find_nearest(array, value): if miscut: angles = dout['Miscut angles (deg)'] nangles = angles.size - power_ratio = np.resize(power_ratio, ( - nlamb, - dout['Power ratio'].shape[0], - dout['Power ratio'].shape[1], - dout['Power ratio'].shape[2], - dout['Power ratio'].shape[3], - ) - ) power_ratio[ll, ...] = dout['Power ratio'] - def find_nearest(array, value): - array = np.asarray(array) - idx = (np.abs(array - value)).argmin() - return idx - id_alpha0 = find_nearest(angles, alpha0) - # Pull the glancing angles 'dth' & the number of points 'ndth' - # depending on the case related to unp & therm_exp, plus - # find the glancing angle related the max power ratio value if miscut and therm_exp: dth = dout['Glancing angles'][0, id_temp0, id_alpha0, :] ndth = dth.size @@ -2254,14 +2391,6 @@ def find_nearest(array, value): ) dth_atprmax = dth[ind_pr_max] - # Resize results arrays - xi_rc = np.resize(xi_rc, (nlamb, ndth, nphi)) - xj_rc = xi_rc.copy() - xi_atprmax = np.resize(xi_atprmax, (nlamb, 1)) - xj_atprmax = xi_atprmax.copy() - bragg_atprmax = xi_atprmax.copy() - lamb_atprmax = xi_atprmax.copy() - # Compute wavelength arcs for each glancing angle to obtain # the shadow of the diffraction pattern on the detector for mm in range(ndth): @@ -2278,8 +2407,24 @@ def find_nearest(array, value): strict=strict, plot=False, ) - xi_atprmax[ll] = xi_rc[ll, ind_pr_max, nphi2] - xj_atprmax[ll] = xj_rc[ll, ind_pr_max, nphi2] + if plot_simu_image: + for nn in range(phi.size): + if np.isfinite(xi_rc[ll, mm, nn]): + idxi = np.nanargmin( + np.abs( + xi_rc[ll, mm, nn] - pix_horiz + ) + ) + idxj = np.nanargmin( + np.abs( + xj_rc[ll, mm, nn] - pix_verti + ) + ) + data[0, idxj, idxi] += 1 + + + xi_atprmax[ll] = xi_rc[ll, ind_pr_max, int(nphi/2)] + xj_atprmax[ll] = xj_rc[ll, ind_pr_max, int(nphi/2)] self.update_miscut(alpha=0., beta=0.) if therm_exp: self.dmat['d'] = d_atom[nn]*1e-10 @@ -2295,6 +2440,20 @@ def find_nearest(array, value): grid=True, return_lamb=True, ) + else: + power_ratio=None + dth=None + ndth=None + nn=None + xi_rc=None + xj_rc=None + xi_atprmax=None + bragg_atprmax=None + lamb_atprmax=None + TD=None + angles=None + data=None + # Reset parameters as at beginning if miscut: @@ -2306,52 +2465,88 @@ def find_nearest(array, value): else: self.dmat['d'] = d_atom[0]*1e-10 - # Plot - if plot: - if merge_rc_data: - return _plot_optics.CrystalBragg_plot_line_tracing_on_det( - cryst=self, dcryst=dcryst, - lamb=lamb, - xi=xi, xj=xj, xi_er=xi_er, xj_er=xj_er, - power_ratio=power_ratio, dth=dth, ndth=ndth, nn=nn, - xi_rc=xi_rc, xj_rc=xj_rc, - xi_atprmax=xi_atprmax, - bragg_atprmax=bragg_atprmax, - lamb_atprmax=lamb_atprmax, - det=det, - johann=johann, rocking=rocking, - miscut=miscut, - therm_exp=therm_exp, - merge_rc_data=merge_rc_data, - alpha0=alpha0, temp0=temp0, - TD=TD, angles=angles, - id_temp0=id_temp0, - ax=ax, dleg=dleg, color=color, - fs=fs, dmargin=dmargin, wintit=wintit, tit=tit, - ) - else: - return _plot_optics.CrystalBragg_plot_line_tracing_on_det( - cryst=self, dcryst=dcryst, - lamb=lamb, xi=xi, xj=xj, xi_er=xi_er, xj_er=xj_er, - alpha0=alpha0, temp0=temp0, - id_temp0=id_temp0, - johann=johann, rocking=rocking, - miscut=miscut, - therm_exp=therm_exp, - merge_rc_data=merge_rc_data, - det=det, - ax=ax, dleg=dleg, color=color, - fs=fs, dmargin=dmargin, wintit=wintit, tit=tit, - ) + if mode == 'raw det': + xi = xi - det['outline'][0, 0] + xj = xj - det['outline'][1, 0] + xi_rc = xi_rc - det['outline'][0, 0] + xj_rc = xj_rc - det['outline'][1, 0] + det['outline'][0] = det['outline'][0] - det['outline'][0, 0] + det['outline'][1] = det['outline'][1] - det['outline'][1, 0] + + if plot_simu_image: + fig0 = plt.figure() + ax0 = fig0.add_subplot() + ax0.set_xlabel(r'Pixel coordinate', fontsize=15) + ax0.set_ylabel(r'Pixel coordinate', fontsize=15) + ax0.set_title(r'Simulated 2D spectra', fontsize=15) + ax0.imshow( + data[0, :, :], + origin='lower', + interpolation='nearest', + aspect='auto', + ) + + dout = { + 'lamb': lamb, + 'xi': xi, + 'xj': xj, + 'xi_rc': xi_rc, + 'xj_rc': xj_rc, + 'xi_atprmax': xi_atprmax, + 'lamb_atprmax': lamb_atprmax, + 'bragg_atprmax': bragg_atprmax, + 'data': data, + } + + if plot_line_tracing: + ax = _plot_optics.CrystalBragg_plot_line_tracing_on_det( + # ------------------------------ + # basic + cryst=self, + dcryst=dcryst, + lamb=lamb, + dlamb=dlamb, + xi=xi, + xj=xj, + xi_er=xi_er, + xj_er=xj_er, + # ----------------------------- + # w/ rocking curves data + merge_rc_data=merge_rc_data, + power_ratio=power_ratio, + dth=dth, + ndth=ndth, + nn=nn, + xi_rc=xi_rc, + xj_rc=xj_rc, + xi_atprmax=xi_atprmax, + bragg_atprmax=bragg_atprmax, + lamb_atprmax=lamb_atprmax, + TD=TD, + angles=angles, + # ----------------------------- + # w/ miscut and/or temp changes + alpha0=alpha0, + temp0=temp0, + id_temp0=id_temp0, + johann=johann, + rocking=rocking, + miscut=miscut, + therm_exp=therm_exp, + det=det, + # ---------------------------- + # plot parameters + ax=ax, + dleg=dleg, + color=color, + fs=fs, + dmargin=dmargin, + wintit=wintit, + tit=tit, + plot_perfect=plot_perfect, + ) + return ax, dout else: - dout = {'lamb': lamb, - 'xi': xi, - 'xj': xj, - 'xi_rc': xi_rc, - 'xj_rc': xj_rc, - 'xi_atprmax': xi_atprmax, - 'lamb_atprmax': lamb_atprmax, - 'bragg_atprmax': bragg_atprmax} return dout def comp_angular_shift_on_det_tracing( @@ -2854,7 +3049,7 @@ def _get_local_coordinates_of_det( det_approx = self.get_detector_ideal( bragg=bragg, lamb=lamb, - tangent_to_rowland=False, + tangent_to_rowland=True,#False, miscut=miscut, ) From 837c4d6b67063f7aa8a92613aa7c6c9e43806130 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 8 Nov 2023 11:07:39 +0100 Subject: [PATCH 2/9] [#870] plot_line_tracing_on_det: legend out of plot, colors for each wavelength, add if already exists name of wavelength (dlamb dict) --- tofu/geom/_plot_optics.py | 173 +++++++++++++++++++++++++++----------- 1 file changed, 125 insertions(+), 48 deletions(-) diff --git a/tofu/geom/_plot_optics.py b/tofu/geom/_plot_optics.py index f501d4bc8..82306bf74 100644 --- a/tofu/geom/_plot_optics.py +++ b/tofu/geom/_plot_optics.py @@ -946,42 +946,66 @@ def CrystalBragg_plot_braggangle_from_xixj(xi=None, xj=None, def CrystalBragg_plot_line_tracing_on_det( - cryst=None, dcryst=None, + cryst=None, + dcryst=None, lamb=None, - xi=None, xj=None, xi_er=None, xj_er=None, - power_ratio=None, dth=None, ndth=None, nn=None, - xi_rc=None, xj_rc=None, + dlamb=None, + xi=None, + xj=None, + xi_er=None, + xj_er=None, + power_ratio=None, + dth=None, + ndth=None, + nn=None, + xi_rc=None, + xj_rc=None, xi_atprmax=None, bragg_atprmax=None, lamb_atprmax=None, det=None, - johann=None, rocking=None, + johann=None, + rocking=None, miscut=None, therm_exp=None, merge_rc_data=None, - alpha0=None, temp0=None, TD=None, angles=None, + alpha0=None, + temp0=None, + TD=None, + angles=None, id_temp0=None, - ax=None, dleg=None, color=None, - fs=None, dmargin=None, wintit=None, tit=None, + ax=None, + dleg=None, + color=None, + fs=None, + dmargin=None, + wintit=None, + tit=None, + plot_perfect=None, ): - # Check inputs # ------------ + # Check inputs if dleg is None: dleg = { 'loc': 'upper left', - 'fontsize': 13, + 'fontsize': 10, + 'bbox_to_anchor': (1., 1.), } if color is None: color = 'k' - if fs is None: fs = (8, 8) if dmargin is None: - dmargin = {'left': 0.15, 'right': 0.95, - 'bottom': 0.08, 'top': 0.92, - 'wspace': None, 'hspace': 0.4} + dmargin = { + 'top': 0.950, + 'bottom': 0.070, + 'left': 0.070, + 'right': 0.750, + 'hspace': 0.200, + 'wspace': 0.200, + } if wintit is None: wintit = _WINTIT @@ -993,11 +1017,19 @@ def CrystalBragg_plot_line_tracing_on_det( tit += " - rocking curve" plot_err = johann is True or rocking is True - markers = ['o', '^', 'D', 's', 'X'] - colors = ['r', 'g', 'c', 'b', 'k'] + markers = [ + 'o', 'v', '8', 's', 'P', + '*', 'H', '+', 'X', 'D', + '^', 'p', 'h', 'x', 'd', + ] + colors = [ + 'b', 'g', 'r', 'c', 'm', + 'y', 'k', 'orange', 'purple', 'brown', + 'pink', 'gray', 'olive', 'coral', 'royalblue', + ] - # Plot # ------------ + # Plot if ax is None: fig = plt.figure(figsize=fs) @@ -1009,6 +1041,7 @@ def CrystalBragg_plot_line_tracing_on_det( fig.suptitle(tit, size=14, weight='bold') ax.set_xlabel(r'Pixel coordinate $x_{i}$ [m]', fontsize=15) ax.set_ylabel(r'Pixel coordinate $x_{j}$ [m]', fontsize=15) + ax.set_xlim( det['outline'][0, :].min() - 0.01, det['outline'][0, :].max() + 0.01, @@ -1017,11 +1050,14 @@ def CrystalBragg_plot_line_tracing_on_det( det['outline'][1, :].min() - 0.01, det['outline'][1, :].max() + 0.01, ) + if det.get('outline') is not None: ax.plot( - det['outline'][0, :], det['outline'][1, :], + det['outline'][0, :], + det['outline'][1, :], ls='-', lw=1., c='k', ) + aa = np.r_[cryst.dmat['alpha']] if therm_exp and merge_rc_data: bb = TD[id_temp0] @@ -1029,43 +1065,84 @@ def CrystalBragg_plot_line_tracing_on_det( bb = temp0 else: bb = 0. - for ll in range(lamb.size): - lab = ( - r'$\lambda$ = {} A'.format(np.round(lamb[ll]*1e10, 6)) + '\n' - + r'$\Delta$T = {} °C, $\alpha$ = {} deg'.format( - bb, aa[0]*(180./np.pi) - ) - ) - l0, = ax.plot( - xi[ll, :], xj[ll, :], - ls='--', lw=1., - marker=markers[ll], ms=4., - c=color, - label=lab, - ) - if plot_err: - ax.plot( - xi_er[ll, ...], xj_er[ll, ...], - ls='None', lw=1., c=l0.get_color(), - ms=4, marker='.', + + if plot_perfect: + for kk in range(lamb.size): + lkey = list(dlamb.keys()) + lval = list(dlamb.values()) + try: + pos = lval.index(lamb[kk]) + key = lkey[pos][:8] + lab = ( + r'$\lambda$ = {} A ({}), $\Delta$T = {}°C, $\alpha$ = {}deg'.format( + np.round(lamb[kk]*1e10, 4), + key, + bb, + aa[0]*(180./np.pi), + ) + ) + except Exception as err: + lab = ( + r'$\lambda$ = {} A, $\Delta$T = {}°C, $\alpha$ = {}deg'.format( + np.round(lamb[kk]*1e10, 4), + bb, + aa[0]*(180./np.pi), + ) + ) + l0, = ax.plot( + xi[kk, :], + xj[kk, :], + ls='--', + lw=1., + marker=markers[kk], + ms=8., + c=colors[kk], + label=lab, ) + if plot_err: + ax.plot( + xi_er[kk, ...], + xj_er[kk, ...], + ls='None', + lw=1., + c=colors[kk], + ms=8, + marker='.', + ) + if merge_rc_data: - for ll in range(lamb.size): + for kk in range(lamb.size): for mm in range(ndth): if mm == int(ndth/2.): - label = r'At $x_j$=0.: $x_i$={}, $\lambda$={}A'.format( - np.round(xi_atprmax[ll], 6), - np.round(lamb_atprmax[ll], 16), - # np.round(bragg_atprmax[ll]*(180./np.pi), 4), - ) + lkey = list(dlamb.keys()) + lval = list(dlamb.values()) + try: + pos = lval.index(lamb[kk]) + key = lkey[pos][:8] + label = ( + r'At $x_j$=0.: $x_i$={}, $\lambda$={} A ({})'.format( + np.round(xi_atprmax[kk][0], 4), + np.round(lamb_atprmax[kk][0], 14), + key, + ) + ) + except Exception as err: + label = ( + r'At $x_j$=0.: $x_i$={}, $\lambda$={} A'.format( + np.round(xi_atprmax[kk][0], 4), + np.round(lamb_atprmax[kk][0], 14), + ) + ) else: label = None - pr1 = power_ratio[ll, 0, 0, 0, mm] - pr2 = power_ratio[ll, 1, 0, 0, mm] + pr1 = power_ratio[kk, 0, 0, 0, mm] + pr2 = power_ratio[kk, 1, 0, 0, mm] ax.plot( - xi_rc[ll, mm, :], xj_rc[ll, mm, :], - ls='-', lw=1., - c=l0.get_color(), + xi_rc[kk, mm, :], + xj_rc[kk, mm, :], + ls='-', + lw=1., + c=colors[kk], alpha=pr1 + pr2, label=label, ) From 01f6b9ae8ecab17f80dbe6fcfc557b1b104e6726 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Thu, 9 Nov 2023 12:27:27 +0100 Subject: [PATCH 3/9] [#870] plot_line_on_det_tracing: add crystal splitting option --- tofu/geom/_core_optics.py | 119 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 113 insertions(+), 6 deletions(-) diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index 1b3c3d7e7..0324f8ae8 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -1959,9 +1959,14 @@ def plot_line_on_det_tracing( # Type of crystal crystal=None, din=None, + split=None, + direction=None, + nb=None, + which=None, # ----------------------- # Wavelength lamb=None, + subarea=None, # ----------------------- # Options of crystal modifications merge_rc_data=None, @@ -1990,6 +1995,7 @@ def plot_line_on_det_tracing( dmargin=None, wintit=None, tit=None, + save=None, ): """ Visualize the de-focusing by ray-tracing of chosen lamb Possibility to plot few wavelength' arcs on the same plot. @@ -2022,6 +2028,102 @@ def plot_line_on_det_tracing( # ----------------------- # Check / format inputs + # ----------------------- + # crystal specific + + if split is None: + split = False + if direction is None and not split: + direction = None + elif direction is None and split: + direction = None + msg = ( + "You choose split = True and direction = None\n" + + " => direction = e1 by default in split() !" + ) + warnings.warn(msg) + direction = 'e1' + + if nb is None and not split: + nb = None + elif nb is None and split: + nb = None + msg = ( + "You choose split = True and nb = None\n" + + " => nb = 2 by default in split() !" + ) + warnings.warn(msg) + nb = 2 + elif nb > 2 and split: + msg = ( + "You cant for now choose more than 2 pieces of crystal, sorry." + ) + raise Exception(msg) + + if split: + self1, self2 = self.split( + direction=direction, + nb=nb, + ) + + if which is None and not split: + which = None + elif which is None and split: + msg = ( + "You must choose if you work w/ both halves or only one !\n" + + "Please choose either :\n" + + "\t - both not alowed for now\n" + + "\t - left or right if direction = e1 & nb = 2\n" + + "\t - up or down if direction = e2 & nb = 2\n" + ) + raise Exception(msg) + + ldir = ['e1', 'e2'] + lwhich = ['both', 'left', 'right', 'up', 'down'] + if split: + if direction == ldir[0] and which == lwhich[0]: + msg = "Choose only 1 halve please !" + raise Exception(msg) + # lself = [self1, self2] + + elif direction == ldir[0] and which == lwhich[1]: + self = self1 + + elif direction == ldir[0] and which == lwhich[2]: + self = self2 + + elif direction == ldir[0] and which == lwhich[3]: + msg = "Choose between left or right along e1 direction !" + raise Exception(msg) + + elif direction == ldir[0] and which == lwhich[4]: + msg = "Choose between left or right along e1 direction !" + raise Exception(msg) + + if direction == ldir[1] and which == lwhich[0]: + msg = "Choose only 1 halve please !" + raise Exception(msg) + # lself = [self1, self2] + + elif direction == ldir[1] and which == lwhich[3]: + self = self2 + + elif direction == ldir[1] and which == lwhich[4]: + self = self1 + + elif direction == ldir[1] and which == lwhich[1]: + msg = "Choose between up or down along e2 direction !" + raise Exception(msg) + + elif direction == ldir[1] and which == lwhich[2]: + msg = "Choose between up or down along e2 direction !" + raise Exception(msg) + else: + self = self + + # -------------------- + # other generic checks + lok = [ k0 for k0 in _rockingcurve_def._DCRYST.keys() if 'xxx' not in k0.lower() @@ -2034,9 +2136,9 @@ def plot_line_on_det_tracing( din = _rockingcurve_def._DCRYST[crystal] dlamb = {} - if lamb is None: + if lamb is None and subarea is None: lamb = self._dbragg['lambref'] - elif lamb == 'ArXVII': + elif lamb is None and subarea == 'ArXVII': dlamb = { 'ArXVII_w_Bruhns': 3.949065e-10, 'ArXV_n3_Adhoc200408': 3.956000e-10, @@ -2054,7 +2156,7 @@ def plot_line_on_det_tracing( } values = list(dlamb.values()) lamb = np.asarray(values) - elif lamb == 'ArXVII-woW': + elif lamb is None and subarea == 'ArXVII-woW': dlamb = { 'ArXVII_w_Bruhns': 3.949065e-10, 'ArXV_n3_Adhoc200408': 3.956000e-10, @@ -2069,7 +2171,7 @@ def plot_line_on_det_tracing( } values = list(dlamb.values()) lamb = np.asarray(values) - elif lamb == 'ArXVII-wxyz': + elif lamb is None and subarea == 'ArXVII-wxyz': dlamb = { 'ArXVII_w_Bruhns': 3.949065e-10, 'ArXVII_x_Adhoc200408':3.965857e-10, @@ -2139,6 +2241,9 @@ def plot_line_on_det_tracing( if nxj is None: nxj = 1467 + if save is None: + save = False + # ----------------------- # Check from args inputs the values of amplitude miscut angle alpha and # inter-reticular spacing @@ -2147,7 +2252,6 @@ def plot_line_on_det_tracing( 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, @@ -2266,7 +2370,6 @@ def find_nearest(array, value): # value at this glancing angle. if merge_rc_data: - # First compute_rockingcurve() for output arrays sizing dout = _rockingcurve.compute_rockingcurve( crystal=crystal, @@ -2498,6 +2601,10 @@ def find_nearest(array, value): 'data': data, } + # save + if save: + None + if plot_line_tracing: ax = _plot_optics.CrystalBragg_plot_line_tracing_on_det( # ------------------------------ From badaf9fe6c42c2c93ea38cdc4c03a5eb3f91bac1 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Tue, 21 Nov 2023 15:53:20 +0100 Subject: [PATCH 4/9] [#870] _core_optics: add Doppler shift option, subarea names changed --- tofu/geom/_core_optics.py | 149 ++++++++++++++++++++++++++++++-------- 1 file changed, 120 insertions(+), 29 deletions(-) diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index 0324f8ae8..6480e3e23 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -1967,6 +1967,8 @@ def plot_line_on_det_tracing( # Wavelength lamb=None, subarea=None, + shift=None, + v=None, # ----------------------- # Options of crystal modifications merge_rc_data=None, @@ -2136,13 +2138,16 @@ def plot_line_on_det_tracing( din = _rockingcurve_def._DCRYST[crystal] dlamb = {} + lsubarea = [ + 'ArXVII-simul', 'ArXVII-woW', 'ArXVII-wxyz' + ] if lamb is None and subarea is None: lamb = self._dbragg['lambref'] - elif lamb is None and subarea == 'ArXVII': + elif lamb is None and subarea == lsubarea[0]: dlamb = { 'ArXVII_w_Bruhns': 3.949065e-10, 'ArXV_n3_Adhoc200408': 3.956000e-10, - 'ArXVII_x_Adhoc200408':3.965857e-10, + 'ArXVII_x_Adhoc200408': 3.965857e-10, 'ArXVII_y_Adhoc200408': 3.969356e-10, 'ArXVI_q_Adhoc200408': 3.981300e-10, 'ArXVI_r_Adhoc200408': 3.983400e-10, @@ -2156,11 +2161,11 @@ def plot_line_on_det_tracing( } values = list(dlamb.values()) lamb = np.asarray(values) - elif lamb is None and subarea == 'ArXVII-woW': + elif lamb is None and subarea == lsubarea[1]: dlamb = { 'ArXVII_w_Bruhns': 3.949065e-10, 'ArXV_n3_Adhoc200408': 3.956000e-10, - 'ArXVII_x_Adhoc200408':3.965857e-10, + 'ArXVII_x_Adhoc200408': 3.965857e-10, 'ArXVII_y_Adhoc200408': 3.969356e-10, 'ArXVI_q_Adhoc200408': 3.981300e-10, 'ArXVI_r_Adhoc200408': 3.983400e-10, @@ -2171,18 +2176,57 @@ def plot_line_on_det_tracing( } values = list(dlamb.values()) lamb = np.asarray(values) - elif lamb is None and subarea == 'ArXVII-wxyz': + elif lamb is None and subarea == lsubarea[2]: dlamb = { 'ArXVII_w_Bruhns': 3.949065e-10, - 'ArXVII_x_Adhoc200408':3.965857e-10, + 'ArXVII_x_Adhoc200408': 3.965857e-10, 'ArXVII_y_Adhoc200408': 3.969356e-10, 'ArXVII_z_Amaro': 3.994130e-10, } values = list(dlamb.values()) lamb = np.asarray(values) + elif lamb is None and subarea not in lsubarea: + msg = ( + "Please choose one of the following subareas :\n" + + "\t - {} for all Ar and W lines \n".format(lsubarea[0]) + + "\t - {} for all Ar lines \n".format(lsubarea[1]) + + "\t - {} for only Ar 16+ lines \n".format(lsubarea[2]) + ) + raise Exception(msg) + lamb = np.atleast_1d(lamb).ravel() nlamb = lamb.size + # --------------------- + # Doppler shift option + + def doppler_shift(lamb, v, c): + newlamb = np.full((lamb.size), np.nan) + for ii in range(lamb.size): + dl = lamb[ii] * (v/c) + newlamb[ii] = lamb[ii] - dl + return newlamb + + if shift is None: + shift = False + elif shift: + msg = ( + "You have to give a velocity in meters per second !" + "You want a velocity of {} m/s !" + ) + warnings.warn(msg) + c = 2.99792458e8 # m/s + newlamb = doppler_shift( + lamb=lamb, + v=v, + c=c, + ) + assert newlamb.size == lamb.size + lamb = newlamb + + # ------------------- + # other checks + if miscut is None: miscut = False if therm_exp is None: @@ -2235,7 +2279,7 @@ def plot_line_on_det_tracing( if strict is None: strict = True if mode is None: - mode = None + mode = None if nxi is None: nxi = 487 if nxj is None: @@ -2387,6 +2431,27 @@ def find_nearest(array, value): plot_cmaps=False, returnas=dict, ) + + TD = np.zeros((na,), dtype=float) + if therm_exp: + TD = dout['Temperature changes (°C)'] + nT = TD.size + angles = np.zeros((na,), dtype=float) + if miscut: + angles = dout['Miscut angles (deg)'] + nangles = angles.size + power_ratio = np.full(( + nlamb, + dout['Power ratio'].shape[0], + dout['Power ratio'].shape[1], + dout['Power ratio'].shape[2], + dout['Power ratio'].shape[3], + ), np.nan) + power_ratio[ll, ...] = dout['Power ratio'] + + id_alpha0 = find_nearest(angles, alpha0) + id_temp0 = find_nearest(TD, temp0) + if miscut and therm_exp: dth = dout['Glancing angles'][0, id_temp0, id_alpha0, :] ndth = dth.size @@ -2426,7 +2491,8 @@ def find_nearest(array, value): det['outline'][1, 2], nxj ) - data = np.full((1, pix_verti.size, pix_horiz.size), 0) + # data = np.full((1, pix_verti.size, pix_horiz.size), 0) + data = np.empty([1, pix_verti.size, pix_horiz.size]) for ll in range(nlamb): dout = _rockingcurve.compute_rockingcurve( @@ -2456,6 +2522,7 @@ def find_nearest(array, value): power_ratio[ll, ...] = dout['Power ratio'] id_alpha0 = find_nearest(angles, alpha0) + id_temp0 = find_nearest(TD, temp0) if miscut and therm_exp: dth = dout['Glancing angles'][0, id_temp0, id_alpha0, :] @@ -2523,16 +2590,18 @@ def find_nearest(array, value): xj_rc[ll, mm, nn] - pix_verti ) ) - data[0, idxj, idxi] += 1 - + pr1 = power_ratio[ll, 0, 0, 0, mm] + pr2 = power_ratio[ll, 1, 0, 0, mm] + val = pr1 + pr2 + data[0, idxj, idxi] += val xi_atprmax[ll] = xi_rc[ll, ind_pr_max, int(nphi/2)] xj_atprmax[ll] = xj_rc[ll, ind_pr_max, int(nphi/2)] - self.update_miscut(alpha=0., beta=0.) - if therm_exp: - self.dmat['d'] = d_atom[nn]*1e-10 - else: - self.dmat['d'] = d_atom[0]*1e-10 + # self.update_miscut(alpha=0., beta=0.) + # if therm_exp: + # self.dmat['d'] = d_atom[nn]*1e-10 + # else: + # self.dmat['d'] = d_atom[0]*1e-10 ( bragg_atprmax[ll], _, lamb_atprmax[ll], ) = self.get_lambbraggphi_from_ptsxixj_dthetapsi( @@ -2544,18 +2613,18 @@ def find_nearest(array, value): return_lamb=True, ) else: - power_ratio=None - dth=None - ndth=None - nn=None - xi_rc=None - xj_rc=None - xi_atprmax=None - bragg_atprmax=None - lamb_atprmax=None - TD=None - angles=None - data=None + power_ratio = None + dth = None + ndth = None + nn = None + xi_rc = None + xj_rc = None + xi_atprmax = None + bragg_atprmax = None + lamb_atprmax = None + TD = None + angles = None + data = None # Reset parameters as at beginning @@ -2603,7 +2672,29 @@ def find_nearest(array, value): # save if save: - None + data2 = { + 'data': data, + } + path = '/Home/AD265925/xics/tofu_west/SpectroX2D/' + if miscut: + aa = str(np.round(alpha0 * (180*60/np.pi), 3)) + else: + aa = 'False' + part2 = '_miscut' + aa + if split: + part1 = '_split' + which + else: + part1 = '_splitFalse' + name = path + subarea + part1 + part2 + '_00001.npz' + np.savez( + name, + **data2, + ) + msg = ( + "Data saved at path : {}".format(path + name) + ) + print(msg) + # None if plot_line_tracing: ax = _plot_optics.CrystalBragg_plot_line_tracing_on_det( @@ -3156,7 +3247,7 @@ def _get_local_coordinates_of_det( det_approx = self.get_detector_ideal( bragg=bragg, lamb=lamb, - tangent_to_rowland=True,#False, + tangent_to_rowland=True, #False miscut=miscut, ) From 913ee7d3d9f848475974c9e96829085fdcebf36b Mon Sep 17 00:00:00 2001 From: AD265925 Date: Tue, 21 Nov 2023 15:55:11 +0100 Subject: [PATCH 5/9] [#870] fit12d_dinput: changed indknots determination for simulated spectra, add mode in arguments --- tofu/spectro/_fit12d_dinput.py | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/tofu/spectro/_fit12d_dinput.py b/tofu/spectro/_fit12d_dinput.py index a54bfe697..8340de60f 100644 --- a/tofu/spectro/_fit12d_dinput.py +++ b/tofu/spectro/_fit12d_dinput.py @@ -1550,6 +1550,7 @@ def _dvalid_checkfocus( def fit12d_dvalid( + mode=None, data=None, lamb=None, phi=None, indok_bool=None, binning=None, valid_nsigma=None, valid_fraction=None, @@ -1629,12 +1630,26 @@ def fit12d_dvalid( for ii in range(knots.size - 1): iphi = (phi >= knots[ii]) & (phi < knots[ii + 1]) fract[:, ii, :] = ( - np.sum(np.sum(indall & iphi[None, ..., None], - axis=1), axis=1) - / np.sum(np.sum(iphi[..., None] & lambok, - axis=0), axis=0) + np.sum( + np.sum( + indall & iphi[None, ..., None], + axis=1 + ), + axis=1 + ) + / np.sum( + np.sum( + iphi[..., None] & lambok, + axis=0 + ), + axis=0 + ) ) - indknots = np.all(fract > valid_fraction, axis=2) + + if mode == 'raw': + indknots = np.all(fract > valid_fraction, axis=2) + elif mode == 'simul': + indknots = np.any(fract > valid_fraction, axis=2) # Deduce ldphi ldphi = [[] for ii in range(nspect)] @@ -1988,6 +2003,7 @@ def fit1d_dinput( def fit2d_dinput( + mode=None, dlines=None, dconstraints=None, dconstants=None, dprepare=None, deg=None, nbsplines=None, knots=None, data=None, lamb=None, phi=None, mask=None, @@ -2107,6 +2123,7 @@ def fit2d_dinput( # S/N threshold indices # ------------------------ dinput['valid'] = fit12d_dvalid( + mode=mode, data=dprepare['data'], lamb=dprepare['lamb'], phi=dprepare['phi'], From bb49e40c3a657eb4886719d82cd03c8e9156cc55 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 17 Jan 2024 09:36:02 +0100 Subject: [PATCH 6/9] [#870] core_optics: plot_line_on_det_tracing deal with crystal splitted, plot RCs vs wavelengths option ok, --- tofu/geom/_core_optics.py | 645 +++++++++++++++++++++----------------- 1 file changed, 351 insertions(+), 294 deletions(-) diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index 6480e3e23..1f7a89343 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -1981,6 +1981,7 @@ def plot_line_on_det_tracing( temp_limits=None, # ----------------------- # Plot + plot_as_wavel=None, plot_rcs=None, plot_line_tracing=None, plot_perfect=None, @@ -2082,17 +2083,18 @@ def plot_line_on_det_tracing( ldir = ['e1', 'e2'] lwhich = ['both', 'left', 'right', 'up', 'down'] + lself = [] if split: if direction == ldir[0] and which == lwhich[0]: - msg = "Choose only 1 halve please !" - raise Exception(msg) - # lself = [self1, self2] + # msg = "Choose only 1 halve please !" + # raise Exception(msg) + lself = [self1, self2] elif direction == ldir[0] and which == lwhich[1]: - self = self1 + lself.append(self1) elif direction == ldir[0] and which == lwhich[2]: - self = self2 + lself.append(self2) elif direction == ldir[0] and which == lwhich[3]: msg = "Choose between left or right along e1 direction !" @@ -2121,7 +2123,7 @@ def plot_line_on_det_tracing( msg = "Choose between up or down along e2 direction !" raise Exception(msg) else: - self = self + lself.append(self) # -------------------- # other generic checks @@ -2266,6 +2268,9 @@ def doppler_shift(lamb, v, c): if plot_perfect is None: plot_perfect = True + if plot_as_wavel is None: + plot_as_wavel = False + if merge_rc_data is None: merge_rc_data = False @@ -2288,90 +2293,120 @@ def doppler_shift(lamb, v, c): if save is None: save = False + if fs is None: + fs = 15 + # ----------------------- # Check from args inputs the values of amplitude miscut angle alpha and # inter-reticular spacing - self.update_miscut(alpha=0., beta=0.) - if miscut: - self.update_miscut(alpha=alpha0, beta=0.) + for ii in range(len(lself)): + lself[ii].update_miscut(alpha=0., beta=0.) + if miscut: + lself[ii].update_miscut(alpha=alpha0, beta=0.) - dout = _rockingcurve.CrystBragg_comp_lattice_spacing( - crystal=crystal, - din=din, - lamb=self.dbragg['lambref']*1e10, - na=na, - nn=nn, - therm_exp=therm_exp, - temp_limits=temp_limits, - plot_therm_exp=False, - ) - 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'] - theta = dout['theta_Bragg (rad)'] - theta_deg = dout['theta_Bragg (deg)'] - - def find_nearest(array, value): - array = np.asarray(array) - idx = (np.abs(array - value)).argmin() - return idx - - id_temp0 = find_nearest(TD, temp0) - self.dmat['d'] = d_atom[id_temp0]*1e-10 + dout = _rockingcurve.CrystBragg_comp_lattice_spacing( + crystal=crystal, + din=din, + lamb=lself[ii].dbragg['lambref']*1e10, + na=na, + nn=nn, + therm_exp=therm_exp, + temp_limits=temp_limits, + plot_therm_exp=False, + ) + 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'] + theta = dout['theta_Bragg (rad)'] + theta_deg = dout['theta_Bragg (deg)'] + + def find_nearest(array, value): + array = np.asarray(array) + idx = (np.abs(array - value)).argmin() + return idx + + id_temp0 = find_nearest(TD, temp0) + lself[ii].dmat['d'] = d_atom[id_temp0]*1e-10 # ----------------------- # Get local basis - nout, e1, e2, miscut = self.get_unit_vectors( - miscut=miscut, - ) - nin = -nout + lnout, le1, le2, lnin, lmiscut = [], [], [], [], [] + for ii in range(len(lself)): + ( + nout, e1, e2, miscut + ) = lself[ii].get_unit_vectors( + miscut=miscut, + ) + nin = -nout + lnout.append(nout) + le1.append(e1) + le2.append(e2) + lmiscut.append(miscut) + lnin.append(nin) # ----------------------- # Compute lamb / phi - _, phi = self.get_lambbraggphi_from_ptsxixj_dthetapsi( - xi=det['outline'][0, :], - xj=det['outline'][1, :], - det=det, - dtheta=0, - psi=0, - miscut=miscut, - n=n, - grid=True, - return_lamb=False, - ) - phimin, phimax = np.nanmin(phi), np.nanmax(phi) - phimin, phimax = phimin-(phimax-phimin)/10, phimax+(phimax-phimin)/10 + lphi = [] + lphimin, lphimax = [], [] + for ii in range(len(lself)): + _, phi = lself[ii].get_lambbraggphi_from_ptsxixj_dthetapsi( + xi=det['outline'][0, :], + xj=det['outline'][1, :], + det=det, + dtheta=0, + psi=0, + miscut=miscut, + n=n, + grid=True, + return_lamb=False, + ) + phimin, phimax = np.nanmin(phi), np.nanmax(phi) + phimin = phimin - (phimax - phimin)/10 + phimax = phimax + (phimax - phimin)/10 + lphi.append(phi) + lphimin.append(phimin) + lphimax.append(phimax) # ----------------------- # Get reference ray-tracing - bragg = self._checkformat_bragglamb(lamb=lamb, n=n) - if nphi is None: - nphi = 100 - phi = np.linspace(phimin, phimax, nphi) - - xi = np.full((nlamb, nphi), np.nan) - xj = np.full((nlamb, nphi), np.nan) - for ll in range(nlamb): - xi[ll, :], xj[ll, :] = self.calc_xixj_from_braggphi( - bragg=np.full(phi.shape, bragg[ll]), - phi=phi, - dtheta=0., - psi=0., - n=n, - det=det, - miscut=miscut, - strict=strict, - plot=False, + lxi, lxj = [], [] + lphi2 = [] + for ii in range(len(lself)): + bragg = lself[ii]._checkformat_bragglamb(lamb=lamb, n=n) + if nphi is None: + nphi = 1467 # 100 + phi = np.linspace( + lphimin[ii], lphimax[ii], nphi ) + xi = np.full((nlamb, nphi), np.nan) + xj = np.full((nlamb, nphi), np.nan) + for ll in range(nlamb): + ( + xi[ll, :], xj[ll, :] + ) = lself[ii].calc_xixj_from_braggphi( + bragg=np.full(phi.shape, bragg[ll]), + phi=phi, + dtheta=0., + psi=0., + n=n, + det=det, + miscut=miscut, + strict=strict, + plot=False, + ) + lphi2.append(phi) + lxi.append(xi) + lxj.append(xj) # ----------------------- # Get johann-error raytracing (multiple positions on crystal) + # TBD w/ both crystals xi_er, xj_er = None, None if johann and not rocking: @@ -2414,91 +2449,17 @@ def find_nearest(array, value): # value at this glancing angle. if merge_rc_data: - # First compute_rockingcurve() for output arrays sizing - dout = _rockingcurve.compute_rockingcurve( - crystal=crystal, - din=din, - lamb=lamb[0]*1e10, - miscut=miscut, - therm_exp=therm_exp, - temp_limits=temp_limits, - plot_therm_exp=plot_rcs, - alpha_limits=alpha_limits, - nn=None, - plot_asf=False, - plot_power_ratio=plot_rcs, - plot_asymmetry=False, - plot_cmaps=False, - returnas=dict, - ) - - TD = np.zeros((na,), dtype=float) - if therm_exp: - TD = dout['Temperature changes (°C)'] - nT = TD.size - angles = np.zeros((na,), dtype=float) - if miscut: - angles = dout['Miscut angles (deg)'] - nangles = angles.size - power_ratio = np.full(( - nlamb, - dout['Power ratio'].shape[0], - dout['Power ratio'].shape[1], - dout['Power ratio'].shape[2], - dout['Power ratio'].shape[3], - ), np.nan) - power_ratio[ll, ...] = dout['Power ratio'] - - id_alpha0 = find_nearest(angles, alpha0) - id_temp0 = find_nearest(TD, temp0) - - if miscut and therm_exp: - dth = dout['Glancing angles'][0, id_temp0, id_alpha0, :] - ndth = dth.size - elif not miscut and not therm_exp: - dth = dout['Glancing angles'][0, 0, 0, :] - ndth = dth.size - elif miscut and not therm_exp: - dth = dout['Glancing angles'][0, 0, id_alpha0, :] - ndth = dth.size - elif not miscut and therm_exp: - dth = dout['Glancing angles'][0, id_temp0, 0, :] - ndth = dth.size - - # For each wavelength, get results dictionnary of the associated - # diffraction pattern - power_ratio = np.full(( - nlamb, - dout['Power ratio'].shape[0], - dout['Power ratio'].shape[1], - dout['Power ratio'].shape[2], - dout['Power ratio'].shape[3], - ), np.nan) - xi_rc = np.full((nlamb, ndth, nphi), np.nan) - xj_rc = xi_rc.copy() - xi_atprmax = np.full((nlamb, 1), np.nan) - xj_atprmax = xi_atprmax.copy() - bragg_atprmax = xi_atprmax.copy() - lamb_atprmax = xi_atprmax.copy() - - pix_horiz = np.linspace( - det['outline'][0, 0], - det['outline'][0, 1], - nxi - ) - pix_verti = np.linspace( - det['outline'][1, 1], - det['outline'][1, 2], - nxj - ) - # data = np.full((1, pix_verti.size, pix_horiz.size), 0) - data = np.empty([1, pix_verti.size, pix_horiz.size]) + lxi_rc, lxj_rc = [], [] + lxi_atprmax = [] + lbragg_atprmax = [] + llamb_atprmax = [] - for ll in range(nlamb): + for ii in range(len(lself)): + # First compute_rockingcurve() for output arrays sizing dout = _rockingcurve.compute_rockingcurve( crystal=crystal, din=din, - lamb=lamb[ll]*1e10, + lamb=lamb[0]*1e10, miscut=miscut, therm_exp=therm_exp, temp_limits=temp_limits, @@ -2509,8 +2470,10 @@ def find_nearest(array, value): plot_power_ratio=plot_rcs, plot_asymmetry=False, plot_cmaps=False, + plot_as_wavel=plot_as_wavel, returnas=dict, ) + TD = np.zeros((na,), dtype=float) if therm_exp: TD = dout['Temperature changes (°C)'] @@ -2519,6 +2482,13 @@ def find_nearest(array, value): if miscut: angles = dout['Miscut angles (deg)'] nangles = angles.size + power_ratio = np.full(( + nlamb, + dout['Power ratio'].shape[0], + dout['Power ratio'].shape[1], + dout['Power ratio'].shape[2], + dout['Power ratio'].shape[3], + ), np.nan) power_ratio[ll, ...] = dout['Power ratio'] id_alpha0 = find_nearest(angles, alpha0) @@ -2527,91 +2497,166 @@ def find_nearest(array, value): if miscut and therm_exp: dth = dout['Glancing angles'][0, id_temp0, id_alpha0, :] ndth = dth.size - ind_pr_max = np.where( - power_ratio[ll, 0, id_temp0, id_alpha0] == np.max( - power_ratio[ll, 0, id_temp0, id_alpha0] - ) - ) - dth_atprmax = dth[ind_pr_max] elif not miscut and not therm_exp: dth = dout['Glancing angles'][0, 0, 0, :] ndth = dth.size - ind_pr_max = np.where( - power_ratio[ll, 0, 0, 0] == np.max( - power_ratio[ll, 0, 0, 0] - ) - ) - dth_atprmax = dth[ind_pr_max] elif miscut and not therm_exp: dth = dout['Glancing angles'][0, 0, id_alpha0, :] ndth = dth.size - ind_pr_max = np.where( - power_ratio[ll, 0, 0, id_alpha0] == np.max( - power_ratio[ll, 0, 0, id_alpha0] - ) - ) - dth_atprmax = dth[ind_pr_max] elif not miscut and therm_exp: dth = dout['Glancing angles'][0, id_temp0, 0, :] ndth = dth.size - ind_pr_max = np.where( - power_ratio[ll, 0, id_temp0, 0] == np.max( - power_ratio[ll, 0, id_temp0, 0] - ) - ) - dth_atprmax = dth[ind_pr_max] - # Compute wavelength arcs for each glancing angle to obtain - # the shadow of the diffraction pattern on the detector - for mm in range(ndth): - ( - xi_rc[ll, mm, :], xj_rc[ll, mm, :], - ) = self.calc_xixj_from_braggphi( - bragg=np.full(phi.shape, dth[mm]), - phi=phi, - dtheta=0., - psi=0., - n=n, - det=det, + # For each wavelength, get results dictionnary of the associated + # diffraction pattern + power_ratio = np.full(( + nlamb, + dout['Power ratio'].shape[0], + dout['Power ratio'].shape[1], + dout['Power ratio'].shape[2], + dout['Power ratio'].shape[3], + ), np.nan) + xi_rc = np.full((nlamb, ndth, nphi), np.nan) + xj_rc = xi_rc.copy() + xi_atprmax = np.full((nlamb, 1), np.nan) + xj_atprmax = xi_atprmax.copy() + bragg_atprmax = xi_atprmax.copy() + lamb_atprmax = xi_atprmax.copy() + + pix_horiz = np.linspace( + det['outline'][0, 0], + det['outline'][0, 1], + nxi + ) + pix_verti = np.linspace( + det['outline'][1, 1], + det['outline'][1, 2], + nxj + ) + if ii == 0: + data = np.empty([1, pix_verti.size, pix_horiz.size]) + else: + None + + for ll in range(nlamb): + dout = _rockingcurve.compute_rockingcurve( + crystal=crystal, + din=din, + lamb=lamb[ll]*1e10, miscut=miscut, - strict=strict, - plot=False, + therm_exp=therm_exp, + temp_limits=temp_limits, + plot_therm_exp=plot_rcs, + alpha_limits=alpha_limits, + nn=None, + plot_as_wavel=plot_as_wavel, + plot_asf=False, + plot_power_ratio=plot_rcs, + plot_asymmetry=False, + plot_cmaps=False, + returnas=dict, ) - if plot_simu_image: - for nn in range(phi.size): - if np.isfinite(xi_rc[ll, mm, nn]): - idxi = np.nanargmin( - np.abs( - xi_rc[ll, mm, nn] - pix_horiz + TD = np.zeros((na,), dtype=float) + if therm_exp: + TD = dout['Temperature changes (°C)'] + nT = TD.size + angles = np.zeros((na,), dtype=float) + if miscut: + angles = dout['Miscut angles (deg)'] + nangles = angles.size + power_ratio[ll, ...] = dout['Power ratio'] + + id_alpha0 = find_nearest(angles, alpha0) + id_temp0 = find_nearest(TD, temp0) + + if miscut and therm_exp: + dth = dout['Glancing angles'][0, id_temp0, id_alpha0, :] + ndth = dth.size + ind_pr_max = np.where( + power_ratio[ll, 0, id_temp0, id_alpha0] == np.max( + power_ratio[ll, 0, id_temp0, id_alpha0] + ) + ) + dth_atprmax = dth[ind_pr_max] + elif not miscut and not therm_exp: + dth = dout['Glancing angles'][0, 0, 0, :] + ndth = dth.size + ind_pr_max = np.where( + power_ratio[ll, 0, 0, 0] == np.max( + power_ratio[ll, 0, 0, 0] + ) + ) + dth_atprmax = dth[ind_pr_max] + elif miscut and not therm_exp: + dth = dout['Glancing angles'][0, 0, id_alpha0, :] + ndth = dth.size + ind_pr_max = np.where( + power_ratio[ll, 0, 0, id_alpha0] == np.max( + power_ratio[ll, 0, 0, id_alpha0] + ) + ) + dth_atprmax = dth[ind_pr_max] + elif not miscut and therm_exp: + dth = dout['Glancing angles'][0, id_temp0, 0, :] + ndth = dth.size + ind_pr_max = np.where( + power_ratio[ll, 0, id_temp0, 0] == np.max( + power_ratio[ll, 0, id_temp0, 0] + ) + ) + dth_atprmax = dth[ind_pr_max] + + # Compute wavelength arcs for each glancing angle to obtain + # the shadow of the diffraction pattern on the detector + for mm in range(ndth): + ( + xi_rc[ll, mm, :], xj_rc[ll, mm, :], + ) = lself[ii].calc_xixj_from_braggphi( + bragg=np.full(phi.shape, dth[mm]), + phi=phi, + dtheta=0., + psi=0., + n=n, + det=det, + miscut=miscut, + strict=strict, + plot=False, + ) + if plot_simu_image: + for nn in range(phi.size): + if np.isfinite(xi_rc[ll, mm, nn]): + idxi = np.nanargmin( + np.abs( + xi_rc[ll, mm, nn] - pix_horiz + ) ) - ) - idxj = np.nanargmin( - np.abs( - xj_rc[ll, mm, nn] - pix_verti + idxj = np.nanargmin( + np.abs( + xj_rc[ll, mm, nn] - pix_verti + ) ) - ) - pr1 = power_ratio[ll, 0, 0, 0, mm] - pr2 = power_ratio[ll, 1, 0, 0, mm] - val = pr1 + pr2 - data[0, idxj, idxi] += val - - xi_atprmax[ll] = xi_rc[ll, ind_pr_max, int(nphi/2)] - xj_atprmax[ll] = xj_rc[ll, ind_pr_max, int(nphi/2)] - # self.update_miscut(alpha=0., beta=0.) - # if therm_exp: - # self.dmat['d'] = d_atom[nn]*1e-10 - # else: - # self.dmat['d'] = d_atom[0]*1e-10 - ( - bragg_atprmax[ll], _, lamb_atprmax[ll], - ) = self.get_lambbraggphi_from_ptsxixj_dthetapsi( - xi=xi_atprmax[ll], xj=xj_atprmax[ll], det=det, - dtheta=0, psi=0, - miscut=miscut, - n=n, - grid=True, - return_lamb=True, - ) + pr1 = power_ratio[ll, 0, 0, 0, mm] + pr2 = power_ratio[ll, 1, 0, 0, mm] + val = (pr1 + pr2) + data[0, idxj, idxi] += val + + xi_atprmax[ll] = xi_rc[ll, ind_pr_max, int(nphi/2)] + xj_atprmax[ll] = xj_rc[ll, ind_pr_max, int(nphi/2)] + ( + bragg_atprmax[ll], _, lamb_atprmax[ll], + ) = lself[ii].get_lambbraggphi_from_ptsxixj_dthetapsi( + xi=xi_atprmax[ll], xj=xj_atprmax[ll], det=det, + dtheta=0, psi=0, + miscut=miscut, + n=n, + grid=True, + return_lamb=True, + ) + lxi_rc.append(xi_rc) + lxj_rc.append(xj_rc) + lxi_atprmax.append(xi_atprmax) + llamb_atprmax.append(lamb_atprmax) + lbragg_atprmax.append(bragg_atprmax) else: power_ratio = None dth = None @@ -2626,47 +2671,58 @@ def find_nearest(array, value): angles = None data = None - # Reset parameters as at beginning - if miscut: - self.update_miscut(alpha=alpha0, beta=0.) - else: - self.update_miscut(alpha=0., beta=0.) - if therm_exp: - self.dmat['d'] = d_atom[id_temp0]*1e-10 - else: - self.dmat['d'] = d_atom[0]*1e-10 + for ii in range(len(lself)): + if miscut: + lself[ii].update_miscut(alpha=alpha0, beta=0.) + else: + lself[ii].update_miscut(alpha=0., beta=0.) + if therm_exp: + lself[ii].dmat['d'] = d_atom[id_temp0]*1e-10 + else: + lself[ii].dmat['d'] = d_atom[0]*1e-10 if mode == 'raw det': - xi = xi - det['outline'][0, 0] - xj = xj - det['outline'][1, 0] - xi_rc = xi_rc - det['outline'][0, 0] - xj_rc = xj_rc - det['outline'][1, 0] - det['outline'][0] = det['outline'][0] - det['outline'][0, 0] - det['outline'][1] = det['outline'][1] - det['outline'][1, 0] + for ii in range(len(lself)): + lxi[ii] = lxi[ii] - det['outline'][0, 0] + lxj[ii] = lxj[ii] - det['outline'][1, 0] + lxi_rc[ii] = lxi_rc[ii] - det['outline'][0, 0] + lxj_rc[ii] = lxj_rc[ii] - det['outline'][1, 0] + det['outline'][0] = det['outline'][0] - det['outline'][0, 0] + det['outline'][1] = det['outline'][1] - det['outline'][1, 0] if plot_simu_image: fig0 = plt.figure() ax0 = fig0.add_subplot() - ax0.set_xlabel(r'Pixel coordinate', fontsize=15) - ax0.set_ylabel(r'Pixel coordinate', fontsize=15) - ax0.set_title(r'Simulated 2D spectra', fontsize=15) - ax0.imshow( + ax0.set_xlabel(r'$\#$ pixel', fontsize=fs) + ax0.set_ylabel(r'$\#$ pixel', fontsize=fs) + #ax0.set_title(r'Simulated 2D spectra', fontsize=fs) + ax0.tick_params(labelsize=fs) + im = ax0.imshow( data[0, :, :], origin='lower', + cmap='viridis', interpolation='nearest', aspect='auto', ) + cbar = plt.colorbar( + im, + orientation='vertical', + ax=ax0, + ) + cbar.ax.tick_params( + labelsize=fs + ) dout = { 'lamb': lamb, - 'xi': xi, - 'xj': xj, - 'xi_rc': xi_rc, - 'xj_rc': xj_rc, - 'xi_atprmax': xi_atprmax, - 'lamb_atprmax': lamb_atprmax, - 'bragg_atprmax': bragg_atprmax, + 'xi': lxi, + 'xj': lxj, + 'xi_rc': lxi_rc, + 'xj_rc': lxj_rc, + 'xi_atprmax': lxi_atprmax, + 'lamb_atprmax': llamb_atprmax, + 'bragg_atprmax': lbragg_atprmax, 'data': data, } @@ -2697,53 +2753,54 @@ def find_nearest(array, value): # None if plot_line_tracing: - ax = _plot_optics.CrystalBragg_plot_line_tracing_on_det( - # ------------------------------ - # basic - cryst=self, - dcryst=dcryst, - lamb=lamb, - dlamb=dlamb, - xi=xi, - xj=xj, - xi_er=xi_er, - xj_er=xj_er, - # ----------------------------- - # w/ rocking curves data - merge_rc_data=merge_rc_data, - power_ratio=power_ratio, - dth=dth, - ndth=ndth, - nn=nn, - xi_rc=xi_rc, - xj_rc=xj_rc, - xi_atprmax=xi_atprmax, - bragg_atprmax=bragg_atprmax, - lamb_atprmax=lamb_atprmax, - TD=TD, - angles=angles, - # ----------------------------- - # w/ miscut and/or temp changes - alpha0=alpha0, - temp0=temp0, - id_temp0=id_temp0, - johann=johann, - rocking=rocking, - miscut=miscut, - therm_exp=therm_exp, - det=det, - # ---------------------------- - # plot parameters - ax=ax, - dleg=dleg, - color=color, - fs=fs, - dmargin=dmargin, - wintit=wintit, - tit=tit, - plot_perfect=plot_perfect, - ) - return ax, dout + for ii in range(len(lself)): + ax = _plot_optics.CrystalBragg_plot_line_tracing_on_det( + # ------------------------------ + # basic + cryst=lself[ii], + dcryst=dcryst, + lamb=lamb, + dlamb=dlamb, + xi=lxi[ii], + xj=lxj[ii], + xi_er=xi_er, + xj_er=xj_er, + # ----------------------------- + # w/ rocking curves data + merge_rc_data=merge_rc_data, + power_ratio=power_ratio, + dth=dth, + ndth=ndth, + nn=nn, + xi_rc=lxi_rc[ii], + xj_rc=lxj_rc[ii], + xi_atprmax=lxi_atprmax[ii], + bragg_atprmax=lbragg_atprmax[ii], + lamb_atprmax=llamb_atprmax[ii], + TD=TD, + angles=angles, + # ----------------------------- + # w/ miscut and/or temp changes + alpha0=alpha0, + temp0=temp0, + id_temp0=id_temp0, + johann=johann, + rocking=rocking, + miscut=miscut, + therm_exp=therm_exp, + det=det, + # ---------------------------- + # plot parameters + ax=ax, + dleg=dleg, + color=color, + fs=fs, + dmargin=dmargin, + wintit=wintit, + tit=tit, + plot_perfect=plot_perfect, + ) + return ax, dout else: return dout From 9c3e1fae4535d6904ff69afc2c9f8e4f09d33c6f Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 17 Jan 2024 09:37:00 +0100 Subject: [PATCH 7/9] [#870] plot_optics: minor plotting option changes --- tofu/geom/_plot_optics.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/tofu/geom/_plot_optics.py b/tofu/geom/_plot_optics.py index 82306bf74..18bb09b4b 100644 --- a/tofu/geom/_plot_optics.py +++ b/tofu/geom/_plot_optics.py @@ -996,7 +996,7 @@ def CrystalBragg_plot_line_tracing_on_det( if color is None: color = 'k' if fs is None: - fs = (8, 8) + fs = 15 if dmargin is None: dmargin = { 'top': 0.950, @@ -1032,15 +1032,16 @@ def CrystalBragg_plot_line_tracing_on_det( # Plot if ax is None: - fig = plt.figure(figsize=fs) + fig = plt.figure() gs = gridspec.GridSpec(1, 1, **dmargin) ax = fig.add_subplot(gs[0, 0], aspect='equal', adjustable='datalim') if wintit is not False: fig.canvas.manager.set_window_title(wintit) if tit is not False: fig.suptitle(tit, size=14, weight='bold') - ax.set_xlabel(r'Pixel coordinate $x_{i}$ [m]', fontsize=15) - ax.set_ylabel(r'Pixel coordinate $x_{j}$ [m]', fontsize=15) + ax.set_xlabel(r'Pixel coordinate $x_{i}$ [m]', fontsize=fs) + ax.set_ylabel(r'Pixel coordinate $x_{j}$ [m]', fontsize=fs) + ax.tick_params(labelsize=fs) ax.set_xlim( det['outline'][0, :].min() - 0.01, @@ -1191,7 +1192,7 @@ def CrystalBragg_plot_angular_shift_on_det_tracing( # Plot # ------------ - fig = plt.figure(figsize=fs) + fig = plt.figure() gs = gridspec.GridSpec(1, 3) # , **dmargin) ax0 = fig.add_subplot(gs[0, 0], aspect='equal', adjustable='datalim') ax0.set_title('Pixel offset [m]', fontsize=20) @@ -1313,7 +1314,7 @@ def CrystalBragg_plot_johannerror( # Plot # ------------ - fig = plt.figure(figsize=fs) + fig = plt.figure() gs = gridspec.GridSpec(1, 3, **dmargin) ax0 = fig.add_subplot(gs[0, 0], aspect='equal') # adjustable='datalim') ax1 = fig.add_subplot( From d3a7ca9bba762a9291b5762af4f204728f26360a Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 17 Jan 2024 09:39:40 +0100 Subject: [PATCH 8/9] [#870] fit12d: remove pdb option --- tofu/spectro/_fit12d.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tofu/spectro/_fit12d.py b/tofu/spectro/_fit12d.py index 3570cedfc..12459f71c 100644 --- a/tofu/spectro/_fit12d.py +++ b/tofu/spectro/_fit12d.py @@ -379,7 +379,6 @@ def reshape_custom(aa, nspect0=nspect0): if dinput['double'] is True: x2[:, dind['dratio']['x']] = sol_x[:, dind['dratio']['x']] x2[:, dind['dshift']['x']] = sol_x[:, dind['dshift']['x']] - import pdb; pdb.set_trace() # DB sol_x = x2 # --------------------------- From 7de2959ac22211998def144f53b5970a055a1516 Mon Sep 17 00:00:00 2001 From: AD265925 Date: Wed, 17 Jan 2024 09:41:37 +0100 Subject: [PATCH 9/9] [#870] _rockingcurve: new modifs about RCs creation, plotting vs wavelength, plot parameter supdated, --- tofu/spectro/_rockingcurve.py | 140 ++++++++++++++++++++++++++-------- 1 file changed, 109 insertions(+), 31 deletions(-) diff --git a/tofu/spectro/_rockingcurve.py b/tofu/spectro/_rockingcurve.py index 657cc405d..b02189b7c 100644 --- a/tofu/spectro/_rockingcurve.py +++ b/tofu/spectro/_rockingcurve.py @@ -13,6 +13,7 @@ # tofu +import tofu as tf from . import _rockingcurve_def as _def @@ -36,6 +37,8 @@ def compute_rockingcurve( therm_exp=None, temp_limits=None, # Plot + plot_as_wavel=None, + dcryst=None, plot_therm_exp=None, plot_asf=None, plot_power_ratio=None, @@ -43,6 +46,7 @@ def compute_rockingcurve( plot_cmaps=None, # Returning dictionnary returnas=None, + fs=None, ): """ The code evaluates, for a given wavelength and Miller indices set, the inter-plane distance d, the Bragg angle of reference and the complex @@ -131,6 +135,7 @@ def compute_rockingcurve( # ------------ ( + dcryst, crystal, din, lamb, # lattice expansion therm_exp, @@ -140,6 +145,7 @@ def compute_rockingcurve( nn, na, # plotting + plot_as_wavel, plot_asf, plot_therm_exp, plot_power_ratio, @@ -147,6 +153,7 @@ def compute_rockingcurve( plot_cmaps, # return returnas, + fs, ) = _checks(**locals()) # Classical electronical radius, in Angstroms, from the NIST Reference on @@ -311,7 +318,7 @@ def compute_rockingcurve( if miscut is False and therm_exp is False: ( alpha, bb, polar, g, y, power_ratio, max_pr, th, dth, - rhg, P_per, P_mos, P_dyn, det_perp, det_para, + rhg, P_per, P_mos, P_dyn, det_perp, det_para, det_tot, ) = CrystBragg_comp_integrated_reflect( lamb=lamb, re=re, Volume=Volume, Zo=din['atoms_Z'][-1], theta=theta, mu=mu, F_re=F_re, psi_re=psi_re, psi0_dre=psi0_dre, psi0_im=psi0_im, @@ -326,7 +333,7 @@ def compute_rockingcurve( alpha, bb, polar, g, y, power_ratio, max_pr, th, dth, rhg, rhg_perp, rhg_para, rhg_perp_norm, rhg_para_norm, P_per, P_mos, P_dyn, - det_perp, det_para, det_perp_norm, det_para_norm, + det_perp, det_para, det_tot, det_perp_norm, det_para_norm, shift_perp, shift_para, ) = CrystBragg_comp_integrated_reflect( lamb=lamb, re=re, Volume=Volume, Zo=din['atoms_Z'][-1], theta=theta, mu=mu, @@ -351,6 +358,7 @@ def compute_rockingcurve( if plot_power_ratio: CrystalBragg_plot_power_ratio_vs_glancing_angle( + dcryst=dcryst, din=din, lamb=lamb, alpha_limits=alpha_limits, theta=theta, theta_deg=theta_deg, @@ -358,6 +366,8 @@ def compute_rockingcurve( bb=bb, polar=polar, alpha=alpha, miscut=miscut, na=na, nn=nn, therm_exp=therm_exp, T0=T0, TD=TD, + plot_as_wavel=plot_as_wavel, + fs=fs, ) # Plot integrated reflect., asymmetry angle & RC width vs glancing angle @@ -424,6 +434,7 @@ def compute_rockingcurve( 'Maximum reflectivity (para. compo)': max_pr[1], 'RC width (perp. compo)': det_perp, 'RC width (para. compo)': det_para, + 'RC width (both compo)': det_tot, # polar 'polar': polar, # y => angles @@ -459,6 +470,7 @@ def compute_rockingcurve( def _checks( # Type of crystal + dcryst=None, crystal=None, din=None, # Wavelength @@ -470,6 +482,7 @@ def _checks( therm_exp=None, temp_limits=None, # Plot + plot_as_wavel=None, plot_therm_exp=None, plot_asf=None, plot_power_ratio=None, @@ -477,6 +490,7 @@ def _checks( plot_cmaps=None, # Returning dictionnary returnas=None, + fs=None, ): # ------------ @@ -510,6 +524,15 @@ def _checks( # -------------- # plotting args + if plot_as_wavel is None: + plot_as_wavel = True + + if plot_as_wavel and dcryst is None: + msg = ( + "Please load the CrystalBragg object !" + ) + raise Exception(msg) + if plot_asf is None: plot_asf = False @@ -532,7 +555,11 @@ def _checks( if returnas is None: returnas = dict + if fs is None: + fs = 15 + return ( + dcryst, crystal, din, lamb, # lattice expansion therm_exp, @@ -542,6 +569,7 @@ def _checks( nn, na, # plotting + plot_as_wavel, plot_asf, plot_therm_exp, plot_power_ratio, @@ -549,6 +577,7 @@ def _checks( plot_cmaps, # return returnas, + fs, ) @@ -1225,11 +1254,11 @@ def half_max_x(x, y): P_dyn = np.full((theta.size, alpha.size), np.nan) ( - rhg_perp, rhg_para, pat_cent_perp, pat_cent_para, - det_perp, det_para, shift_perp, shift_para, + rhg_perp, rhg_para, pat_cent_perp, pat_cent_para, pat_cent_tot, + det_perp, det_para, det_tot, shift_perp, shift_para, ) = ( P_dyn.copy(), P_dyn.copy(), P_dyn.copy(), P_dyn.copy(), P_dyn.copy(), - P_dyn.copy(), P_dyn.copy(), P_dyn.copy(), + P_dyn.copy(), P_dyn.copy(), P_dyn.copy(), P_dyn.copy(), P_dyn.copy(), ) for i in range(theta.size): @@ -1265,12 +1294,16 @@ def half_max_x(x, y): # Coordinates of full width at mid high sides of FWHM hmx_perp = half_max_x(dth[0, i, j, :], power_ratio[0, i, j, :]) hmx_para = half_max_x(dth[1, i, j, :], power_ratio[1, i, j, :]) + pr_tot = power_ratio[0, i, j, :] + power_ratio[1, i, j, :] + hmx_tot = half_max_x(dth[0, i, j, :], pr_tot) # Center of FWMH pat_cent_perp[i, j] = (hmx_perp[1] + hmx_perp[0])/2. pat_cent_para[i, j] = (hmx_para[1] + hmx_para[0])/2. + pat_cent_tot[i, j] = (hmx_tot[1] + hmx_tot[0])/2. # Width of FWMH det_perp[i, j] = hmx_perp[1] - hmx_perp[0] det_para[i, j] = hmx_para[1] - hmx_para[0] + det_tot[i, j] = hmx_tot[1] - hmx_tot[0] # Normalization for DeltaT=0 & alpha=0 and # computation of the shift in glancing angle corresponding to @@ -1308,7 +1341,7 @@ def half_max_x(x, y): alpha, bb, polar, g, y, power_ratio, max_pr, th, dth, rhg, P_per, P_mos, P_dyn, - det_perp, det_para, + det_perp, det_para, det_tot, ) else: return ( @@ -1316,7 +1349,8 @@ def half_max_x(x, y): power_ratio, max_pr, th, dth, rhg, rhg_perp, rhg_para, rhg_perp_norm, rhg_para_norm, P_per, P_mos, P_dyn, - det_perp, det_para, det_perp_norm, det_para_norm, + det_perp, det_para, det_tot, + det_perp_norm, det_para_norm, shift_perp, shift_para, ) @@ -1398,6 +1432,7 @@ def CrystalBragg_plot_atomic_scattering_factor( def CrystalBragg_plot_power_ratio_vs_glancing_angle( + dcryst=None, # Lattice parameters din=None, lamb=None, # Lattice modifications @@ -1408,6 +1443,8 @@ 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, + plot_as_wavel=None, + fs=None, ): """ All plots of rocking curve is done, not with respect to the glancing angle (theta - thetaBragg) where thetaBragg may vary if the temperature @@ -1451,6 +1488,10 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( col2 = {'black': dd2[0]} name = din['name'] + if 'Quartz' in name: + name2 = 'Quartz' + elif 'Ge' in name: + name2 = 'Ge' miller = np.r_[ int(din['miller'][0]), int(din['miller'][1]), @@ -1467,65 +1508,85 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( 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} $\AA$', fontsize=15, + f'{name2}, ' + f'Miller indices: ({miller[0]},{miller[1]},{miller[2]})' + + fr', $\lambda$={lamb[0]} $\AA$', + fontsize=fs, ) - ax.set_xlabel(r'Diffracting angle $\theta$ (rad)', fontsize=15) - ax.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=15) + if plot_as_wavel: + ax.set_xlabel(r'Wavelength $\lambda$ [$\AA$]', fontsize=fs) + else: + ax.set_xlabel(r'Diffracting angle $\theta$ [rad]', fontsize=fs) + ax.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=fs) + ax.set_ylim(-0.05, 1.) + ax.tick_params(labelsize=fs) + ax.grid() 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]) + ax00.tick_params(labelsize=fs) + ax01 = fig1.add_subplot(gs[0, 1]) + ax01.tick_params(labelsize=fs) ax02 = fig1.add_subplot(gs[0, 2]) + ax02.tick_params(labelsize=fs) ax11 = fig1.add_subplot(gs[1, 1]) + ax11.tick_params(labelsize=fs) ax10 = fig1.add_subplot(gs[1, 0]) + ax10.tick_params(labelsize=fs) ax12 = fig1.add_subplot(gs[1, 2]) + ax12.tick_params(labelsize=fs) ax21 = fig1.add_subplot(gs[2, 1]) + ax21.tick_params(labelsize=fs) ax20 = fig1.add_subplot(gs[2, 0]) + ax20.tick_params(labelsize=fs) ax22 = fig1.add_subplot(gs[2, 2]) + ax22.tick_params(labelsize=fs) fig1.suptitle( - f'{name}, ' + f'({miller[0]},{miller[1]},{miller[2]})' + - fr', $\lambda$={lamb} $\AA$' + + f'{name2}, ' + f'({miller[0]},{miller[1]},{miller[2]})' + + fr', $\lambda$={lamb[0]} $\AA$' + r', $\theta_{B}$=' + fr'{np.round(theta[nn], 5)} rad', - fontsize=15, + fontsize=fs, ) ax00.set_title( r'$\alpha$=({}) arcmin'.format(np.round(alpha_deg[0]*60, 3)), - fontsize=15 + fontsize=fs ) ax01.set_title( r'$\alpha$=({}) arcmin'.format(np.round(alpha_deg[nn]*60, 3)), - fontsize=15 + fontsize=fs ) ax02.set_title( r'$\alpha$=({}) arcmin'.format(np.round(alpha_deg[na-1]*60, 3)), - fontsize=15 + fontsize=fs ) ax022 = ax02.twinx() ax022.set_ylabel( - r'$\Delta T$=({})°C'.format(TD[0]), fontsize=15 + r'$\Delta T$=({})°C'.format(TD[0]), fontsize=fs ) + ax022.tick_params(labelsize=fs) ax122 = ax12.twinx() ax122.set_ylabel( - r'$\Delta T$=({})°C'.format(TD[nn]), fontsize=15 + r'$\Delta T$=({})°C'.format(TD[nn]), fontsize=fs ) + ax122.tick_params(labelsize=fs) ax222 = ax22.twinx() ax222.set_ylabel( - r'$\Delta T$=({})°C'.format(TD[na-1]), fontsize=15 + r'$\Delta T$=({})°C'.format(TD[na-1]), fontsize=fs ) - 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) + ax222.tick_params(labelsize=fs) + ax20.set_xlabel(r'$\theta$ (rad)', fontsize=fs) + ax21.set_xlabel(r'$\theta$ (rad)', fontsize=fs) + ax22.set_xlabel(r'$\theta$ (rad)', fontsize=fs) + ax00.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=fs) + ax10.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=fs) + ax20.set_ylabel('Power ratio P$_H$/P$_0$', fontsize=fs) # Plot # ---- lc = [miscut is True, miscut is False] + """ if not therm_exp and any(lc): for j in range(na): if any(j == dd): @@ -1577,19 +1638,36 @@ def CrystalBragg_plot_power_ratio_vs_glancing_angle( # Plot the sum of both polarizations lc = [miscut is True, miscut is False] if not therm_exp and any(lc): + if plot_as_wavel: + dth2 = tf.geom.CrystalBragg.get_lamb_from_bragg( + self=dcryst, + bragg=dth[0, 0, 0, :] + ) + x0 = dth2 + ll = tf.geom.CrystalBragg.get_lamb_from_bragg( + self=dcryst, + bragg=theta, + ) + x1 = ll + else: + x0 = dth[0, 0, 0, :] + x1 = theta + ax.plot( - dth[0, 0, 0, :], + # dth[0, 0, 0, :], + x0 * 1e10, power_ratio[0, 0, 0] + power_ratio[1, 0, 0], '-', c='black', ) ax.axvline( - theta, color='black', linestyle='-.', + # theta, + x1 * 1e10, + color='black', linestyle='-.', label=r'$\theta_B$= {} rad'.format( - np.round(theta, 6) + np.round(x1, 6) ), ) - """ if not miscut and therm_exp is True: colors = ['blue', 'black', 'red']