diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index 957984a71..1f7a89343 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -1946,29 +1946,59 @@ 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, + split=None, + direction=None, + nb=None, + which=None, + # ----------------------- # Wavelength lamb=None, + subarea=None, + shift=None, + v=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_as_wavel=None, 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, + save=None, ): """ Visualize the de-focusing by ray-tracing of chosen lamb Possibility to plot few wavelength' arcs on the same plot. @@ -1998,7 +2028,106 @@ def plot_line_on_det_tracing( By default to 10°C so temp0=35 """ + # ----------------------- # 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'] + lself = [] + 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]: + lself.append(self1) + + elif direction == ldir[0] and which == lwhich[2]: + lself.append(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: + lself.append(self) + + # -------------------- + # other generic checks + lok = [ k0 for k0 in _rockingcurve_def._DCRYST.keys() if 'xxx' not in k0.lower() @@ -2010,20 +2139,96 @@ 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 = {} + lsubarea = [ + 'ArXVII-simul', 'ArXVII-woW', 'ArXVII-wxyz' + ] + if lamb is None and subarea 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 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_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 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_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 is None and subarea == lsubarea[2]: + 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) + 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: @@ -2033,9 +2238,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,91 +2251,163 @@ 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 plot_as_wavel is None: + plot_as_wavel = False + + 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 + 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.) - # 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, - 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 + 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=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 - # Get reference ray-tracing - bragg = self._checkformat_bragglamb(lamb=lamb, n=n) - if nphi2 is None: - nphi2 = 50 - nphi = 2*nphi2 - 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, + 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, - strict=strict, - plot=False, + 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 + + 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: if lpsi is None: @@ -2156,39 +2436,44 @@ 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) - xj_rc = xi_rc.copy() - power_ratio = xi_rc.copy() - xi_atprmax = xi_rc.copy() - xj_atprmax = xi_rc.copy() - bragg_atprmax = xi_atprmax.copy() - lamb_atprmax = xi_atprmax.copy() - - # For each wavelength, get results dictionnary of the associated - # diffraction pattern - for ll in range(nlamb): + lxi_rc, lxj_rc = [], [] + lxi_atprmax = [] + lbragg_atprmax = [] + llamb_atprmax = [] + + 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, + 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, + alpha_limits=alpha_limits, + nn=None, + plot_asf=False, + 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)'] @@ -2197,161 +2482,326 @@ def find_nearest(array, value): if miscut: angles = dout['Miscut angles (deg)'] nangles = angles.size - power_ratio = np.resize(power_ratio, ( + 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'] - def find_nearest(array, value): - array = np.asarray(array) - idx = (np.abs(array - value)).argmin() - return idx - id_alpha0 = find_nearest(angles, alpha0) + id_temp0 = find_nearest(TD, temp0) - # 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 - 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] - # Resize results arrays - xi_rc = np.resize(xi_rc, (nlamb, ndth, nphi)) + # 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.resize(xi_atprmax, (nlamb, 1)) + xi_atprmax = np.full((nlamb, 1), np.nan) 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): + 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, + 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, + ) + 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 + ) + ) + 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)] ( - 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, + 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, - strict=strict, - plot=False, + n=n, + grid=True, + return_lamb=True, ) - 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: - 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, - ) + 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 + 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: - 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': + 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', 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 + ) - # Plot - if plot: - if merge_rc_data: - return _plot_optics.CrystalBragg_plot_line_tracing_on_det( - cryst=self, dcryst=dcryst, + dout = { + 'lamb': lamb, + '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, + } + + # save + if save: + 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: + for ii in range(len(lself)): + ax = _plot_optics.CrystalBragg_plot_line_tracing_on_det( + # ------------------------------ + # basic + cryst=lself[ii], + 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, + 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, - 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, + 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, + 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, + # ---------------------------- + # 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 +3304,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, ) diff --git a/tofu/geom/_plot_optics.py b/tofu/geom/_plot_optics.py index f501d4bc8..18bb09b4b 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) + fs = 15 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,22 +1017,32 @@ 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) + 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, det['outline'][0, :].max() + 0.01, @@ -1017,11 +1051,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 +1066,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, ) @@ -1114,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) @@ -1236,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( 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 # --------------------------- 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'], 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']