def mob_klassen(Nd, Na, Δn=1, T=298.16):
    """Return the mobility (cm2/Vs)
    given the dopings and excess carrier concentration
    DOI: 10.1109/IEDM.1990.237157"""
    s1 = 0.89233
    s2 = 0.41372
    s3 = 0.19778
    s4 = 0.28227
    s5 = 0.005978
    s6 = 1.80618
    s7 = 0.72169
    r1 = 0.7643
    r2 = 2.2999
    r3 = 6.5502
    r4 = 2.367
    r5 = -0.01552
    r6 = 0.6478
    fCW = 2.459
    fBH = 3.828
    mh_me = 1.258
    me_m0 = 1
    T = 298.16
    n0 = Nd
    p0 = Nd / 10000000000.0 ** 2
    n_i = 8310000000.0
    cA = 0.5
    cD = 0.21
    Nref_A = 7.2e+20
    Nref_D = 4e+20
    p = p0 + Δn
    n = n0 + Δn
    cc = p + n
    Za_Na = 1 + 1 / (cA + (Nref_A / Na) ** 2)
    Zd_Nd = 1 + 1 / (cD + (Nref_D / Nd) ** 2)
    Na_h = Za_Na * Na
    Nd_h = Zd_Nd * Nd
    phosphorus = [1414, 68.5, 9.2e+16, 0.711, 2.285]
    boron = [470.5, 44.9, 2.23e+17, 0.719, 2.247]
    boron_μmax = 470.5
    boron_μmin = 44.9
    boron_Nref_1 = 2.23e+17
    boron_α = 0.719
    boron_θ = 2.247
    phosphorus_μmax = 1414
    phosphorus_μmin = 68.5
    phosphorus_Nref_1 = 9.2e+16
    phosphorus_α = 0.711
    phosphorus_θ = 2.285
    μ_eN = phosphorus_μmax ** 2 / (phosphorus_μmax - phosphorus_μmin) * (T /
        300) ** (3 * phosphorus_α - 1.5)
    μ_hN = boron_μmax ** 2 / (boron_μmax - boron_μmin) * (T / 300) ** (3 *
        boron_α - 1.5)
    μ_ec = phosphorus_μmax * phosphorus_μmin / (phosphorus_μmax -
        phosphorus_μmin) * (300 / T) ** 0.5
    μ_hc = boron_μmax * boron_μmin / (boron_μmax - boron_μmin) * (300 / T
        ) ** 0.5
    μ_eD = μ_eN * (phosphorus_Nref_1 / Nd) ** phosphorus_α + μ_ec * (cc / Nd)
    μ_hA = μ_hN * (boron_Nref_1 / Na) ** boron_α + μ_hc * (cc / Na)
    Ne_sc = Na_h + Nd_h + p
    Nh_sc = Na_h + Nd_h + n
    PBHe = 1.36e+20 / cc * me_m0 * (T / 300) ** 2
    PBHh = 1.36e+20 / cc * mh_me * (T / 300) ** 2
    PCWe = 39700000000000.0 * (1 / (Zd_Nd ** 3 * (Nd_h + Na_h + p)) * (T / 
        300) ** 3) ** (2 / 3)
    PCWh = 39700000000000.0 * (1 / (Za_Na ** 3 * (Nd_h + Na_h + n)) * (T / 
        300) ** 3) ** (2 / 3)
    Pe = 1 / (fCW / PCWe + fBH / PBHe)
    Ph = 1 / (fCW / PCWh + fBH / PBHh)
    G_Pe = 1 - s1 / (s2 + (1 / me_m0 * 300 / T) ** s4 * Pe) ** s3 + s5 / ((
        me_m0 * 300 / T) ** s7 * Pe) ** s6
    G_Ph = 1 - s1 / (s2 + (1 / (me_m0 * mh_me) * T / 300) ** s4 * Ph
        ) ** s3 + s5 / ((me_m0 * mh_me * 300 / T) ** s7 * Ph) ** s6
    F_Pe = (r1 * Pe ** r6 + r2 + r3 / mh_me) / (Pe ** r6 + r4 + r5 / mh_me)
    F_Ph = (r1 * Ph ** r6 + r2 + r3 * mh_me) / (Ph ** r6 + r4 + r5 * mh_me)
    Ne_sc_eff = Nd_h + G_Pe * Na_h + p / F_Pe
    Nh_sc_eff = Na_h + G_Ph * Nd_h + n / F_Ph
    μ_eL = phosphorus_μmax * (300 / T) ** phosphorus_θ
    μ_hL = boron_μmax * (300 / T) ** boron_θ
    μe_Dah = μ_eN * Ne_sc / Ne_sc_eff * (phosphorus_Nref_1 / Ne_sc
        ) ** phosphorus_α + μ_ec * ((p + n) / Ne_sc_eff)
    μh_Dae = μ_hN * Nh_sc / Nh_sc_eff * (boron_Nref_1 / Nh_sc
        ) ** boron_α + μ_hc * ((p + n) / Nh_sc_eff)
    μe = 1 / (1 / μ_eL + 1 / μe_Dah)
    μh = 1 / (1 / μ_hL + 1 / μh_Dae)
    return μe, μh
