function [p,te] = p_WKB(E_list, W)
% caclulate ionization rate via WKB approximation; Input: work function W [eV], electric field E_list [V/Å]
load_const;

if isempty(E_list)
    p = [];
    te = [];
else
    for E_index = 1:numel(E_list)
        E = E_list(E_index);
        x_te = W/E;

        x = linspace(0, x_te - 1e-10,1e3);

        V = -x*E;

        p_loc = sqrt(2 * const.m_e * const.e * (V+W));

        gamma = trapz(x, 2 / const.hbar * p_loc * 1e-10);

        p(E_index) = exp(-gamma);
        if isnan(p(E_index))
            p(E_index) = 0;
        end
        te(E_index) = x_te;
    end
end
end