Issue1166 physicselectrons

tofu/data/_class08_move3d_check.py

@@ -227,16 +227,15 @@ def _add_asis(
         if isinstance(kk, tuple):
             for ii, ki in enumerate(kk):
                 k0 = ki.split('_')[0]
-                if dgeom0.get(k0) is not None:
+                if isinstance(dgeom0.get(k0), tuple):
                     dgeom[ki] = coll.ddata[dgeom0[k0][ii]]['data']
 
         elif dgeom0.get(kk) is not None:
 
             if isinstance(dgeom0[kk], str):
                 dgeom[kk] = coll.ddata[kk]['data']
 
-            else:
-                if dgeom0.get(kk) is not None:
-                    dgeom[kk] = dgeom0[kk]
+            elif not isinstance(dgeom0[kk], tuple):
+                dgeom[kk] = dgeom0[kk]
 
     return

tofu/data/_class08_move_rotate3d_by.py

@@ -268,26 +268,27 @@ def _rotate_camera(
     # unit vects
     lk_vect = ['nin', 'e0', 'e1']
     for kk in lk_vect:
-        if dgeom0.get(kk) is not None:
-            dgeom[kk] = np.r_[_rotate_pts(
+
+        if isinstance(dgeom0.get(kk), tuple):
+            kx, ky, kz = f"{kk}_x", f"{kk}_y", f"{kk}_z"
+            dgeom[kx], dgeom[ky], dgeom[kz] = _rotate_pts(
                 axis_pt,
                 axis_vect,
                 angle,
-                *dgeom0[kk],
+                coll.ddata[dgeom0[kk][0]]['data'],
+                coll.ddata[dgeom0[kk][1]]['data'],
+                coll.ddata[dgeom0[kk][2]]['data'],
                 isvect=True,
-            )]
+            )
 
-        if dgeom0.get(f"{kk}_x") is not None:
-            kx, ky, kz = f"{kk}_x", f"{kk}_y", f"{kk}_z"
-            dgeom[kx], dgeom[ky], dgeom[kz] = _rotate_pts(
+        elif dgeom0.get(kk) is not None:
+            dgeom[kk] = np.r_[_rotate_pts(
                 axis_pt,
                 axis_vect,
                 angle,
-                dgeom0[kx],
-                dgeom0[ky],
-                dgeom0[kz],
+                *dgeom0[kk],
                 isvect=True,
-            )
+            )]
 
     # ----------------
     # add as-is

tofu/data/_class08_move_translate3d_by.py

@@ -212,9 +212,9 @@ def _translate_camera(
     # asis
     lk_asis = [
         'nin', 'e0', 'e1',
-        'nin_x', 'nin_y', 'nin_z',
-        'e0_x', 'e0_y', 'e0_z',
-        'e1_x', 'e1_y', 'e1_z',
+        ('nin_x', 'nin_y', 'nin_z'),
+        ('e0_x', 'e0_y', 'e0_z'),
+        ('e1_x', 'e1_y', 'e1_z'),
         ('outline_x0', 'outline_x1'),
     ]
 

tofu/physics_tools/electrons/distribution/__init__.py

@@ -1,5 +1,6 @@
 from ._runaway_growth import get_RE_critical_dreicer_electric_fields
 from ._runaway_growth import get_RE_growth_source_terms
 from ._distribution import main as get_distribution
+from ._distribution import get_dist1d_E
 from ._distribution_plot import main as plot_distribution
 from ._distribution_study import study_RE_vs_Maxwellian_distribution

tofu/physics_tools/electrons/distribution/_distribution.py

@@ -7,7 +7,9 @@
 import numpy as np
 import scipy.integrate as scpinteg
 import scipy.constants as scpct
+import scipy.stats as scpstats
 import astropy.units as asunits
+import datastock as ds
 
 
 from . import _distribution_check as _check
@@ -161,13 +163,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 +199,7 @@ def main(
             ddist=ddist,
             kdist=kdist,
             dcoords=dcoords,
-            version=version,
+            version=dfunc[kdist]['version'],
         )
 
         # store
@@ -355,7 +369,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)
@@ -372,6 +389,175 @@ def _get_velocity_par(ddist, kdist):
     return velocity_par
 
 
+# #####################################################
+# #####################################################
+#           Get Energy
+# #####################################################
+
+
+def get_dist1d_E(ddist, kdist, nbins=None):
+
+    # -----------
+    # check inputs
+    # -----------
+
+    nbins = int(ds._generic_check._check_var(
+        nbins, 'nbins',
+        types=(float, int),
+        default=ddist['coords']['x0']['data'].size,
+        sign='>3',
+    ))
+
+    # -----------
+    # prepare
+    # -----------
+
+    kcoords = tuple([
+        ddist['coords'][kk]['key'] for kk in ['x0', 'x1']
+        if ddist['coords'].get(kk) is not None
+    ])
+    shape = ddist['dist'][kdist]['dist']['data'].shape
+
+    # -----------
+    # (E, theta)
+    # -----------
+
+    if kcoords == ('E_eV', 'theta'):
+
+        assert ddist['coords']['x0']['units'] == 'eV'
+        E = ddist['coords']['x0']['data']
+
+        dist1d = scpinteg.trapezoid(
+            ddist['dist'][kdist]['dist']['data'],
+            x=ddist['coords']['x1']['data'],
+            axis=-1,
+        )
+        units = (
+            asunits.Unit(ddist['dist'][kdist]['dist']['units'])
+            * asunits.Unit(ddist['coords']['x1']['units'])
+        )
+
+    # -----------
+    # (E, pitch)
+    # -----------
+
+    elif kcoords == ('E_eV', 'pitch'):
+
+        assert ddist['coords']['x0']['units'] == 'eV'
+        E = ddist['coords']['x0']['data']
+
+        dist1d = scpinteg.trapezoid(
+            ddist['dist'][kdist]['dist']['data'],
+            x=ddist['coords']['x1']['data'],
+            axis=-1,
+        )
+        units = (
+            asunits.Unit(ddist['dist'][kdist]['dist']['units'])
+            * asunits.Unit(ddist['coords']['x1']['units'])
+        )
+
+    # -----------
+    # (p_par, p_perp)
+    # -----------
+
+    elif kcoords == ('p_par_norm', 'p_perp_norm'):
+
+        sli = (None,)*(len(shape)-2) + (slice(None),)*2
+        pnorm = np.sqrt(
+            ddist['coords']['x0']['data'][:, None]**2
+            + ddist['coords']['x1']['data'][None, :]**2
+        )[sli]
+        pnorm = np.broadcast_to(pnorm, shape)
+
+        # abs(velocity)
+        energy = _convert.convert_momentum_velocity_energy(
+            momentum_normalized=pnorm,
+        )['energy_kinetic_eV']
+
+        # nbins
+        E = np.linspace(np.min(energy)-1e-14, np.max(energy)+1e-14, nbins + 1)
+        dE = E[1] - E[0]
+
+        # binning
+        dist1d = scpstats.binned_statictics(
+            energy,
+            ddist['dist'][kdist]['dist']['data'],
+            bins=E,
+            statistic='sum',
+        ).statistic / dE
+        E = 0.5*(E[1:] + E[:-1])
+
+        units = (
+            asunits.Unit(ddist['dist'][kdist]['dist']['units'])
+            / asunits.unit('eV')
+        )
+
+    # -----------
+    # (v_par, v_perp)
+    # -----------
+
+    elif kcoords == ('v_par_norm', 'v_perp_norm'):
+
+        sli = (None,)*(len(shape)-2) + (slice(None),)*2
+        vv = np.sqrt(
+            ddist['coords']['x0']['data'][:, None]**2
+            + ddist['coords']['x1']['data'][None, :]**2
+        )[sli]
+        vv = np.broadcast_to(vv, shape)
+
+        # abs(velocity)
+        energy = _convert.convert_momentum_velocity_energy(
+            velocity_ms=vv,
+        )['energy_kinetic_eV']
+
+        # nbins
+        E = np.linspace(np.min(energy)-1e-14, np.max(energy)+1e-14, nbins + 1)
+        dE = E[1] - E[0]
+
+        # binning
+        dist1d = scpstats.binned_statictics(
+            energy,
+            ddist['dist'][kdist]['dist']['data'],
+            bins=E,
+            statistic='sum',
+        ).statistic / dE
+        E = 0.5*(E[1:] + E[:-1])
+
+        units = (
+            asunits.Unit(ddist['dist'][kdist]['dist']['units'])
+            / asunits.unit('eV')
+        )
+
+    # -----------
+    # (E,)
+    # -----------
+
+    elif kcoords == ('E_eV',):
+
+        E = ddist['coords']['x0']['data']
+        dist1d = ddist['dist'][kdist]['dist']['data']
+        units = ddist['dist'][kdist]['dist']['units']
+
+    else:
+        raise NotImplementedError(kcoords)
+
+    # ---------------
+    # abs() => v_par
+    # ---------------
+
+    ddist1d = {
+        'E': {
+            'data': E,
+            'units': 'eV',
+        },
+        'dist': {
+            'data': dist1d,
+            'units': units,
+        },
+    }
+
+    return ddist1d
+
 # #####################################################
 # #####################################################
 #           Integrate numerically
@@ -401,6 +587,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