diff --git a/tofu/physics_tools/electrons/emission/__init__.py b/tofu/physics_tools/electrons/emission/__init__.py index 5deac83c0..324733945 100644 --- a/tofu/physics_tools/electrons/emission/__init__.py +++ b/tofu/physics_tools/electrons/emission/__init__.py @@ -13,4 +13,5 @@ from ._xray_thin_target_integrated_plot import plot_xray_thin_d2cross_ei_anisotropy from ._xray_thin_target_integrated_dist import get_xray_thin_integ_dist from ._xray_thin_target_integrated_d2crossphi import get_d2cross_phi +from ._xray_thin_target_integrated_dist_responsivity_plot import plot_xray_thin_integ_dist_filter_anisotropy from ._xray_thin_target_integrated_dist_plot import plot_xray_thin_integ_dist diff --git a/tofu/physics_tools/electrons/emission/_xray_thin_target.py b/tofu/physics_tools/electrons/emission/_xray_thin_target.py index 155497344..215cea37a 100644 --- a/tofu/physics_tools/electrons/emission/_xray_thin_target.py +++ b/tofu/physics_tools/electrons/emission/_xray_thin_target.py @@ -2,6 +2,7 @@ import os import warnings +from typing import Any, Optional # Dict import numpy as np @@ -12,6 +13,9 @@ import datastock as ds +TupleDict = tuple[dict] + + # #################################################### # #################################################### # DEFAULT @@ -21,6 +25,10 @@ _PATH_HERE = os.path.dirname(__file__) +# version +_VERSION = 'EH' # most accurate, but slow + + # #################################################### # #################################################### # Differential cross-section @@ -28,8 +36,9 @@ def get_xray_thin_d3cross_ei( - # inputs - Z=None, + # target ion charge + Z: Optional[int] = None, + # Energy E_e0_eV=None, E_e1_eV=None, # directions @@ -37,19 +46,19 @@ def get_xray_thin_d3cross_ei( theta_e=None, dphi=None, # hypergeometric parameter - ninf=None, - source=None, + ninf: Optional[int] = None, + source: Optional[str] = None, # output customization - per_energy_unit=None, + per_energy_unit: Optional[str] = None, # version - version=None, + version: Optional[str] = None, # debug - debug=None, -): + debug: Optional[bool] = None, +) -> dict: """ Return a differential cross-section for thin-target bremsstrahlung - Allowws several formulas (version): - - 'BE': Elwert-Haug [1] + Allows several formulas (version): + - 'EH': Elwert-Haug [1] (default) . most general and accurate . Uses Sommerfield-Maue eigenfunctions . eq. (30) in [1] @@ -74,14 +83,15 @@ def get_xray_thin_d3cross_ei( Physics Reports, vol. 243, p. 317—353, 1994. - Inputs: + Inputs: (all angles in rad) E_e0_eV = kinetic energy of incident electron in eV E_e1_eV = kinetic energy of scattered electron in eV theta_e = (spherical) theta angle of scattered e vs incident e theta_ph = (spherical) theta angle of photon vs incident e phi_e = (spherical) phi angle of scattered e vs incident e phi_ph = (spherical) theta angle of photon vs incident e - (all angles in rad) + version = 'EH', 'BH' or 'BHE' + Limitations: - 'EH' implementation currently stalled because: @@ -354,7 +364,7 @@ def _check_cross( # ------------ if version is None: - version = 'EH' + version = _VERSION if isinstance(version, str): version = [version] @@ -891,7 +901,13 @@ def _angle_dependent_internediates( sintp = np.sin(theta_ph) cosdphi = np.cos(dphi) - cossindphi = costp*coste + sintp*sinte*cosdphi + costpcoste = costp * coste + sintpsintecosdphi = sintp * sinte * cosdphi + + cossindphi = costpcoste + sintpsintecosdphi + + sintecostp = sinte * costp + costesintp = coste * sintp # ------------- # Vectors / scalar product @@ -899,15 +915,15 @@ def _angle_dependent_internediates( # scalar sca_kp0 = kk * p0 * costp - sca_p01 = p1 * p0 * coste + sca_p01 = p0 * p1 * coste sca_kp1 = kk * p1 * cossindphi # vect{q} = vect{p0 - p1 - k} q2 = ( p0**2 + p1**2 + k2 - - 2.*p0*p1*coste - - 2.*p0*kk*costp - + 2.*p1*kk*cossindphi + - 2. * sca_p01 + - 2. * sca_kp0 + + 2. * sca_kp1 ) # ------------- @@ -945,9 +961,9 @@ def _angle_dependent_internediates( # + (sinte*sintp)**2 * sin(phi_p - phi_e)**2 # ) eta12 = p12 * ( - (sinte*costp)**2 - + (coste*sintp)**2 - - 2*(coste*costp)*(sinte*sintp*cosdphi) + sintecostp**2 + + costesintp**2 + - 2 * costpcoste * sintpsintecosdphi + (sinte*sintp)**2 * (1 - cosdphi**2) ) @@ -966,7 +982,7 @@ def _angle_dependent_internediates( # coste*sintp*(sinpp**2 + cospp**2) # - sinte*costp*(sinpe*sinpp + cospe*cospp) # ) - sca_eta01 = p0 * p1 * sintp * (coste*sintp - sinte*costp*cosdphi) + sca_eta01 = p0 * p1 * sintp * (costesintp - sintecostp * cosdphi) # --------------- # Intermediates 1 @@ -1326,8 +1342,8 @@ def _hyp2F1( bb=None, cc=None, zz=None, - ninf=None, - source=None, + ninf: Optional[int] = None, + source: Optional[str] = None, ): """ Hypergeometric function 2F1 with complex arguments @@ -1389,8 +1405,12 @@ def _hyp2F1( # Number of terms # ---------- - if ninf is None: - ninf = 50 + ninf = ds._generic_check._check_var( + ninf, 'ninf', + types=int, + default=50, + sign='>0', + ) nn = np.arange(0, ninf)[None, :] @@ -1495,11 +1515,11 @@ def _hyp2F1( def plot_xray_thin_d3cross_ei_vs_Literature( - version=None, - ninf=None, - source=None, - dax=None, -): + version: Optional[str] = None, + ninf: Optional[int] = None, + source: Optional[str] = None, + dax: Optional[dict[str, Any]] = None, +) -> TupleDict: """ Compare computed cross-sections vs literature values from Elwert-Haug Triply differential cross-section @@ -1509,6 +1529,19 @@ def plot_xray_thin_d3cross_ei_vs_Literature( [2] W. Nakel, Physics Reports, 243, p. 317—353, 1994 + Return: + ------- + dax: dict + Dict of axes + ddata_iso: dict + Dict of + ddata_ph_dist: dict + Dict of + ddata_ph_dist_nakel: dict + Dict of + ddata_ph_spect_nakel: dict + Dict of + """ # -------------- diff --git a/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated.py b/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated.py index 2cb60228c..ff9fcd339 100644 --- a/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated.py +++ b/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated.py @@ -1,6 +1,9 @@ import os +import copy +import warnings +from typing import Optional # Dict, Any import numpy as np @@ -14,6 +17,9 @@ from . import _xray_thin_target +TupleDict = tuple[dict] + + # #################################################### # #################################################### # DEFAULT @@ -32,6 +38,48 @@ _NTHETAE = 31 _NDPHI = 51 +# VERSION +_VERSION = 'BHE' # good compromise + + +# Default naming +_DSCALE = { + 12: 'T', + 9: 'G', + 6: 'M', + 3: 'k', + 0: '', + -3: 'm', + -6: 'u', + -9: 'n', +} + + +_DFORMAT = { + # energies + 'E_e0': { + 'data': None, + 'units': 'eV', + }, + 'E_ph': { + 'data': None, + 'units': 'eV', + }, + # angles + 'theta_ph': { + 'data': None, + 'units': 'rad', + }, + 'theta_e': { + 'data': None, + 'units': 'rad', + }, + 'dphi': { + 'data': None, + 'units': 'rad', + }, +} + # #################################################### # #################################################### @@ -40,25 +88,39 @@ def get_xray_thin_d2cross_ei_integrated_thetae_dphi( - # inputs - Z=None, + # optional input d2cross file + d2cross: Optional[str | dict] = None, + # target ion charge + Z: Optional[int] = None, + # energies E_e0_eV=None, E_ph_eV=None, theta_ph=None, # hypergeometric parameter - ninf=None, - source=None, + ninf: Optional[int] = None, + source: Optional[str] = None, # integration parameters nthetae=None, ndphi=None, # output customization - per_energy_unit=None, + per_energy_unit: Optional[str] = None, # version - version=None, + version: Optional[str] = None, # verb - verb=None, - verb_tab=None, -): + verb: Optional[bool] = None, + verb_tab: Optional[str] = None, + # saving + save: Optional[bool] = None, + pfe_save: Optional[str] = None, + overwrite: Optional[bool] = None, +) -> dict: + """ Compute d2cross, which is d3cross integrated over dphi + + Optionally loads / checks formatting of a pre-existing d2cross + - d2cross = str, should be a pfe to a local .npz + - d2cross = dict, will use as-is + + """ # ------------ # inputs @@ -67,120 +129,36 @@ def get_xray_thin_d2cross_ei_integrated_thetae_dphi( ( E_e0_eV, E_ph_eV, theta_ph, nthetae, ndphi, - shape, shape_theta_e, shape_dphi, + version, + shape, verb, verb_tab, - ) = _check( - # inputs - E_e0_eV=E_e0_eV, - E_ph_eV=E_ph_eV, - theta_ph=theta_ph, - # integration parameters - nthetae=nthetae, - ndphi=ndphi, - # verb - verb=verb, - verb_tab=verb_tab, - ) + save, pfe_save, overwrite, + ) = _check(**locals()) # ------------------ - # Derive angles + # compute # ------------------ - # E_e1_eV - E_e1_eV = E_e0_eV - E_ph_eV + if d2cross is None: - # angles - theta_e = np.pi * np.linspace(0, 1, nthetae) - dphi = np.pi * np.linspace(-1, 1, ndphi) - theta_ef = theta_e.reshape(shape_theta_e) - dphif = dphi.reshape(shape_dphi) + d2cross = _compute(**locals()) - # derived - sinte = np.sin(theta_ef) - - # ------------------ - # get d3cross - # ------------------ - - if verb >= 1: - msg = f"{verb_tab}Computing d3cross for shape {shape}... " - print(msg) - - d3cross = _xray_thin_target.get_xray_thin_d3cross_ei( - # inputs - Z=Z, - E_e0_eV=E_e0_eV[..., None, None], - E_e1_eV=E_e1_eV[..., None, None], - # directions - theta_ph=theta_ph[..., None, None], - theta_e=theta_ef, - dphi=dphif, - # hypergeometric parameter - ninf=ninf, - source=source, - # output customization - per_energy_unit=per_energy_unit, - # version - version=version, - # debug - debug=False, - ) - - # ------------------ - # prepare output - # ------------------ - - d2cross = { - # energies - 'E_e0': { - 'data': E_e0_eV, - 'units': 'eV', - }, - 'E_ph': { - 'data': E_ph_eV, - 'units': 'eV', - }, - # angles - 'theta_ph': { - 'data': theta_ph, - 'units': 'rad', - }, - 'theta_e': { - 'data': theta_e, - 'units': 'rad', - }, - 'dphi': { - 'data': dphi, - 'units': 'rad', - }, - # cross-section - 'cross': { - vv: { - 'data': np.full(shape, 0.), - 'units': asunits.Unit(vcross['units']) * asunits.Unit('sr'), - } - for vv, vcross in d3cross['cross'].items() - }, - } + # optional save + if save is True: + _save( + d2cross, + pfe_save=pfe_save, + overwrite=overwrite, + verb=verb, + ) # ------------------ - # integrate + # load # ------------------ - if verb >= 1: - msg = f"{verb_tab}Integrating..." - print(msg) + else: - for vv, vcross in d3cross['cross'].items(): - d2cross['cross'][vv]['data'][...] = scpinteg.simpson( - scpinteg.simpson( - vcross['data'] * sinte, - x=theta_e, - axis=-1, - ), - x=dphi, - axis=-1, - ) + d2cross = _load(d2cross) return d2cross @@ -199,9 +177,16 @@ def _check( # integration parameters nthetae=None, ndphi=None, + version=None, # verb verb=None, verb_tab=None, + # saving + save=None, + pfe_save=None, + overwrite=None, + # unused + **kwdargs, ): # ----------- @@ -240,14 +225,6 @@ def _check( theta_ph=theta_ph, ) - # ----------- - # shapes - # ----------- - - shape = np.broadcast_shapes(E_e0_eV.shape, E_ph_eV.shape, theta_ph.shape) - shape_theta_e = (1,) * (len(shape)+1) + (-1,) - shape_dphi = (1,) * len(shape) + (-1, 1) - # ----------- # integers # ----------- @@ -268,6 +245,22 @@ def _check( default=_NDPHI, ) + # ------------ + # version + # ------------ + + if version is None: + version = _VERSION + if isinstance(version, str): + version = [version] + + version = ds._generic_check._check_var_iter( + version, 'version', + types=(list, tuple), + types_iter=str, + allowed=['EH', 'BH', 'BHE'], + ) + # ----------- # verb # ----------- @@ -292,21 +285,425 @@ def _check( ) verb_tab = '\t'*verb_tab + # ----------- + # save + # ----------- + + savedef = pfe_save not in [None, False] + save = ds._generic_check._check_var( + save, 'save', + types=bool, + default=savedef, + ) + + # ----------- + # pfe_save + # ----------- + + if save is True: + if pfe_save is not None: + pfe_save = int(ds._generic_check._check_var( + pfe_save, 'pfe_save', + types=str, + )) + + try: + pfe_save = os.path.abspath(pfe_save) + except Exception as err: + msg = ( + "Arg 'pfe_save' must point to an valid path/file.ext\n" + f"Provided: {pfe_save}\n" + ) + raise Exception(msg) from err + + if not pfe_save.endswith('.npz'): + msg = ( + "Arg 'pfe_save' must point to an valid path/file.npz\n" + f"Provided: {pfe_save}\n" + ) + raise Exception(msg) + + else: + pfe_save = None + + # ----------- + # overwrite + # ----------- + + overwrite = ds._generic_check._check_var( + overwrite, 'overwrite', + types=bool, + default=False, + ) + return ( E_e0_eV, E_ph_eV, theta_ph, nthetae, ndphi, - shape, shape_theta_e, shape_dphi, + version, + shape, verb, verb_tab, + save, pfe_save, overwrite, ) +# #################################################### +# #################################################### +# compute +# #################################################### + + +def _compute( + Z=None, + E_e0_eV=None, + E_e1_eV=None, + E_ph_eV=None, + theta_ph=None, + # shapes + nthetae=None, + ndphi=None, + # parameters + ninf=None, + source=None, + version=None, + per_energy_unit=None, + # misc + shape=None, + daxis=None, + verb=None, + verb_tab=None, + # unused + **kwdargs, +): + + # ------------------ + # get angles and shape + # ------------------ + + # E_e1_eV + E_e1_eV = E_e0_eV - E_ph_eV + + # angles + theta_e = np.pi * np.linspace(0, 1, nthetae) + dphi = np.pi * np.linspace(-1, 1, ndphi) + + # --------------------- + # determine how to loop + # --------------------- + + # loop on largest dimension + iloop = np.argmax(shape) + sli = np.array([slice(None)]*len(shape) + [None, None]) + sli_Ee0 = np.copy(sli) + sli_Ee1 = np.copy(sli) + sli_theta = np.copy(sli) + slistr = [':'] * (len(shape) + 2) + + sli_ang = (None,) * (len(shape) - 1) + (slice(None),)*2 + theta_ef = theta_e[None, :][sli_ang] + dphif = dphi[:, None][sli_ang] + + # derived + sinte = np.sin(theta_ef) + + # ------------------ + # prepare output + # ------------------ + + d2cross = copy.deepcopy(_DFORMAT) + for kk, vv in _DFORMAT.items(): + kv = f"{kk}_{vv['units']}" if vv['units'] == 'eV' else kk + d2cross[kk]['data'] = eval(kv) + + # cross-sections + d2cross['cross'] = {vv: {'data': np.full(shape, 0.)} for vv in version} + srunits = asunits.Unit('sr') + + # ------------------ + # get d3cross + # ------------------ + + if verb >= 1: + msg = f"{verb_tab}Computing d3cross for shape {shape}... " + print(msg) + + # ------------------------------- + # loop on all but phi and theta_e + + size = shape[iloop] + for ii in range(size): + + # sli + sli[iloop] = ii + sli_Ee0[iloop] = min(ii, E_e0_eV.shape[iloop]-1) + sli_Ee1[iloop] = min(ii, E_e1_eV.shape[iloop]-1) + sli_theta[iloop] = min(ii, theta_ph.shape[iloop]-1) + + # verb + if verb >= 2: + slistr[iloop] = str(ii) + str0 = ', '.join(slistr) + str1 = str(shape + (None, None)) + end = '\n' if ii == size - 1 else '\r' + msg = f"\t{ii+1} / {size}, index ({str0}) / {str1}" + print(msg, end=end) + + # d3cross + d3cross = _xray_thin_target.get_xray_thin_d3cross_ei( + # inputs + Z=Z, + E_e0_eV=E_e0_eV[tuple(sli_Ee0)], + E_e1_eV=E_e1_eV[tuple(sli_Ee1)], + # directions + theta_ph=theta_ph[tuple(sli_theta)], + theta_e=theta_ef, + dphi=dphif, + # hypergeometric parameter + ninf=ninf, + source=source, + # output customization + per_energy_unit=per_energy_unit, + # version + version=version, + # debug + debug=False, + ) + + # integrate of theta_e + for vv, vcross in d3cross['cross'].items(): + d2cross['cross'][vv]['data'][tuple(sli[:-2])] = scpinteg.trapezoid( + scpinteg.trapezoid( + vcross['data'] * sinte, + x=theta_e, + axis=-1, + ), + x=dphi, + axis=-1, + ) + d2cross['cross'][vv]['units'] = ( + srunits * asunits.Unit(vcross['units']) + ) + + return d2cross + + +# #################################################### +# #################################################### +# save +# #################################################### + + +def _save( + d2cross=None, + pfe_save=None, + overwrite=None, + verb=None, +): + + # ---------- + # pfe_save + # ---------- + + if pfe_save is None: + + # extract sizes + ntheta = d2cross['theta_ph']['data'].size + Eph = _format_vect2str( + d2cross['E_ph']['data'], + base=d2cross['E_ph']['units'], + ) + Ee0 = _format_vect2str( + d2cross['E_e0']['data'], + base=d2cross['E_e0']['units'], + ) + + # extract boundaries + path = os.path.abspath(_PATH_HERE) + fname = f"d2cross_Ee0{Ee0}_Eph{Eph}_ntheta{ntheta}" + pfe_save = os.path.join(path, f"{fname}.npz") + + # ---------- + # overwrite + # ---------- + + if os.path.isfile(pfe_save): + if overwrite is True: + if verb is True: + msg = f"Overwritting file {pfe_save}\n" + warnings.warn(msg) + else: + msg = ( + "File {pfe_save} already exists!\n" + "\t=> use overwrite=True to overwrite\n" + ) + raise Exception(msg) + + # ---------- + # save + # ---------- + + np.savez(pfe_save, **d2cross) + + # ---------- + # verb + # ---------- + + if verb >= 1: + msg = f"Saved d2cross in:\n\t{pfe_save}\n" + print(msg) + + return + + +# #################################################### +# #################################################### +# Built str vector +# #################################################### + + +def _format_vect2str(vect, base='eV'): + + # ------------- + # extract + # ------------- + + v0 = vect.min() + v1 = vect.max() + nv = vect.size + + # ------------- + # test scale + # ------------- + + if np.max(vect.shape) == vect.size: + dlog = np.diff(np.log(vect)) + dlin = np.diff(vect) + islog = np.allclose(dlog, dlog[0]) + islin = np.allclose(dlin, dlin[0]) + + if islog is True: + scale = 'log' + elif islin is True: + scale = 'lin' + else: + scale = 'pts' + else: + scale = 'pts' + + # ------------- + # format + # ------------- + + v0 = _format(v0, base=base) + v1 = _format(v1, base=base) + + return f"{v0}-{v1}-{nv}{scale}" + + +def _format(vv, base='eV'): + + ls = sorted(_DSCALE.keys()) + ind = np.searchsorted(ls, np.log10(vv), side='right') - 1 + key = ls[ind] + factor = 10**(-key) + + return f"{vv*factor:3.0f}{_DSCALE[key]}{base}".strip() + + +# #################################################### +# #################################################### +# load +# #################################################### + + +def _load( + d2cross=None, +): + + # --------- + # str + # --------- + + if isinstance(d2cross, str): + + if not (os.path.isfile(d2cross) and d2cross.endswith('.npz')): + msg = ( + "Arg 'd2cross', if a str, should be valid path/file.npz\n" + f"Provided: {d2cross}\n" + ) + raise Exception(msg) + + d2cross = { + kk: vv.tolist() + for kk, vv in np.load(d2cross, allow_pickle=True).items() + } + + # --------- + # dict + # --------- + + if not isinstance(d2cross, dict): + msg = "Arg 'd2cross' must be a dict!\nProvided: {type(d2cross)}\n" + raise Exception(msg) + + # ---------------- + # inner formatting + # ---------------- + + dfail = {} + for kk in _DFORMAT.keys(): + if not isinstance(d2cross.get(kk), dict): + dfail[kk] = 'not a key or value not a dict ()' + elif str(d2cross[kk].get('units')) != _DFORMAT[kk]['units']: + dfail[kk] = 'no or wrong units ()' + elif not isinstance(d2cross[kk]['data'], np.ndarray): + dfail[kk] = "data is not a np.ndarray (type(d2cross[kk]['data']))" + else: + dfail[kk] = 'ok' + + if any([vv != 'ok' for vv in dfail.values()]) > 0: + lstr = [f"\t- {kk}: {vv}" for kk, vv in dfail.items()] + msg = ( + "Arg 'd2cross' must be a dict of subdicts " + "{'data': np.ndarray, 'units': str}, with keys:\n" + + "\n".join(lstr) + ) + raise Exception(msg) + + # ---------------- + # cross formatting + # ---------------- + + lok = ['BHE', 'BH', 'EH'] + typunits = (str, asunits.Unit, asunits.CompositeUnit) + c0 = ( + isinstance(d2cross.get('cross'), dict) + and all([ + kk in lok + and isinstance(vv, dict) + and isinstance(vv.get('data'), np.ndarray) + and isinstance(vv.get('units'), typunits) + for kk, vv in d2cross['cross'].items() + ]) + ) + if not c0: + msg = ( + "Arg d2cross['cross'] must be a dict " + "with {'version': dict} subdict with:\n" + f"\t- 'version' in {lok}\n" + f"Provided:\n{d2cross.get('cross')}\n" + ) + raise Exception(msg) + + return d2cross + + # #################################################### # #################################################### # plot vs litterature # #################################################### -def plot_xray_thin_d2cross_ei_vs_literature(): +def plot_xray_thin_d2cross_ei_vs_literature() -> TupleDict: """ Plot electron-angle-integrated cross section vs [1] G. Elwert and E. Haug, Phys. Rev., 183, pp. 90–105, 1969 diff --git a/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated_cases.py b/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated_cases.py new file mode 100644 index 000000000..59cb7c78b --- /dev/null +++ b/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated_cases.py @@ -0,0 +1,108 @@ + +# ANISOTROPY CASES +_DCASES_PRE = { + + # ------------- + # standard span + + 'standard': { + 0: { + 'E_e0_eV': 20e3, + 'E_ph_eV': 10e3, + 'color': 'r', + 'marker': '*', + 'ms': 14, + }, + 1: { + 'E_e0_eV': 100e3, + 'E_ph_eV': 50e3, + 'color': 'c', + 'marker': '*', + 'ms': 14, + }, + 2: { + 'E_e0_eV': 100e3, + 'E_ph_eV': 10e3, + 'color': 'm', + 'marker': '*', + 'ms': 14, + }, + 3: { + 'E_e0_eV': 1000e3, + 'E_ph_eV': 10e3, + 'color': (0.8, 0.8, 0), + 'marker': '*', + 'ms': 14, + }, + 4: { + 'E_e0_eV': 10000e3, + 'E_ph_eV': 10e3, + 'color': (0., 0.8, 0.8), + 'marker': '*', + 'ms': 14, + }, + 5: { + 'E_e0_eV': 1000e3, + 'E_ph_eV': 50e3, + 'color': (0.8, 0., 0.8), + 'marker': '*', + 'ms': 14, + }, + }, + + # ------------------ + # large span + + 'span_10eV_10MeV': { + 0: { + 'E_e0_eV': 100, + 'E_ph_eV': 50, + 'color': (1, 0, 0), + }, + 1: { + 'E_e0_eV': 10e3, + 'E_ph_eV': 50, + 'color': (0.8, 0, 0), + }, + 2: { + 'E_e0_eV': 10e3, + 'E_ph_eV': 5e3, + 'color': (0, 1, 0), + }, + 3: { + 'E_e0_eV': 600e3, + 'E_ph_eV': 50, + 'color': (0.6, 0, 0), + }, + 4: { + 'E_e0_eV': 600e3, + 'E_ph_eV': 5e3, + 'color': (0, 0.75, 0), + }, + 5: { + 'E_e0_eV': 600e3, + 'E_ph_eV': 450e3, + 'color': (0, 0, 1), + }, + 6: { + 'E_e0_eV': 7e6, + 'E_ph_eV': 50, + 'color': (0.4, 0, 0), + }, + 7: { + 'E_e0_eV': 7e6, + 'E_ph_eV': 5e3, + 'color': (0, 0.5, 0), + }, + 8: { + 'E_e0_eV': 7e6, + 'E_ph_eV': 450e3, + 'color': (0, 0, 0.5), + }, + 9: { + 'E_e0_eV': 7e6, + 'E_ph_eV': 5.5e6, + 'color': 'y', + }, + }, +} diff --git a/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated_dist_responsivity_plot.py b/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated_dist_responsivity_plot.py new file mode 100644 index 000000000..63eb89679 --- /dev/null +++ b/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated_dist_responsivity_plot.py @@ -0,0 +1,581 @@ + + +import copy +from typing import Optional # Any, Dict + + +import numpy as np +import astropy.units as asunits +import matplotlib.pyplot as plt +# import matplotlib.lines as mlines +import matplotlib.gridspec as gridspec +import datastock as ds + + +# from . import _xray_thin_target_integrated as _mod +from . import _xray_thin_target_integrated_plot as _mod_plot +# from ..distribution import get_distribution + + +TupleDict = tuple[dict] + + +# #################################################### +# #################################################### +# DEFAULT +# #################################################### + + +# Energy vectors +_E_E0_EV = np.logspace(3, 6, 51) +_E_PH_EV = np.linspace(1, 100, 25) * 1e3 +_THETA_PH = np.linspace(0, np.pi, 41) +_VERSION = 'BHE' + + +# DSCALES +_DSCALES = { + 'E_ph': 'log', + 'E_e0': 'log', + 'theta': 'linear', + 'dist': 'log', + 'resp': 'log', +} + + +# DCASES +_DCASES = { + 'cvd no filter maxwell': { + 'E_ph': {}, + 'responsivity': {}, + 'dist': {}, + 'E_e0': {}, + }, +} + + +# #################################################### +# #################################################### +# plot integrated filtered anisotropy +# #################################################### + + +def plot_xray_thin_integ_dist_filter_anisotropy( + # optional input d2cross file + d2cross: Optional[str | dict] = None, + # target ion charge + Z: Optional[int] = None, + # Energy + E_e0_eV=None, + E_ph_eV=None, + theta_ph=None, + # hypergeometric parameter + ninf: Optional[int] = None, + source: Optional[str] = None, + # output customization + per_energy_unit: Optional[str] = None, + # version + version: Optional[str] = None, + # selected cases + dcases_cross: Optional[dict[int, dict]] = None, + dcases_dist_resp: Optional[dict[int, dict]] = None, + # verb + verb: Optional[bool] = None, + # plot + dax: Optional[dict] = None, + fs: Optional[tuple] = None, + fontsize: Optional[int] = None, + E_e0_scale: Optional[str] = None, + E_ph_scale: Optional[str] = None, + dist_scale: Optional[str] = None, + resp_scale: Optional[str] = None, + theta_scale: Optional[str] = None, + dplot_forbidden=None, + dplot_peaking=None, + dplot_thetamax=None, + dplot_mean=None, +) -> TupleDict: + """ Compute and plot a (E_e0, E_ph) countour map of the d2cross section + + Where d2cross is the fully differentiated cross-section (d3cross), + integrated over one of the two the emission angle (dphi) + + Actually 3 overlayed contour plots with: + - integral of of the cross-section (over photon emission angle) + - angle of max cross-section + - peaking of the cross-section (std vs angle) + + Can overlay a few selected cases and plot them vs angle of emission + In normalized-linear and log scales + + """ + + # --------------- + # check inputs + # --------------- + + dcases_dist_resp, dscales, fs, fontsize = _check(**locals()) + + # --------------- + # prepare data + # --------------- + + # -------------- + # prepare axes + # -------------- + + if dax is None: + dax = _dax( + Z=Z, + version=version, + fs=fs, + fontsize=fontsize, + dscales=dscales, + ) + + dax = ds._generic_check._check_dax(dax) + + # ------------------- + # plot cross-section + # ------------------- + + dax_cross, d2cross_cross = _mod_plot.plot_xray_thin_d2cross_ei_anisotropy( + # optional input d2cross file + d2cross=d2cross, + # target ion charge + Z=Z, + # Energy + E_e0_eV=E_e0_eV, + E_ph_eV=E_ph_eV, + theta_ph=theta_ph, + # hypergeometric parameter + ninf=ninf, + source=source, + # output customization + per_energy_unit=per_energy_unit, + # version + version=version, + # selected cases + dcases=dcases_cross, + # verb + verb=verb, + # plot + dax=dax, + dplot_forbidden=dplot_forbidden, + dplot_peaking=dplot_peaking, + dplot_thetamax=dplot_thetamax, + dplot_mean=dplot_mean, + ) + + # ------------------- + # Compute integrand + # ------------------- + + dinteg = _integrand( + d2cross=d2cross, + dcases=dcases, + ) + + # ------------------- + # plot responsivity + # ------------------- + + kax = 'responsivity' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + for kk, vv in dcases_dist_resp.items(): + + l0, = ax.plot( + vv['responsivity']['data'], + vv['E_ph']['data']*1e-3, + c=vv.get('color'), + ls=vv.get('ls', '-'), + marker=vv.get('marker'), + lw=vv.get('lw', 1.), + label=kk, + ) + dcases_dist_resp[kk]['color'] = l0.get_color() + + ax.set_xlabel( + vv['responsivity']['units'], + fontsize=fontsize, + fontweight='bold', + ) + ax.invert_xaxis() + + # ------------------- + # plot distribution + # ------------------- + + kax = 'dist' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + for kk, vv in dcases_dist_resp.items(): + + ax.semilogy( + vv['E_e0']['data']*1e-3, + vv['dist']['data'], + c=vv['color'], + ls=vv.get('ls', '-'), + marker=vv.get('marker'), + lw=vv.get('lw', 1.), + label=kk, + ) + + ax.set_ylabel( + vv['dist']['units'], + fontsize=fontsize, + fontweight='bold', + ) + + return dax + + +# ############################################# +# ############################################# +# Integrand +# ############################################# + + +def _integrand( + d2cross=None, + dcases=None, +): + + # -------------- + # loop on cases + # -------------- + + for kcase, vcase in dcases.items(): + + pass + + return + + +# ############################################# +# ############################################# +# Axes for anisotropy +# ############################################# + + +def _check( + dcases_dist_resp=None, + fs=None, + fontsize=None, + E_e0_scale=None, + E_ph_scale=None, + dist_scale=None, + resp_scale=None, + theta_scale=None, + # unused + **kwdargs, +): + + # ------------ + # dcases + # ------------ + + ddef = copy.deepcopy(_DCASES) + if dcases_dist_resp in [None, False]: + dcases_dist_resp = {} + else: + for k0, v0 in dcases_dist_resp.items(): + dcases_dist_resp[k0] = _check_case( + v0, + f"dcases['{k0}']", + ddef[list(ddef.keys())[0]], + ) + + # ------------ + # fs + # ------------ + + if fs is None: + fs = (15, 12) + + # ------------ + # fontsize + # ------------ + + fontsize = ds._generic_check._check_var( + fontsize, 'fontsize', + types=int, + default=14, + sign='>0', + ) + + # ------------ + # scales + # ------------ + + dscales = { + 'E_ph': E_ph_scale, + 'E_e0': E_e0_scale, + 'theta': theta_scale, + 'dist': dist_scale, + 'resp': resp_scale, + } + + for kk, vv in dscales.items(): + dscales[kk] = ds._generic_check._check_var( + vv, f'{kk}_scale', + types=str, + allowed=['log', 'linear'], + default=_DSCALES[kk], + ) + + return dcases_dist_resp, dscales, fs, fontsize + + +def _check_case( + case=None, + key=None, + ddef=None, +): + + # -------------- + # general structure + # -------------- + + dfail = {} + lok = list(ddef.keys()) + ltunits = (str, asunits.Unit, asunits.CompositeUnit) + for kk in lok: + if not isinstance(case.get(kk), dict): + typ = type(case.get(kk)) + dfail[kk] = f'absent or not a dict ({typ})' + elif not isinstance(case[kk].get('data'), np.ndarray): + typ = type(case[kk].get('data')) + dfail[kk] = f'data key not a np.ndarray ({typ})' + elif not isinstance(case[kk].get('units'), ltunits): + typ = type(case[kk].get('units')) + dfail[kk] = f"units not a str ({typ})" + else: + dfail[kk] = 'ok' + + if any([vv != 'ok' for vv in dfail.values()]): + lstr = [f"\t- {kk}: {vv}" for kk, vv in dfail.items()] + msg = ( + f"Arg {key} must be a dict with keys {lok}, " + "where each is a {'data': np.ndarray, 'units': str} subdict!\n" + + "\n".join(lstr) + ) + raise Exception(msg) + + # -------------- + # shape consistency + # -------------- + + shape_Eph = case['E_ph']['data'].shape + shape_resp = case['responsivity']['data'].shape + if shape_Eph != shape_resp: + msg = ( + "The 2 fields below must have the same shape:\n" + f"{key}['E_ph']['data'].shape = {shape_Eph}\n" + f"{key}['responsivity']['data'].shape = {shape_resp}\n" + ) + raise Exception(msg) + + shape_Ee0 = case['E_e0']['data'].shape + shape_dist = case['dist']['data'].shape + if shape_Ee0 != shape_dist: + msg = ( + "The 2 fields below must have the same shape:\n" + f"{key}['E_e0']['data'].shape = {shape_Ee0}\n" + f"{key}['dist']['data'].shape = {shape_dist}\n" + ) + raise Exception(msg) + + return case + + +# ############################################# +# ############################################# +# Axes for anisotropy +# ############################################# + + +def _dax( + Z=None, + version=None, + fs=None, + fontsize=None, + dscales=None, +): + + # -------------- + # prepare figure + # -------------- + + tit = ( + "Thin-target Bremsstrahlung emission anisotropy" + ) + + dmargin = { + 'left': 0.05, 'right': 0.98, + 'bottom': 0.05, 'top': 0.92, + 'wspace': 0.30, 'hspace': 0.20, + } + + fig = plt.figure(figsize=(15, 12)) + fig.suptitle(tit, size=fontsize+2, fontweight='bold') + + nh0, nh1, nh2 = 2, 6, 4 + nv0, nv1, nv2 = 3, 1, 2 + gs = gridspec.GridSpec( + ncols=nh0+nh1+2*(nh2 + 1), + nrows=nv0 + nv1, + **dmargin, + ) + dax = {} + + # -------------- + # prepare axes + # -------------- + + # -------------- + # ax - isolines + + ax = fig.add_subplot( + gs[:nv0, nh0:nh0+nh1], + xscale=dscales['E_e0'], + yscale=dscales['E_ph'], + ) + ax.set_title( + r"$d^2\sigma(E_{e0}, E_{ph}, \theta_{ph}, Z)$" + + f"\n Z = {Z}, version = {version}", + size=fontsize, + fontweight='bold', + ) + + # store + dax['map'] = {'handle': ax, 'type': 'isolines'} + + # -------------- + # ax - responsivity + + ax = fig.add_subplot( + gs[:nv0, :nh0], + sharey=dax['map']['handle'], + xscale=dscales['resp'], + ) + ax.set_title( + "responsivity", + size=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + r"$E_{ph}$ (keV)", + size=fontsize, + fontweight='bold', + ) + ax.grid(True) + + # store + dax['responsivity'] = {'handle': ax, 'type': 'isolines'} + + # -------------- + # ax - dist + + ax = fig.add_subplot( + gs[nv0:, nh0:nh0 + nh1], + sharex=dax['map']['handle'], + yscale=dscales['dist'], + ) + ax.set_xlabel( + r"$E_{e,0}$ (keV)", + size=fontsize, + fontweight='bold', + ) + ax.grid(True) + + # store + dax['dist'] = {'handle': ax, 'type': 'isolines'} + + # -------------- + # ax - theta_cross - norm + + ax = fig.add_subplot( + gs[:nv2, nh0 + nh1 + 1:nh0 + nh1 + 1 + nh2], + xscale=dscales['theta'], + yscale='linear', + ) + ax.set_xlabel( + r"$\theta_{ph}^{e0}$ (deg)", + size=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + "cross-section norm.", + size=fontsize, + fontweight='bold', + ) + + # store + dax['theta_norm'] = {'handle': ax, 'type': 'isolines'} + + # -------------- + # ax - theta_cross - abs + + ax = fig.add_subplot( + gs[nv2:, nh0 + nh1 + 1:nh0 + nh1 + 1 + nh2], + sharex=dax['theta_norm']['handle'], + yscale='log', + ) + ax.set_xlabel( + r"$\theta_{ph}^{e0}$ (deg)", + size=fontsize, + fontweight='bold', + ) + + # store + dax['theta_abs'] = {'handle': ax, 'type': 'isolines'} + + # -------------- + # ax - theta_emiss - norm + + ax = fig.add_subplot( + gs[:nv2, nh0 + nh1 + 2 + nh2:], + xscale=dscales['theta'], + yscale='linear', + ) + ax.set_xlabel( + r"$\theta_{ph}^B$ (deg)", + size=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + "emiss norm.", + size=fontsize, + fontweight='bold', + ) + + # store + dax['theta_emiss_norm'] = {'handle': ax, 'type': 'isolines'} + + # -------------- + # ax - theta_emiss - abs + + ax = fig.add_subplot( + gs[nv2:, nh0 + nh1 + 2 + nh2:], + sharex=dax['theta_norm']['handle'], + yscale='log', + ) + ax.set_xlabel( + r"$\theta_{ph}^B$ (deg)", + size=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + "emiss (ph/sr/s/m3)", + size=fontsize, + fontweight='bold', + ) + + # store + dax['theta_emiss_abs'] = {'handle': ax, 'type': 'isolines'} + + return dax diff --git a/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated_plot.py b/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated_plot.py index c967ae16b..8688d968f 100644 --- a/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated_plot.py +++ b/tofu/physics_tools/electrons/emission/_xray_thin_target_integrated_plot.py @@ -1,6 +1,7 @@ import copy +from typing import Optional # Any, Dict import numpy as np @@ -12,7 +13,11 @@ import datastock as ds -from . import _xray_thin_target_integrated +from . import _xray_thin_target_integrated as _mod +from ._xray_thin_target_integrated_cases import _DCASES_PRE + + +TupleDict = tuple[dict] # #################################################### @@ -21,51 +26,36 @@ # #################################################### -# ANISOTROPY CASES -_DCASES = { - 0: { - 'E_e0_eV': 20e3, - 'E_ph_eV': 10e3, - 'color': 'r', - 'marker': '*', - 'ms': 14, - }, - 1: { - 'E_e0_eV': 100e3, - 'E_ph_eV': 50e3, - 'color': 'c', - 'marker': '*', - 'ms': 14, - }, - 2: { - 'E_e0_eV': 100e3, - 'E_ph_eV': 10e3, - 'color': 'm', - 'marker': '*', - 'ms': 14, - }, - 3: { - 'E_e0_eV': 1000e3, - 'E_ph_eV': 10e3, - 'color': (0.8, 0.8, 0), - 'marker': '*', - 'ms': 14, - }, - 4: { - 'E_e0_eV': 10000e3, - 'E_ph_eV': 10e3, - 'color': (0., 0.8, 0.8), - 'marker': '*', - 'ms': 14, - }, - 5: { - 'E_e0_eV': 1000e3, - 'E_ph_eV': 50e3, - 'color': (0.8, 0., 0.8), - 'marker': '*', - 'ms': 14, - }, +# Energy vectors +_E_E0_EV = np.logspace(3, 6, 51) +_E_PH_EV = np.linspace(1, 100, 25) * 1e3 +_THETA_PH = np.linspace(0, np.pi, 41) +_VERSION = 'BHE' + + +# DSCALES +_DSCALES = { + 'E_ph': 'log', + 'E_e0': 'log', + 'theta': 'linear', } + + +# ANISOTROPY CASES FORMATTING +_DCASES_FORMAT = { + 'E_e0_eV': {'type': (int, float), 'val': 1e3}, + 'E_ph_eV': {'type': (int, float), 'val': 10e3}, + 'color': {'type': (str, tuple), 'val': 'k'}, + 'marker': {'type': str, 'val': '*'}, + 'ms': {'type': (int, float), 'val': 18}, + 'ls': {'type': str, 'val': '-'}, +} + + +# ANISOTROPY CASES DEFAULT +_DCASES_CASE = 'standard' + + # #################################################### # #################################################### # plot anisotropy @@ -73,26 +63,50 @@ def plot_xray_thin_d2cross_ei_anisotropy( - # compute - Z=None, + # optional input d2cross file + d2cross: Optional[str | dict] = None, + # target ion charge + Z: Optional[int] = None, + # Energy E_e0_eV=None, E_ph_eV=None, theta_ph=None, - per_energy_units=None, - version=None, - # hypergeometrc - ninf=None, - source=None, + # hypergeometric parameter + ninf: Optional[int] = None, + source: Optional[str] = None, + # output customization + per_energy_unit: Optional[str] = None, + # version + version: Optional[str] = None, # selected cases - dcases=None, + dcases: Optional[dict[int, dict]] = None, + # verb + verb: Optional[bool] = None, # plot - dax=None, - fontsize=None, + dax: Optional[dict] = None, + fontsize: Optional[int] = None, + E_e0_scale: Optional[str] = None, + E_ph_scale: Optional[str] = None, + theta_scale: Optional[str] = None, dplot_forbidden=None, dplot_peaking=None, dplot_thetamax=None, - dplot_integ=None, -): + dplot_mean=None, +) -> TupleDict: + """ Compute and plot a (E_e0, E_ph) countour map of the d2cross section + + Where d2cross is the fully differentiated cross-section (d3cross), + integrated over one of the two the emission angle (dphi) + + Actually 3 overlayed contour plots with: + - integral of of the cross-section (over photon emission angle) + - angle of max cross-section + - peaking of the cross-section (std vs angle) + + Can overlay a few selected cases and plot them vs angle of emission + In normalized-linear and log scales + + """ # --------------- # check inputs @@ -100,44 +114,49 @@ def plot_xray_thin_d2cross_ei_anisotropy( ( E_e0_eV, E_ph_eV, theta_ph, - dcases, + version, + dscales, + verb, fontsize, - dplot_forbidden, dplot_peaking, dplot_thetamax, dplot_integ, - ) = _check_anisotropy( - E_e0_eV=E_e0_eV, - E_ph_eV=E_ph_eV, - theta_ph=theta_ph, - version=version, - # selected cases - dcases=dcases, - # plotting - fontsize=fontsize, - dplot_forbidden=dplot_forbidden, - dplot_peaking=dplot_peaking, - dplot_thetamax=dplot_thetamax, - dplot_integ=dplot_integ, - ) + dplot_forbidden, dplot_peaking, dplot_thetamax, dplot_mean, + ) = _check_anisotropy(**locals()) # --------------- # prepare data # --------------- - mod = _xray_thin_target_integrated - d2cross = mod.get_xray_thin_d2cross_ei_integrated_thetae_dphi( + d2cross = _mod.get_xray_thin_d2cross_ei_integrated_thetae_dphi( + d2cross=d2cross, # inputs Z=Z, E_e0_eV=E_e0_eV[None, :, None], E_ph_eV=E_ph_eV[None, None, :], theta_ph=theta_ph[:, None, None], # output customization - per_energy_unit=per_energy_units, + per_energy_unit=per_energy_unit, # version version=version, # hypergeometric ninf=ninf, source=source, # verb - verb=False, + verb=verb, + ) + + # ------------------- + # update from d2cross + # ------------------- + + # if d2cross was provided + theta_ph = d2cross['theta_ph']['data'].ravel() + E_ph_eV = d2cross['E_ph']['data'].ravel() + E_e0_eV = d2cross['E_e0']['data'].ravel() + + # dcases + dcases = _check_dcases( + dcases=dcases, + E_e0_eV=E_e0_eV, + E_ph_eV=E_ph_eV, ) # -------------- @@ -145,12 +164,15 @@ def plot_xray_thin_d2cross_ei_anisotropy( # -------------- if dax is None: - dax = _get_axes_anisotropy( + dax = _dax( Z=Z, version=version, fontsize=fontsize, + dscales=dscales, ) + dax = ds._generic_check._check_dax(dax) + # --------------- # plot - map # --------------- @@ -162,20 +184,25 @@ def plot_xray_thin_d2cross_ei_anisotropy( for iv, (kk, vv) in enumerate(d2cross['cross'].items()): # compute integral and peaking - integ, peaking = _get_peaking( + mean, peaking = _get_peaking( vv['data'], theta_ph*180/np.pi, axis=0, ) + mean_log10 = np.full(mean.shape, np.nan) + iok = np.isfinite(mean) + iok[iok] = mean[iok] > 0. + mean_log10[iok] = np.log10(mean[iok]) + mean_units = vv['units'] # integral - if dplot_integ is not False: + if dplot_mean is not False: im0 = ax.contour( E_e0_eV * 1e-3, E_ph_eV * 1e-3, - np.log10(integ).T, - levels=dplot_integ['levels'], - colors=dplot_integ['colors'], + mean_log10.T, + levels=dplot_mean['levels'], + colors=dplot_mean['colors'], ) # clabels @@ -240,24 +267,27 @@ def plot_xray_thin_d2cross_ei_anisotropy( ax.add_patch(patch) # legend - lh = [ - mlines.Line2D( + lh = [] + if dplot_mean is not False: + lh.append(mlines.Line2D( [], [], - c=dplot_integ['colors'], - label='log10(integral)', - ), - mlines.Line2D( + c=dplot_mean['colors'], + label=f'log10() (log10({mean_units}))', + )) + if dplot_peaking is not False: + lh.append(mlines.Line2D( [], [], c=dplot_peaking['colors'], label='peaking (1/std)', - ), - mlines.Line2D( + )) + if dplot_thetamax is not False: + lh.append(mlines.Line2D( [], [], c=dplot_thetamax['colors'], label='theta_max (deg)', - ), - ] - ax.legend(handles=lh, loc='upper left') + )) + if len(lh) > 0: + ax.legend(handles=lh, loc='upper left') # add cases for ic, (kcase, vcase) in enumerate(dcases.items()): @@ -273,7 +303,8 @@ def plot_xray_thin_d2cross_ei_anisotropy( ) # limits - ax.set_ylim(0, ymax*1e-3) + if dscales['E_ph'] == 'linear': + ax.set_ylim(0, ymax*1e-3) # --------------- # plot - cases @@ -287,8 +318,8 @@ def plot_xray_thin_d2cross_ei_anisotropy( yy = vv['data'][:, vcase['ie'], vcase['iph']] if np.any(yy > 0): - # normalized - kax = 'norm' + # theta_norm + kax = 'theta_norm' if dax.get(kax) is not None: ax = dax[kax]['handle'] @@ -296,11 +327,12 @@ def plot_xray_thin_d2cross_ei_anisotropy( theta_ph * 180/np.pi, yy / np.max(yy), c=vcase['color'], + ls=vcase['ls'], label=labi, ) - # abs - kax = 'log' + # theta_abs + kax = 'theta_abs' if dax.get(kax) is not None: ax = dax[kax]['handle'] @@ -308,11 +340,12 @@ def plot_xray_thin_d2cross_ei_anisotropy( theta_ph * 180/np.pi, yy*1e28, c=vcase['color'], + ls=vcase['ls'], label=labi, ) # normalized - kax = 'norm' + kax = 'theta_norm' if dax.get(kax) is not None: ax = dax[kax]['handle'] ax.legend(prop={'size': 12}) @@ -320,7 +353,7 @@ def plot_xray_thin_d2cross_ei_anisotropy( ax.set_xlim(0, 180) # normalized - kax = 'abs' + kax = 'theta_abs' if dax.get(kax) is not None: ax = dax[kax]['handle'] ax.legend(prop={'size': 12}) @@ -347,73 +380,61 @@ def _check_anisotropy( E_ph_eV=None, theta_ph=None, version=None, - # selected cases - dcases=None, + # verb + verb=None, + # scales + E_e0_scale=None, + E_ph_scale=None, + theta_scale=None, # plotting fontsize=None, dplot_forbidden=None, dplot_peaking=None, dplot_thetamax=None, - dplot_integ=None, + dplot_mean=None, + # unused + **kwdargs, ): # E_e0_eV if E_e0_eV is None: - E_e0_eV = np.logspace(3, 6, 51) + E_e0_eV = _E_E0_EV E_e0_eV = ds._generic_check._check_flat1darray( E_e0_eV, 'E_e0_eV', dtype=float, sign='>0', + unique=True, ) # E_ph_eV if E_ph_eV is None: - E_ph_eV = np.linspace(1, 100, 25) * 1e3 + E_ph_eV = _E_PH_EV E_ph_eV = ds._generic_check._check_flat1darray( E_ph_eV, 'E_ph_eV', dtype=float, sign='>0', + unique=True, ) # theta_ph if theta_ph is None: - theta_ph = np.linspace(0, np.pi, 41) + theta_ph = _THETA_PH # version if version is None: - version = 'BHE' - - # ------------ - # dcases - # ------------ + version = _VERSION - ddef = copy.deepcopy(_DCASES) - if dcases in [None, True]: - dcases = ddef + # ----------- + # verb + # ----------- - if dcases is not False: - for k0, v0 in dcases.items(): - dcases[k0] = _check_anisotropy_dplot( - v0, - f'dcases[{k0}]', - ddef[0], - ) - - # update with indices - ie = np.argmin(np.abs(E_e0_eV - dcases[k0]['E_e0_eV'])) - iph = np.argmin(np.abs(E_ph_eV - dcases[k0]['E_ph_eV'])) - dcases[k0].update({'ie': ie, 'iph': iph}) - - # update with label - ee0 = E_e0_eV[ie] - eph = E_ph_eV[iph] - dcases[k0]['lab'] = ( - r"$E_{e0} / E_{ph}$ = " - + f"{ee0*1e-3:3.0f} / {eph*1e-3:3.0f} keV = " - + f"{round(ee0 / eph, ndigits=1)}" - ) + verb = ds._generic_check._check_var( + verb, 'verb', + types=(bool, int), + default=False, + ) # ------------ # plotting @@ -452,22 +473,105 @@ def _check_anisotropy( ddef, ) - # dplot_integ + # dplot_mean ddef = {'colors': 'g', 'levels': 20} - dplot_integ = _check_anisotropy_dplot( - dplot_integ, - 'dplot_integ', + dplot_mean = _check_anisotropy_dplot( + dplot_mean, + 'dplot_mean', ddef, ) + # ------------ + # scales + # ------------ + + dscales = { + 'E_ph': E_ph_scale, + 'E_e0': E_e0_scale, + 'theta': theta_scale, + } + + for kk, vv in dscales.items(): + dscales[kk] = ds._generic_check._check_var( + vv, f'{kk}_scale', + types=str, + allowed=['log', 'linear'], + default=_DSCALES[kk], + ) + return ( E_e0_eV, E_ph_eV, theta_ph, - dcases, + version, + dscales, + verb, fontsize, - dplot_forbidden, dplot_peaking, dplot_thetamax, dplot_integ, + dplot_forbidden, dplot_peaking, dplot_thetamax, dplot_mean, ) +def _check_dcases( + dcases=None, + E_e0_eV=None, + E_ph_eV=None, +): + + # -------------- + # dcases default + # -------------- + + ddef = copy.deepcopy(_DCASES_FORMAT) + if dcases in [None, True]: + dcases = _DCASES_CASE + + # -------------- + # dcases from predefined + # -------------- + + if isinstance(dcases, str): + lok = sorted(_DCASES_PRE.keys()) + if dcases not in lok: + lstr = [f"\t- {kk}" for kk in lok] + msg = ( + "Arg 'dcases' must be either:\n" + "\t- dict of cases\n" + "\t- a key to a predefined dict of cases\n" + "Available predefined dcases:\n" + + "\n".join(lstr) + ) + raise Exception(msg) + dcases = copy.deepcopy(_DCASES_PRE[dcases]) + + # -------------- + # generic check + # -------------- + + if dcases is not False: + for k0, v0 in dcases.items(): + dcases[k0] = _check_anisotropy_dplot( + v0, + f'dcases[{k0}]', + ddef, + ) + + # update with indices + ie = np.argmin(np.abs(E_e0_eV - dcases[k0]['E_e0_eV'])) + iph = np.argmin(np.abs(E_ph_eV - dcases[k0]['E_ph_eV'])) + dcases[k0].update({'ie': ie, 'iph': iph}) + + # update with label + ee0 = E_e0_eV[ie] + eph = E_ph_eV[iph] + dcases[k0]['lab'] = ( + r"$E_{e0} / E_{ph}$ = " + + f"{ee0*1e-3:3.0f} / {eph*1e-3:3.0f} keV = " + + f"{round(ee0 / eph, ndigits=1)}" + ) + else: + dcases = {} + + return dcases + + def _check_anisotropy_dplot(din, dname, ddef): # ------------- @@ -484,10 +588,18 @@ def _check_anisotropy_dplot(din, dname, ddef): if din is not False: c0 = ( isinstance(din, dict) - and all([kk in ddef.keys() for kk in din.keys()]) + and all([ + kk in ddef.keys() + and ( + isinstance(ddef[kk], dict) + and ddef[kk].get('type') is not None + and isinstance(din[kk], ddef[kk]['type']) + ) + for kk in din.keys() + ]) ) if not c0: - lstr = [f"\t- ''{k0}': {v0}" for k0, v0 in ddef.items()] + lstr = [f"\t- '{k0}': {v0['type']}" for k0, v0 in ddef.items()] msg = ( f"Arg '{dname}' must be either False or a dict with:\n" + "\n".join(lstr) @@ -501,7 +613,11 @@ def _check_anisotropy_dplot(din, dname, ddef): if din is not False: for k0, v0 in ddef.items(): - din[k0] = din.get(k0, v0) + if isinstance(v0, dict): + vv = v0['val'] + else: + vv = v0 + din[k0] = din.get(k0, vv) return din @@ -519,11 +635,18 @@ def _get_peaking(data, x, axis=None): # ---------- integ = scpinteg.trapezoid(data, x=x, axis=axis) - shape_integ = tuple([ - 1 if ii == axis else ss - for ii, ss in enumerate(data.shape) - ]) - data_n = data / integ.reshape(shape_integ) + shape_integ = list(data.shape) + shape_integ[axis] = 1 + + data_n = np.full(data.shape, np.nan) + iok = np.isfinite(integ) + iok[iok] = integ[iok] > 0 + iokn = iok.nonzero() + sli0 = list(iokn) + sli1 = list(iokn) + sli0.insert(axis, None) + sli1.insert(axis, slice(None)) + data_n[tuple(sli1)] = data[tuple(sli1)] / integ[tuple(sli0)] # ---------- # get average @@ -531,10 +654,14 @@ def _get_peaking(data, x, axis=None): shape_x = tuple([-1 if ii == axis else 1 for ii in range(data.ndim)]) xf = x.reshape(shape_x) - x_avf = scpinteg.simpson(data_n * xf, x=x, axis=axis).reshape(shape_integ) + x_avf = scpinteg.trapezoid( + data_n * xf, + x=x, + axis=axis, + ).reshape(shape_integ) std = np.sqrt(scpinteg.simpson(data_n * (xf - x_avf)**2, x=x, axis=axis)) - return integ, 1/std + return integ/180, 1/std # ############################################# @@ -543,20 +670,22 @@ def _get_peaking(data, x, axis=None): # ############################################# -def _get_axes_anisotropy( +def _dax( Z=None, version=None, fontsize=None, + dax=None, + dscales=None, ): tit = ( - "Anisotropy" + "Thin-target Bremsstrahlung cross-section anisotropy" ) dmargin = { - 'left': 0.08, 'right': 0.95, - 'bottom': 0.06, 'top': 0.85, - 'wspace': 0.2, 'hspace': 0.40, + 'left': 0.06, 'right': 0.95, + 'bottom': 0.06, 'top': 0.90, + 'wspace': 0.20, 'hspace': 0.20, } fig = plt.figure(figsize=(15, 12)) @@ -572,7 +701,11 @@ def _get_axes_anisotropy( # -------------- # ax - isolines - ax = fig.add_subplot(gs[:, 0], xscale='log') + ax = fig.add_subplot( + gs[:, 0], + xscale=dscales['E_e0'], + yscale=dscales['E_ph'], + ) ax.set_xlabel( r"$E_{e,0}$ (keV)", size=fontsize, @@ -594,9 +727,12 @@ def _get_axes_anisotropy( dax['map'] = {'handle': ax, 'type': 'isolines'} # -------------- - # ax - norm + # ax - theta_norm - ax = fig.add_subplot(gs[0, 1]) + ax = fig.add_subplot( + gs[0, 1], + xscale=dscales['theta'], + ) ax.set_xlabel( r"$\theta_{ph}$ (deg)", size=fontsize, @@ -614,22 +750,20 @@ def _get_axes_anisotropy( ) # store - dax['norm'] = {'handle': ax, 'type': 'isolines'} + dax['theta_norm'] = {'handle': ax, 'type': 'isolines'} # -------------- - # ax - log + # ax - theta_abs - ax = fig.add_subplot(gs[1, 1], sharex=dax['norm']['handle']) + ax = fig.add_subplot( + gs[1, 1], + sharex=dax['theta_norm']['handle'], + ) ax.set_xlabel( r"$\theta_{ph}$ (deg)", size=fontsize, fontweight='bold', ) - ax.set_ylabel( - r"$\frac{d^2\sigma_{ei}}{dkd\Omega_{ph}}$ ()", - size=fontsize, - fontweight='bold', - ) ax.set_title( "Absolute cross-section vs photon emission angle", size=fontsize, @@ -637,6 +771,6 @@ def _get_axes_anisotropy( ) # store - dax['log'] = {'handle': ax, 'type': 'isolines'} + dax['theta_abs'] = {'handle': ax, 'type': 'isolines'} return dax