import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.special import expit, gamma
from functools import lru_cache
import argparse
from openpyxl import Workbook

try:
    from tklib.tktransport.tkfermi_integral import FermiIntegral_fast
except:
    FermiIntegral_fast = None

"""
熱電輸送特性計算モジュール

概要:
ゼーベック係数、ローレンツ因子、電気伝導度、移動度、電子熱伝導度、パワーファクター、ZT因子を計算し、
キャリア濃度依存性をプロット・Excelファイルに保存します。

詳細説明:
このモジュールは、3次元放物線バンドモデルと緩和時間近似（tau(E) ∝ E^(r-1/2)）に基づき、
フェルミ・ディラック積分を用いた厳密な輸送積分によって熱電輸送特性を評価します。
コマンドライン引数を通じて、温度、有効質量、参照移動度、格子熱伝導率、キャリア種別、散乱因子rのリストなどを設定できます。

関連リンク:
:doc:`Seebeck_ZT_usage`
"""

# =========================
# constants (SI)
# =========================
kB = 1.380649e-23          # J/K
e  = 1.602176634e-19       # C
h  = 6.62607015e-34        # J s
hbar = h / (2.0 * np.pi)
m0 = 9.1093837015e-31      # kg


r"""
散乱因子 r は tau(E) ∝ E^(r-1/2) で定義
"""



figsize = (10, 10)

# ============================================================
# Fermi–Dirac integral (NO 1/Gamma(j+1))
#   Fj(j, xi) = ∫_0^∞ x^j / (1+exp(x-xi)) dx
# ============================================================
@lru_cache(maxsize=40000)
def Fj_buildin(j, xi):
    """
    概要:
    フェルミ・ディラック積分 `F_j(xi) = ∫_0^∞ x^j / (1+exp(x-xi)) dx` を計算します。

    詳細説明:
    SciPyの`quad`関数を使用して数値積分を実行します。結果は`lru_cache`でキャッシュされ、
    同じ引数での再計算を回避し、パフォーマンスを向上させます。
    この関数は、ガンマ関数による正規化（1/Gamma(j+1)）は含まれていません。

    :param j: 積分の次数。
    :type j: float
    :param xi: 還元フェルミ準位 (reduced Fermi level)。
    :type xi: float
    :returns: フェルミ・ディラック積分 `F_j(xi)` の値。
    :rtype: float
    """
    j = float(j)
    xi = float(np.round(xi, 12))
    integrand = lambda x: (x**j) * expit(xi - x)
    val, _ = quad(integrand, 0.0, np.inf, epsabs=1e-10, epsrel=1e-10, limit=400)
    return val

fermi_integral_type = "build-in"
Fj = Fj_buildin
if FermiIntegral_fast:
    fermi_integral_type = "tkfermi_integral.FermiIntegral_fast"
    Fj = lambda j, xj: FermiIntegral_fast(xj, j, epsabs=1.0e-10, epsrel=1.0e-8, limit=50, use_cache = True)


def Nc_3D(T, m_eff=1.0):
    """
    概要:
    3次元放物線バンドにおける有効状態密度 `Nc` を計算します。

    詳細説明:
    以下の公式を用いて有効状態密度を計算します。
    `Nc = 2 * (2π m* kB T / h^2)^(3/2)`
    ここで、`m*` は有効質量、`kB` はボルツマン定数、`T` は温度、`h` はプランク定数です。

    :param T: 温度 [K]。
    :type T: float
    :param m_eff: 有効質量の比率 (m*/m0)。`m0` は自由電子質量です。デフォルトは1.0。
    :type m_eff: float, optional
    :returns: 有効状態密度 `Nc` [m^-3]。
    :rtype: float
    """
    mstar = m_eff * m0
    return 2.0 * (2.0 * np.pi * mstar * kB * T / (h**2))**1.5


def electron_density_from_xi(xi, T, m_eff=1.0):
    """
    概要:
    還元フェルミ準位 `xi` から電子濃度 `n` を計算します。

    詳細説明:
    電子濃度は以下の式で計算されます。
    `n [m^-3] = Nc * F_{1/2}^{std}(xi)`
    ここで、`Nc` は有効状態密度、`F^{std}_{1/2}(xi) = Fj(1/2, xi) / Gamma(3/2)` は
    標準化されたフェルミ・ディラック積分です。

    :param xi: 還元フェルミ準位。
    :type xi: float
    :param T: 温度 [K]。
    :type T: float
    :param m_eff: 有効質量の比率 (m*/m0)。デフォルトは1.0。
    :type m_eff: float, optional
    :returns: 電子濃度 `n` [m^-3]。
    :rtype: float
    """
    Nc = Nc_3D(T, m_eff=m_eff)
    F12_std = Fj(0.5, xi) / gamma(1.5)
    return Nc * F12_std


# ============================================================
# Transport integrals for tau(E) ∝ E^(r-1/2)
#   I0 = (r+1) F_r
#   I1 = (r+2) F_{r+1} - xi (r+1) F_r
#   I2 = (r+3) F_{r+2} - 2 xi (r+2) F_{r+1} + xi^2 (r+1) F_r
# ============================================================
def transport_I0(xi, r):
    """
    概要:
    輸送積分 `I0` を計算します。

    詳細説明:
    `tau(E) ∝ E^(r-1/2)` の緩和時間に対する輸送積分の一つである `I0` を計算します。
    `I0 = (r + 1) * F_r(xi)`

    :param xi: 還元フェルミ準位。
    :type xi: float
    :param r: 散乱因子。緩和時間のエネルギー依存性 `E^(r-1/2)` を定義します。
    :type r: float
    :returns: 輸送積分 `I0` の値。
    :rtype: float
    """
    return (r + 1.0) * Fj(r, xi)


def transport_I1(xi, r):
    """
    概要:
    輸送積分 `I1` を計算します。

    詳細説明:
    `tau(E) ∝ E^(r-1/2)` の緩和時間に対する輸送積分の一つである `I1` を計算します。
    `I1 = (r + 2) * F_{r+1}(xi) - xi * (r + 1) * F_r(xi)`

    :param xi: 還元フェルミ準位。
    :type xi: float
    :param r: 散乱因子。緩和時間のエネルギー依存性 `E^(r-1/2)` を定義します。
    :type r: float
    :returns: 輸送積分 `I1` の値。
    :rtype: float
    """
    return (r + 2.0) * Fj(r + 1.0, xi) - xi * (r + 1.0) * Fj(r, xi)


def transport_I2(xi, r):
    """
    概要:
    輸送積分 `I2` を計算します。

    詳細説明:
    `tau(E) ∝ E^(r-1/2)` の緩和時間に対する輸送積分の一つである `I2` を計算します。
    `I2 = (r + 3) * F_{r+2}(xi) - 2 * xi * (r + 2) * F_{r+1}(xi) + xi^2 * (r + 1) * F_r(xi)`

    :param xi: 還元フェルミ準位。
    :type xi: float
    :param r: 散乱因子。緩和時間のエネルギー依存性 `E^(r-1/2)` を定義します。
    :type r: float
    :returns: 輸送積分 `I2` の値。
    :rtype: float
    """
    return (
        (r + 3.0) * Fj(r + 2.0, xi)
        - 2.0 * xi * (r + 2.0) * Fj(r + 1.0, xi)
        + (xi**2) * (r + 1.0) * Fj(r, xi)
    )


def transport_prefactor_sigma(T, m_eff=1.0):
    """
    概要:
    電気伝導度 `sigma` の計算における前因子 `A_sigma` を計算します。

    詳細説明:
    電気伝導度は `sigma = A_sigma * tau_pref * (kB T)^(r+1) * I0` の形で表されます。
    この関数は、`tau(E) = tau_pref * E^(r-1/2)` と定義される緩和時間において、
    以下の物理定数を含む前因子 `A_sigma` を計算します。
    `A_sigma [SI] = e^2 * 2^(3/2) * sqrt(m*) / (3 π^2 ħ^3)`
    ここで、`e` は素電荷、`m*` は有効質量、`ħ` は換算プランク定数です。

    :param T: 温度 [K]。
    :type T: float
    :param m_eff: 有効質量の比率 (m*/m0)。デフォルトは1.0。
    :type m_eff: float, optional
    :returns: 電気伝導度の前因子 `A_sigma` [SI単位]。
    :rtype: float
    """
    mstar = m_eff * m0
    return e**2 * (2.0**1.5) * np.sqrt(mstar) / (3.0 * np.pi**2 * hbar**3)


def sigma_from_xi_tau_pref(xi, r, T, tau_pref, m_eff=1.0):
    """
    概要:
    還元フェルミ準位 `xi`、散乱因子 `r`、温度 `T`、緩和時間の前因子 `tau_pref` から
    電気伝導度 `sigma` を計算します。

    詳細説明:
    電気伝導度は以下の式で計算されます。
    `sigma = A_sigma * tau_pref * (kB T)^(r + 1.0) * I0`
    ここで、`A_sigma` は `transport_prefactor_sigma` で計算される前因子、
    `I0` は `transport_I0` で計算される輸送積分です。

    :param xi: 還元フェルミ準位。
    :type xi: float
    :param r: 散乱因子。
    :type r: float
    :param T: 温度 [K]。
    :type T: float
    :param tau_pref: 緩和時間の前因子 `tau_pref` [s / J^(r-1/2)]。
    :type tau_pref: float
    :param m_eff: 有効質量の比率 (m*/m0)。デフォルトは1.0。
    :type m_eff: float, optional
    :returns: 電気伝導度 `sigma` [S/m]。
    :rtype: float
    """
    A = transport_prefactor_sigma(T, m_eff=m_eff)
    I0 = transport_I0(xi, r)
    return A * tau_pref * (kB * T)**(r + 1.0) * I0


def seebeck_from_xi_transport(xi, r, carrier="electron"):
    """
    概要:
    還元フェルミ準位 `xi` と散乱因子 `r` からゼーベック係数 `S` を計算します。

    詳細説明:
    ゼーベック係数は厳密な輸送積分を用いて以下の式で計算されます。
    `S = (kB/q) * I1/I0`
    `q` はキャリアの種類によって異なり、電子の場合は `-e`、正孔の場合は `+e` となります。

    :param xi: 還元フェルミ準位。
    :type xi: float
    :param r: 散乱因子。`tau(E) ∝ E^(r-1/2)` で定義されます。
    :type r: float
    :param carrier: キャリアの種類。"electron"または"hole"を指定します。デフォルトは"electron"。
    :type carrier: str, optional
    :returns: ゼーベック係数 `S` [V/K]。
    :rtype: float
    :raises ValueError: `carrier` が "electron" または "hole" 以外の場合。
    """
    if carrier.lower().startswith("e"):
        q = -e
    elif carrier.lower().startswith("h"):
        q = +e
    else:
        raise ValueError("carrier must be 'electron' or 'hole'")

    I0 = transport_I0(xi, r)
    I1 = transport_I1(xi, r)
    return (kB / q) * (I1 / I0)


def lorenz_from_xi_transport(xi, r):
    """
    概要:
    還元フェルミ準位 `xi` と散乱因子 `r` からローレンツ因子 `L` を計算します。

    詳細説明:
    ローレンツ因子は厳密な輸送積分を用いて以下の式で計算されます。
    `L = (kB/e)^2 * ( I2/I0 - (I1/I0)^2 )`

    :param xi: 還元フェルミ準位。
    :type xi: float
    :param r: 散乱因子。
    :type r: float
    :returns: ローレンツ因子 `L` [V^2/K^2 (または W ohm / K^2)]。
    :rtype: float
    """
    I0 = transport_I0(xi, r)
    I1 = transport_I1(xi, r)
    I2 = transport_I2(xi, r)
    return (kB / e)**2 * (I2 / I0 - (I1 / I0)**2)


def mobility_from_sigma_n(sigma, n):
    """
    概要:
    電気伝導度 `sigma` と電子濃度 `n` から移動度 `mu` を計算します。

    詳細説明:
    電気伝導度 `sigma`、電子濃度 `n`、素電荷 `e`、移動度 `mu` の間に成り立つ関係式
    `sigma = n * e * mu` を用いて移動度を計算します。
    したがって、`mu = sigma / (n * e)` となります。

    :param sigma: 電気伝導度 [S/m]。
    :type sigma: float
    :param n: 電子濃度 [m^-3]。
    :type n: float
    :returns: 移動度 `mu` [m^2/V/s]。
    :rtype: float
    """
    return sigma / (n * e)


def tau_pref_from_mu_ref(mu_ref_cm2_Vs, xi_ref, r, T, m_eff=1.0):
    """
    概要:
    参照移動度 `mu_ref` と参照還元フェルミ準位 `xi_ref` に基づいて、
    緩和時間の前因子 `tau_pref` を決定します。

    詳細説明:
    緩和時間 `tau(E)` が `tau(E) = tau_pref * E^(r-1/2)` で定義されるとき、
    指定された `xi_ref` と `T` における移動度が `mu_ref` となるように
    `tau_pref` を逆算します。
    参照移動度 `mu_ref_cm2_Vs` は cm^2/V/s 単位で与えられますが、計算ではSI単位に変換されます。

    :param mu_ref_cm2_Vs: 参照移動度 [cm^2/V/s]。
    :type mu_ref_cm2_Vs: float
    :param xi_ref: 参照還元フェルミ準位。この `xi` で参照移動度が再現されます。
    :type xi_ref: float
    :param r: 散乱因子。
    :type r: float
    :param T: 温度 [K]。
    :type T: float
    :param m_eff: 有効質量の比率 (m*/m0)。デフォルトは1.0。
    :type m_eff: float, optional
    :returns: 緩和時間の前因子 `tau_pref` [s / J^(r-1/2)]。
    :rtype: float
    """
    mu_ref_SI = mu_ref_cm2_Vs * 1e-4  # cm^2/V/s -> m^2/V/s
    n_ref = electron_density_from_xi(xi_ref, T=T, m_eff=m_eff)
    sigma_ref = n_ref * e * mu_ref_SI

    A = transport_prefactor_sigma(T, m_eff=m_eff)
    I0_ref = transport_I0(xi_ref, r)

    tau_pref = sigma_ref / (A * (kB * T)**(r + 1.0) * I0_ref)
    return tau_pref


def tau_at_energy(E_J, tau_pref, r):
    """
    概要:
    特定のエネルギー `E_J` における緩和時間 `tau(E_J)` を計算します。

    詳細説明:
    緩和時間は以下の式で定義されます。
    `tau(E) = tau_pref * E^(r-1/2)`
    `E_J` が0以下の場合は `np.nan` を返します。

    :param E_J: エネルギー [J]。
    :type E_J: float
    :param tau_pref: 緩和時間の前因子 [s / J^(r-1/2)]。
    :type tau_pref: float
    :param r: 散乱因子。
    :type r: float
    :returns: エネルギー `E_J` における緩和時間 `tau(E_J)` [s]。
    :rtype: float
    """
    if E_J <= 0.0:
        return np.nan
    return tau_pref * (E_J**(r - 0.5))


def equivalent_l0_from_tau_pref(tau_pref, r, m_eff=1.0):
    """
    概要:
    緩和時間の前因子 `tau_pref` から等価な平均自由行程 `l0` を計算します。

    詳細説明:
    一部のスライドでの慣例的な定義
    `tau(E) = sqrt(m*/2) * [ l0(T) * e_chg^(-r) ] * E^(r-1/2)`
    から、`l0(T)` を逆算します。
    したがって、`l0(T) = tau_pref * sqrt(2/m*) * e_chg^r` となります。

    :param tau_pref: 緩和時間の前因子 [s / J^(r-1/2)]。
    :type tau_pref: float
    :param r: 散乱因子。
    :type r: float
    :param m_eff: 有効質量の比率 (m*/m0)。デフォルトは1.0。
    :type m_eff: float, optional
    :returns: 等価な平均自由行程 `l0` [m]。
    :rtype: float
    """
    mstar = m_eff * m0
    return tau_pref * np.sqrt(2.0 / mstar) * (e**r)


def save_results_to_excel(outfile, meta_rows, data_rows):
    """
    概要:
    計算結果とメタデータをExcelファイルに保存します。

    詳細説明:
    指定された出力ファイル名で新しいExcelワークブックを作成し、
    "metadata"と"transport_vs_Ne"の2つのシートにデータを書き込みます。
    "metadata"シートには計算条件などのメタデータが保存され、
    "transport_vs_Ne"シートにはキャリア濃度に対する輸送特性の計算結果が保存されます。

    :param outfile: 出力するExcelファイル名。例: "results.xlsx"。
    :type outfile: str
    :param meta_rows: メタデータを含む行のリスト。各要素はExcelの1行に対応するリストです。
    :type meta_rows: list[list]
    :param data_rows: 計算結果データを含む行のリスト。各要素はExcelの1行に対応するリストです。
    :type data_rows: list[list]
    :returns: なし
    :rtype: None
    """
    wb = Workbook()
    ws_meta = wb.active
    ws_meta.title = "metadata"

    for row in meta_rows:
        ws_meta.append(row)

    ws = wb.create_sheet("transport_vs_Ne")
    ws.append([
        "r",
        "xi",
        "tau_pref_s_per_J^(r-1/2)",
        "tau_at_kBT_s",
        "l0_equiv_m",
        "Ne_m^-3",
        "Ne_cm^-3",
        "S_V_per_K",
        "S_uV_per_K",
        "L_Wohm_per_K2",
        "sigma_S_per_m",
        "mu_m2_per_Vs",
        "mu_cm2_per_Vs",
        "kappa_e_W_per_mK",
        "PF_W_per_mK2",
        "PF_mW_per_mK2",
        "ZT",
    ])
    for row in data_rows:
        ws.append(row)

    wb.save(outfile)


def main():
    """
    概要:
    スクリプトの主要な実行フローを定義します。コマンドライン引数を解析し、
    熱電輸送特性を計算・プロット・保存します。

    詳細説明:
    1.  `argparse` モジュールを使用してコマンドライン引数（温度、有効質量、参照移動度、
        格子熱伝導率、キャリア種別、還元フェルミ準位の範囲と点数、散乱因子`r`のリスト、
        出力ファイル名など）を解析します。
    2.  指定されたパラメータに基づいて、各散乱因子 `r` に対して輸送特性（電子濃度、
        ゼーベック係数、ローレンツ因子、電気伝導度、移動度、電子熱伝導度、パワーファクター、ZT因子）を計算します。
    3.  計算されたデータは `save_results_to_excel` 関数を使用してExcelファイルに保存されます。
    4.  計算結果は matplotlib を用いてキャリア濃度に対する各特性のグラフとしてプロットされ、表示されます。

    :returns: なし
    :rtype: None
    """
    parser = argparse.ArgumentParser(
        description=(
            "S, L, sigma, mu, kappa_e, PF, ZT vs Ne using exact transport integrals "
            "for tau(E) ∝ E^(r-1/2) in a 3D parabolic band."
        )
    )
    parser.add_argument("--T", type=float, default=300.0, help="Temperature [K]")
    parser.add_argument("--m_eff", type=float, default=1.0, help="Effective mass m*/m0")
    parser.add_argument("--mu_ref", type=float, default=10.0,
                        help="Reference mobility at xi_ref [cm^2/V/s]")
    parser.add_argument("--xi_ref", type=float, default=0.0,
                        help="Reference reduced Fermi level xi where mu_ref is reproduced")
    parser.add_argument("--klatt", type=float, default=5.0,
                        help="Lattice thermal conductivity [W/m/K] (default: 5)")
    parser.add_argument("--carrier", type=str, default="electron",
                        choices=["electron", "hole"], help="Carrier type")
    parser.add_argument("--xi_min", type=float, default=-5.0, help="Minimum xi")
    parser.add_argument("--xi_max", type=float, default=40.0, help="Maximum xi")
    parser.add_argument("--nxi", type=int, default=181, help="Number of xi points")
    parser.add_argument("--r_list", type=str, default="2.0,1.5,1.0,0.5,0.0,-0.5",
                        help="Comma-separated r list")
    parser.add_argument("--outfile", type=str, default="transport_vs_Ne.xlsx",
                        help="Excel output filename")
    args = parser.parse_args()

    T = args.T
    m_eff = args.m_eff
    mu_ref_cm2_Vs = args.mu_ref
    xi_ref = args.xi_ref
    klatt = args.klatt
    carrier = args.carrier
    r_list = [float(x.strip()) for x in args.r_list.split(",")]
    xi_range = np.linspace(args.xi_min, args.xi_max, args.nxi)
    outfile = args.outfile

    print("Convention: tau(E) ∝ E^(r-1/2)")
    print(f"T         = {T:g} K")
    print(f"m_eff     = {m_eff:g} m0")
    print(f"mu_ref    = {mu_ref_cm2_Vs:g} cm^2/V/s")
    print(f"xi_ref    = {xi_ref:g}")
    print(f"klatt     = {klatt:g} W/m/K")
    print(f"carrier   = {carrier}")
    print(f"outfile   = {outfile}")
    print()

    all_results = []

    for r in r_list:
        tau_pref = tau_pref_from_mu_ref(mu_ref_cm2_Vs, xi_ref, r, T, m_eff=m_eff)
        tau_kBT = tau_at_energy(kB * T, tau_pref, r)
        l0_equiv = equivalent_l0_from_tau_pref(tau_pref, r, m_eff=m_eff)

        print(f"r = {r:g}")
        print(f"  tau_pref        = {tau_pref:.6e}  [s / J^({r:g}-1/2)]")
        print(f"  tau(E=kBT)      = {tau_kBT:.6e}  [s]")
        print(f"  equivalent l0   = {l0_equiv:.6e}  [m]  (slide convention)")

        n_list = []
        S_list = []
        L_list = []
        sig_list = []
        mu_list = []
        ke_list = []
        PF_list = []
        ZT_list = []

        for xi in xi_range:
            n = electron_density_from_xi(xi, T=T, m_eff=m_eff)                # [m^-3]
            sigma = sigma_from_xi_tau_pref(xi, r, T, tau_pref, m_eff=m_eff)   # [S/m]
            mu_SI = mobility_from_sigma_n(sigma, n)                            # [m^2/V/s]
            mu_cm2_Vs = mu_SI * 1e4

            S = seebeck_from_xi_transport(xi, r, carrier=carrier)              # [V/K]
            L = lorenz_from_xi_transport(xi, r)                                # [V^2/K^2]
            kappa_e = L * sigma * T                                            # [W/m/K]
            PF = S**2 * sigma                                                   # [W/m/K^2]
            ZT = PF * T / (kappa_e + klatt)

            n_list.append(n)
            S_list.append(S)
            L_list.append(L)
            sig_list.append(sigma)
            mu_list.append(mu_cm2_Vs)
            ke_list.append(kappa_e)
            PF_list.append(PF)
            ZT_list.append(ZT)

            all_results.append([
                r,
                xi,
                tau_pref,
                tau_kBT,
                l0_equiv,
                n,
                n / 1e6,
                S,
                S * 1e6,
                L,
                sigma,
                mu_SI,
                mu_cm2_Vs,
                kappa_e,
                PF,
                PF * 1e3,
                ZT,
            ])

        # 配列化して保持
        # この行は、ループ内でリストに追加された要素を再度リストに代入するもので、
        # 結果的にはall_resultsのリスト構造に影響を与えません。
        # 論理的な意味は薄いですが、既存コードの変更禁止ルールに従い保持します。
        all_results[-len(xi_range):] = all_results[-len(xi_range):]

    # ---------- save Excel before plotting ----------
    meta_rows = [
        ["Convention", "tau(E) ∝ E^(r-1/2)"],
        ["T_K", T],
        ["m_eff_over_m0", m_eff],
        ["mu_ref_cm2_per_Vs", mu_ref_cm2_Vs],
        ["xi_ref", xi_ref],
        ["klatt_W_per_mK", klatt],
        ["carrier", carrier],
        ["r_list", ",".join(str(r) for r in r_list)],
        ["xi_min", args.xi_min],
        ["xi_max", args.xi_max],
        ["nxi", args.nxi],
    ]
    save_results_to_excel(outfile, meta_rows, all_results)
    print(f"\nSaved calculation results to [{outfile}]\n")

    # ---------- plot ----------
    fig, axes = plt.subplots(4, 2, figsize=figsize)
    axS   = axes[0, 0]
    axL   = axes[0, 1]
    axSig = axes[1, 0]
    axMu  = axes[1, 1]
    axKe  = axes[2, 0]
    axPF  = axes[2, 1]
    axZT  = axes[3, 0]
    axBlank = axes[3, 1]

    for r in r_list:
        tau_pref = tau_pref_from_mu_ref(mu_ref_cm2_Vs, xi_ref, r, T, m_eff=m_eff)

        n_list = []
        S_list = []
        L_list = []
        sig_list = []
        mu_list = []
        ke_list = []
        PF_list = []
        ZT_list = []

        for xi in xi_range:
            n = electron_density_from_xi(xi, T=T, m_eff=m_eff)
            sigma = sigma_from_xi_tau_pref(xi, r, T, tau_pref, m_eff=m_eff)
            mu_SI = mobility_from_sigma_n(sigma, n)
            S = seebeck_from_xi_transport(xi, r, carrier=carrier)
            L = lorenz_from_xi_transport(xi, r)
            kappa_e = L * sigma * T
            PF = S**2 * sigma
            ZT = PF * T / (kappa_e + klatt)

            n_list.append(n)
            S_list.append(S)
            L_list.append(L)
            sig_list.append(sigma)
            mu_list.append(mu_SI * 1e4)
            ke_list.append(kappa_e)
            PF_list.append(PF)
            ZT_list.append(ZT)

        n_arr_cm3 = np.array(n_list) / 1e6
        S_arr_uVK = np.array(S_list) * 1e6
        L_arr     = np.array(L_list)
        sig_arr   = np.array(sig_list)
        mu_arr    = np.array(mu_list)
        ke_arr    = np.array(ke_list)
        PF_arr    = np.array(PF_list) * 1e3
        ZT_arr    = np.array(ZT_list)

        label = f"r = {r:g}"

        axS.plot(n_arr_cm3, S_arr_uVK, label=label)
        axL.plot(n_arr_cm3, L_arr, label=label)
        axSig.plot(n_arr_cm3, sig_arr, label=label)
        axMu.plot(n_arr_cm3, mu_arr, label=label)
        axKe.plot(n_arr_cm3, ke_arr, label=label)
        axPF.plot(n_arr_cm3, PF_arr, label=label)
        axZT.plot(n_arr_cm3, ZT_arr, label=label)

    for ax in [axS, axL, axSig, axMu, axKe, axPF, axZT]:
        ax.set_xscale("log")
        ax.grid(True, which="both", alpha=0.25)
        ax.set_xlabel(r"$N_e$ (cm$^{-3}$)")

    axS.axhline(0.0, linewidth=1.0, alpha=0.4)
    axS.set_ylabel(r"$S$ ($\mu$V/K)")
    axS.set_title("Seebeck coefficient")

    axL.set_ylabel(r"$L$ (V$^2$/K$^2$)")
    axL.set_title("Lorenz factor")

    axSig.set_ylabel(r"$\sigma$ (S/m)")
    axSig.set_title("Electrical conductivity")

    axMu.set_ylabel(r"$\mu$ (cm$^2$/V/s)")
    axMu.set_title(r"Mobility  ($\tau(E)\propto E^{r-1/2}$)")

    axKe.set_ylabel(r"$\kappa_e$ (W/m/K)")
    axKe.set_title("Electronic thermal conductivity")

    axPF.set_ylabel(r"$PF=S^2\sigma$ (mW/m/K$^2$)")
    axPF.set_title("Power factor")

    axZT.set_ylabel(r"$ZT$")
    axZT.set_title(rf"$ZT$  ($\kappa_{{latt}}$={klatt:g} W/m/K)")

    axBlank.axis("off")

    handles, labels = axS.get_legend_handles_labels()
    fig.legend(handles, labels, loc="upper center", ncol=min(len(r_list), 6))
    fig.suptitle(
        f"Transport vs Ne  (T={T:g} K, m*/m0={m_eff:g}, "
        f"mu_ref={mu_ref_cm2_Vs:g} cm$^2$/V/s at xi_ref={xi_ref:g}, "
        r"$\tau(E)\propto E^{r-1/2}$)",
        y=0.995
    )
    plt.tight_layout(rect=[0, 0, 1, 0.96])
    plt.show()


if __name__ == "__main__":
    main()