Tmfp

cstool/compile.py

@@ -147,3 +147,11 @@ def compute_tcs_icdf(f, a, b, P, sampling=100000):
     if cf[-1]==0:
         return 0 * x.units*y.units, np.zeros_like(P) * x.units
     return cf[-1] * x.units*y.units, np.interp(P, cf/cf[-1], x) * x.units
+
+def compute_tcs(f, a, b, P, sampling=100000):
+    x = np.linspace(a, b.to(a.units), sampling) * a.units
+    y = f(x)
+    cf = np.r_[0, np.cumsum((x[1:] - x[:-1]) * (y[:-1] + y[1:]) / 2.0)]
+    if cf[-1]==0:
+        return 0 * x.units*y.units
+    return cf[-1] * x.units*y.units

cstool/inelastic.py

@@ -39,6 +39,37 @@ def L_Kieft(K, w0, F):
     return np.maximum(0, (w0 < 50 * units.eV) * L1
                       + (w0 > 50 * units.eV) * L2)
 
+def L_dv1(K, w0, F):
+    """Computes electron cross-sections for inelastic scattering from
+    optical data. Model is conceptually somewhere between L_Kieft and L_Ashley:
+    L1 is a Fermi-corrected version of Ashley without the factor 3/2
+    rescale by Kieft; L2 is the same as in Kieft.
+
+    :param K:
+        Kinetic energy of electron.
+    :param w0:
+        ω₀ - transition energy
+    :param F:
+        Fermi energy
+    """
+
+    a = (w0 / K).magnitude
+    b = (F / K).magnitude
+    s = sqrt(1 - 2*a, where = (a <= .5), out = np.zeros(a.shape))
+
+    L1_range = (a > 0) * (a < .5) * (a - s < 1 - 2*b)
+    L2_range = (a > 0) * (a < 1 - b)
+
+    # Calculate L1
+    wm = (1 + a - s)/2
+    wp = np.minimum((1 + a + s)/2, 1 - b)
+    L1 = log((wp - a) * wm / (wp * (wm - a)), where = L1_range, out = np.zeros(a.shape))
+
+    # Calculate L2
+    L2 = -log(a, where = L2_range, out = np.zeros(a.shape))
+
+    return np.maximum(0, (w0 < 50 * units.eV) * L1
+                      + (w0 > 50 * units.eV) * L2)
 
 def L_Ashley_w_ex(K, w0, _):
     a = w0 / K
@@ -54,7 +85,8 @@ def L_Ashley_wo_ex(K, w0, _):
 methods = {
     'Kieft': L_Kieft,
     'Ashley_w_ex': L_Ashley_w_ex,
-    'Ashley_wo_ex': L_Ashley_wo_ex
+    'Ashley_wo_ex': L_Ashley_wo_ex,
+    'dv1': L_dv1
 }
 
 

examples/cs.py

@@ -6,6 +6,7 @@
 from cstool.phonon import phonon_cs_fn
 from cstool.inelastic import inelastic_cs_fn
 from cstool.compile import compute_tcs_icdf
+from cstool.compile import compute_tcs
 from cstool.ionization import ionization_shells, outer_shell_energies, \
                               loglog_interpolate as ion_loglog_interp
 from cslib.noodles import registry
@@ -23,6 +24,12 @@ def integrant(theta):
 
     return compute_tcs_icdf(integrant, 0*units('rad'), np.pi*units('rad'), P)
 
+def compute_elastic_tcs(dcs, P):
+    def integrant(theta):
+        return dcs(theta) * 2 * np.pi * np.sin(theta)
+
+    return compute_tcs(integrant, 0*units('rad'), np.pi*units('rad'), P)
+
 
 def compute_inelastic_tcs_icdf(dcs, P, K0, K, max_interval):
     def integrant(w):
@@ -33,6 +40,15 @@ def integrant(w):
             int(np.ceil((K - K0) / max_interval))
             ]))
 
+def compute_inelastic_tcs(dcs, P, K0, K, max_interval):
+    def integrant(w):
+        return dcs(w)
+
+    return compute_tcs(integrant, K0, K, P,
+        sampling = np.min([100000,
+            int(np.ceil((K - K0) / max_interval))
+            ]))
+
 
 if __name__ == "__main__":
     import argparse
@@ -193,4 +209,74 @@ def dcs(w):
     ionization_osi = ionization_grp.create_dataset("outer_shells", data=s.elf_file.get_outer_shells().to('eV'))
     ionization_osi.attrs['units'] = 'eV'
 
+    # electron range
+    e_ran = np.logspace(-2, 3, 129) * units.eV
+    p_ran = np.linspace(0.0, 1.0, 1024)
+
+    electron_range_grp = outfile.create_group("electron_range")
+    ran_energies = electron_range_grp.create_dataset("energy", data=e_ran.to('eV'))
+    ran_energies.attrs['units'] = 'eV'
+    ran_range = electron_range_grp.create_dataset("range", e_ran.shape)
+    ran_range.attrs['units'] = 'm'
+
+    ran_tmfp = np.zeros(e_ran.shape) * units.m;
+    tmfpmax = 0. * units.m
+    print("# Computing transport mean free path.")
+    for i, K in enumerate(e_ran):
+        def dcs(theta):
+            return elastic_cs_fn(theta, K) * (1 - np.cos(theta)) * s.rho_n
+        inv_tmfp = compute_elastic_tcs(dcs, p_ran)
+        tmfp = 1./inv_tmfp
+        if ( K >= s.band_structure.barrier):
+            # make sure the tmfp is monotonously increasing for energies above
+            # the vacuum potential barrier
+            if (tmfp < tmfpmax):
+                tmfp = tmfpmax
+            else:
+                tmfpmax = tmfp
+        else:
+            tmfpmax = tmfp
+        ran_tmfp[i] = tmfp
+        print('.', end='', flush=True)
+    print()
+
+    rangemax = 0. * units.m
+    tmprange = np.zeros(e_ran.shape) * units.m
+    print("# Computing inelastic electron range.")
+    for i, K in enumerate(e_ran):
+        ran_energies[i] = K.to('eV')
+        if (K < s.band_structure.barrier):
+            ran_range[i] = 0. * units.m
+        else:
+            w0_max = K-s.band_structure.fermi # it is not possible to lose so
+            # much energy that the primary electron ends up below the Fermi
+            # level in an inelastic scattering event
+
+            def dcs(w):
+                return inelastic_cs_fn(s)(K, w) * s.rho_n
+            tcs = compute_inelastic_tcs(dcs, p_inel,
+                s.elf_file.get_min_energy(), w0_max,
+                s.elf_file.get_min_energy_interval())
+            tmprange[i] = 1/tcs
+            for j, omega in enumerate(e_ran):
+                if (omega < K/2 and omega < K-s.band_structure.fermi):
+                    if j==0:
+                        omega_old = 0. * units.eV
+                    else:
+                        omega_old = e_ran[j-1]
+                    cdcsvalue = inelastic_cs_fn(s)(K,omega) * s.rho_n * (omega - omega_old)
+                    index = sum(e_ran <= (K - omega)) - 1 # get last energy index smaller
+                    # than the energy remaining after one energy loss event
+                    tmprange[i] += tmprange[index] * cdcsvalue / tcs
+                else:
+                    break
+            rangevalue = tmprange[i]
+            if (2 * ran_tmfp[i] < tmprange[i]):
+                rangevalue = np.sqrt(2 * ran_tmfp[i] * tmprange[i])
+            if (rangevalue > rangemax):
+                rangemax = rangevalue
+            ran_range[i] = rangemax.to('m')
+        print('.', end='', flush=True)
+    print()
+
     outfile.close()