Issue1166 physicselectrons

tofu/physics_tools/electrons/distribution/_distribution.py

@@ -161,13 +161,25 @@ def main(
         inan = np.isnan(ddist['dist'][kdist]['dist']['data'])
         ddist['dist'][kdist]['dist']['data'][inan] = 0.
 
+        # neg => error
+        ineg = ddist['dist'][kdist]['dist']['data'] < 0.
+        if np.any(ineg):
+            msg = (
+                "Electron dist has negative values!\n"
+                f"\t- dist = '{kdist}'\n"
+                f"\t- version = '{dfunc[kdist]['version']}'\n"
+                f"\t- module = {dfunc[kdist]['func'].__module__}\n"
+                f"\t- func = {dfunc[kdist]['func'].__name__}\n"
+            )
+            raise Exception(msg)
+
         # scale
         ne_re = _scale(
             din=din,
             ddist=ddist,
             kdist=kdist,
             dcoords=dcoords,
-            version=version,
+            version=dfunc[kdist]['version'],
         )
 
     # --------------
@@ -185,7 +197,7 @@ def main(
             ddist=ddist,
             kdist=kdist,
             dcoords=dcoords,
-            version=version,
+            version=dfunc[kdist]['version'],
         )
 
         # store
@@ -355,7 +367,10 @@ def _get_velocity_par(ddist, kdist):
             energy_kinetic_eV=ddist['coords']['x0']['data'],
         )['velocity_ms']
         units = v_par_ms['units']
-        v_par_ms = v_par_ms['data']
+        # v_par_ms = v_par_ms['data']
+
+        # assume 0 drift velocity => average v_par = 0
+        v_par_ms = np.zeros(v_par_ms['data'].shape)
 
     else:
         raise NotImplementedError(kcoords)
@@ -401,6 +416,7 @@ def _integrate(
         )
         ne = ddist['dist'][kdist]['dist']['data']
         x0 = dcoords['x0']['data']
+
     else:
         current = scpinteg.trapezoid(
             scpct.e

tofu/physics_tools/electrons/distribution/_distribution_check.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import astropy.units as asunits
+import scipy.constants as scpct
 import datastock as ds
 import tofu as tf
 
@@ -397,6 +398,40 @@ def _coords(
         )
         raise Exception(msg)
 
+    # --------------
+    # check values
+    # --------------
+
+    for k0, v0 in dcoords.items():
+
+        # check validity
+        if k0 == 'theta':
+            iout = (v0['data'] < 0.) | (v0['data'] > np.pi)
+            msg = "must be in [0; pi]"
+        elif k0 == 'pitch':
+            iout = (v0['data'] < -1) | (v0['data'] > 1)
+            msg = "must be in [-1; 1]"
+        elif k0.startswith('v_par'):
+            iout = (v0['data'] > scpct.c)
+            msg = "must be <= c"
+        elif k0.startswith('v_perp'):
+            iout = (v0['data'] > scpct.c) | (v0['data'] < 0.)
+            msg = "must be in [0; c]"
+        elif k0.startswith('p_perp') or k0 == 'E_eV':
+            iout = (v0['data'] < 0.)
+            msg = "must be >= 0"
+        else:
+            msg = ''
+            continue
+
+        # Exception
+        if np.any(iout):
+            msg = (
+                f"Some values in dcoords['{k0}'] are invalid:\n"
+                + msg + "\n"
+            )
+            raise Exception(msg)
+
     return dcoords
 
 

tofu/physics_tools/electrons/distribution/_distribution_maxwell.py

@@ -209,6 +209,7 @@ def f3d_E_theta(
         v0_par_ms=v0_par_ms,
     )
 
+    # caution: sin(theta) => negative values !
     dist = np.sin(theta) * dist0 / (2.*np.pi)
     units = units0 * asunits.Unit('1/rad^2')
 

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

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,35 +25,40 @@
 _PATH_HERE = os.path.dirname(__file__)
 
 
+# version
+_VERSION = 'EH'  # most accurate, but slow
+
+
 # ####################################################
 # ####################################################
 #        Differential cross-section
 # ####################################################
 
 
 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
     theta_ph=None,
     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,23 +901,29 @@ 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
     # -------------
 
     # 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
+
     """
 
     # --------------