"""
概要:
このスクリプトは、フェルミ準位、キャリア濃度、温度、有効質量を様々にスキャンしながら、熱電材料のゼーベック係数を計算しプロットします。

詳細説明:
散乱因子 r の異なる値に対して、以下のスキャンモードをサポートしています。
- EF (フェルミ準位): フェルミ準位をスキャンし、ゼーベック係数を計算します。
- Ne (キャリア濃度): キャリア濃度をスキャンし、ゼーベック係数を計算します。
- T (温度): 温度をスキャンし、ゼーベック係数を計算します。このモードでは、フェルミ準位、キャリア濃度、または還元フェルミ準位のいずれかを固定できます。
- m_eff (有効質量): 有効質量をスキャンし、ゼーベック係数を計算します。このモードでも、フェルミ準位、キャリア濃度、または還元フェルミ準位のいずれかを固定できます。

計算はFermi–Dirac積分に基づき、キャリアタイプ（電子または正孔）を指定できます。
結果はCSVファイルとして保存することも可能です。

散乱因子 r はキャリアの緩和時間 tau(E) がエネルギー E の (r-1/2)乗に比例することで定義されます。

関連リンク:
Seebeck_scan_usage
"""
import argparse
import csv
from functools import lru_cache

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq
from scipy.special import expit, gamma

# =========================
# constants (SI)
# =========================
kB_SI = 1.380649e-23           # J/K
e_SI  = 1.602176634e-19        # C
h_SI  = 6.62607015e-34         # J s
m0_SI = 9.1093837015e-31       # kg
kB_eV = 8.617333262145e-5      # eV/K


# ============================================================
# Fermi–Dirac integral (NO 1/Gamma(j+1) prefactor)
#   Fj(j,xi) = ∫_0^∞ x^j / (1+exp(x-xi)) dx
# ============================================================
@lru_cache(maxsize=20000)
def Fj(j, xi):
    """
    概要:
        フェルミ・ディラック積分を計算します（ガンマ関数による前因子なし）。
    詳細説明:
        Fj(j, xi) = ∫_0^∞ x^j / (1+exp(x-xi)) dx を数値的に評価します。
        この関数は、引数 j と xi に対してフェルミ・ディラック積分を計算します。
        標準的な定義に含まれる 1/Gamma(j+1) の前因子は含まれません。
        計算結果は lru_cache デコレータによってキャッシュされ、同じ引数での再計算を防ぎ、
        パフォーマンスを向上させます。
        j が -1.0 以下の場合、ValueError を発生させます。
    引数:
        :param j: フェルミ・ディラック積分の次数。j は -1 より大きい必要があります。
        :type j: float
        :param xi: 還元フェルミ準位 (xi = E_F / (k_B T))。
        :type xi: float
    戻り値:
        :returns: 計算されたフェルミ・ディラック積分の値。
        :rtype: float
    例外:
        :raises ValueError: j が -1.0 以下の場合。
    """
    j = float(j)
    xi = float(np.round(xi, 12))
    if j <= -1.0:
        raise ValueError(f"Fj(j, xi) requires j > -1, got j={j}")

    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


def Nc_3D(T, m_eff=1.0):
    """
    概要:
        3次元における有効状態密度 Nc を計算します。
    詳細説明:
        温度 T と有効質量 m_eff に基づいて、3次元バンドの有効状態密度 Nc を計算します。
        Nc は以下の式で計算されます。
        Nc = 2 * (2 * pi * m_eff * m0_SI * kB_SI * T / (h_SI**2)) ** 1.5
    引数:
        :param T: 温度 [K]。
        :type T: float
        :param m_eff: 有効質量比 m*/m0。m0 は自由電子質量です。
        :type m_eff: float
    戻り値:
        :returns: 計算された有効状態密度 Nc [m^-3]。
        :rtype: float
    """
    mstar = m_eff * m0_SI
    return 2.0 * (2.0 * np.pi * mstar * kB_SI * T / (h_SI**2)) ** 1.5


def electron_density_from_xi(xi, T, m_eff=1.0):
    """
    概要:
        還元フェルミ準位 xi から電子濃度を計算します。
    詳細説明:
        還元フェルミ準位 xi、温度 T、有効質量 m_eff を用いて電子濃度 n を計算します。
        電子濃度 n は以下の式で与えられます:
        n = Nc * Fj(0.5, xi) / Gamma(1.5)
        ここで、Nc は有効状態密度、Fj(0.5, xi) はガンマ関数による前因子を含まないフェルミ・ディラック積分です。
    引数:
        :param xi: 還元フェルミ準位。
        :type xi: float
        :param T: 温度 [K]。
        :type T: float
        :param m_eff: 有効質量比 m*/m0。
        :type m_eff: float
    戻り値:
        :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


def xi_from_electron_density(n_m3, T, m_eff=1.0, xi_low=-40.0, xi_high=120.0):
    """
    概要:
        与えられた電子濃度から還元フェルミ準位 xi を計算します。
    詳細説明:
        この関数は、electron_density_from_xi 関数を解いて xi を見つけるために、
        Brentのアルゴリズム (scipy.optimize.brentq) を使用します。
        初期の探索範囲 [xi_low, xi_high] で解が見つからない場合、
        自動的に探索範囲を拡張しようと試みます。
        電子濃度が0以下の場合、ValueError を発生させます。
    引数:
        :param n_m3: 電子濃度 [m^-3]。正の値である必要があります。
        :type n_m3: float
        :param T: 温度 [K]。
        :type T: float
        :param m_eff: 有効質量比 m*/m0。
        :type m_eff: float
        :param xi_low: brentqの初期探索範囲の下限。
        :type xi_low: float
        :param xi_high: brentqの初期探索範囲の上限。
        :type xi_high: float
    戻り値:
        :returns: 計算された還元フェルミ準位 xi。
        :rtype: float
    例外:
        :raises ValueError: 電子濃度が正でない場合。
        :raises RuntimeError: 指定された電子濃度に対するxiの探索範囲を見つけられなかった場合。
    """
    if n_m3 <= 0.0:
        raise ValueError(f"Electron density must be positive, got {n_m3}")

    def f(xi):
        return electron_density_from_xi(xi, T, m_eff) - n_m3

    fl = f(xi_low)
    fh = f(xi_high)

    # 探索範囲を動的に拡張
    n_expand = 0
    while fl > 0.0 and n_expand < 20:
        xi_low -= 20.0
        fl = f(xi_low)
        n_expand += 1

    n_expand = 0
    while fh < 0.0 and n_expand < 20:
        xi_high += 20.0
        fh = f(xi_high)
        n_expand += 1

    if fl * fh > 0.0:
        raise RuntimeError(
            "Failed to bracket xi for the requested electron density. "
            f"n={n_m3:.6e} m^-3, T={T}, m_eff={m_eff}, "
            f"f({xi_low})={fl:.6e}, f({xi_high})={fh:.6e}"
        )

    return brentq(f, xi_low, xi_high, xtol=1e-12, rtol=1e-12, maxiter=200)


def seebeck_from_xi(xi, r, carrier="electron"):
    """
    概要:
        還元フェルミ準位 xi からゼーベック係数を計算します。
    詳細説明:
        還元フェルミ準位 xi、散乱因子 r、キャリアタイプに基づいてゼーベック係数 S を計算します。
        ゼーベック係数は以下の式で計算されます。
        S = (k_B / q) * [((r + 2) / (r + 1)) * (Fj(r+1, xi) / Fj(r, xi)) - xi]
        ここで q はキャリアの電荷（電子の場合は -e_SI、正孔の場合は +e_SI）です。
        散乱因子 r は -1.0 より大きい必要があります。
    引数:
        :param xi: 還元フェルミ準位。
        :type xi: float
        :param r: 散乱因子。r は -1 より大きい必要があります。
        :type r: float
        :param carrier: キャリアのタイプ ("electron" または "hole")。
        :type carrier: str
    戻り値:
        :returns: 計算されたゼーベック係数 S [V/K]。
        :rtype: float
    例外:
        :raises ValueError: キャリアタイプが 'electron' または 'hole' ではない場合、または r が -1.0 以下の場合。
    """
    if carrier.lower().startswith("e"):
        q = -e_SI
    elif carrier.lower().startswith("h"):
        q = +e_SI
    else:
        raise ValueError("carrier must be 'electron' or 'hole'")

    if r <= -1.0:
        raise ValueError(f"r must be > -1, got r={r}")

    Fr = Fj(r, xi)
    Fr1 = Fj(r + 1.0, xi)
    bracket = ((r + 2.0) / (r + 1.0)) * (Fr1 / Fr) - xi
    return (kB_SI / q) * bracket


def parse_r_list(text):
    """
    概要:
        カンマ区切りの文字列から散乱因子 r のリストをパースします。
    詳細説明:
        入力文字列を解析し、浮動小数点数の散乱因子 r のリストを返します。
        各散乱因子は -1 より大きい必要があります。
    引数:
        :param text: カンマで区切られた散乱因子の文字列（例: "2,1.5,1"）。
        :type text: str
    戻り値:
        :returns: パースされた散乱因子 r のリスト。
        :rtype: list[float]
    例外:
        :raises ValueError: r のリストが空の場合、またはリスト内のいずれかの r が -1.0 以下の場合。
    """
    vals = [float(x.strip()) for x in text.split(",") if x.strip()]
    if not vals:
        raise ValueError("r-list is empty")
    for r in vals:
        if r <= -1.0:
            raise ValueError(f"Each r must be > -1, got {r}")
    return vals


def default_scan_settings(mode):
    """
    概要:
        指定されたスキャンモードのデフォルト設定を返します。
    詳細説明:
        スキャンモードに基づいて、スキャン範囲の最小値、最大値、点数、および対数スケールフラグのデフォルト値を決定します。
    引数:
        :param mode: スキャンモード ("EF", "Ne", "T", "m_eff")。
        :type mode: str
    戻り値:
        :returns: (scan_min, scan_max, npts, logscan) のタプル。
        :rtype: tuple
    例外:
        :raises ValueError: 未知のスキャンモードが指定された場合。
    """
    if mode == "EF":
        return -0.2, 0.5, 121, False
    if mode == "Ne":
        return 1e15, 1e21, 121, True
    if mode == "T":
        return 100.0, 800.0, 121, False
    if mode == "m_eff":
        return 0.1, 3.0, 121, False
    raise ValueError(f"Unknown mode: {mode}")


def make_scan_values(mode, scan_min, scan_max, npts, logscan):
    """
    概要:
        スキャンする値の配列を生成します。
    詳細説明:
        指定された範囲 (scan_min, scan_max) と点数 npts に基づいて、
        線形または対数スケールのスキャン値配列を生成します。
        対数スケールの場合、scan_min と scan_max は正の値である必要があります。
    引数:
        :param mode: スキャンモード。この引数は主にエラーメッセージの表示に使用されます。
        :type mode: str
        :param scan_min: スキャン範囲の最小値。
        :type scan_min: float
        :param scan_max: スキャン範囲の最大値。
        :type scan_max: float
        :param npts: スキャン点の数。
        :type npts: int
        :param logscan: スキャン値を対数スケールで生成するかどうか。
        :type logscan: bool
    戻り値:
        :returns: 生成されたスキャン値の配列。
        :rtype: numpy.ndarray
    例外:
        :raises ValueError: 対数スキャンで scan_min または scan_max が0以下の場合。
    """
    if logscan:
        if scan_min <= 0.0 or scan_max <= 0.0:
            raise ValueError("log scan requires scan_min > 0 and scan_max > 0")
        return np.logspace(np.log10(scan_min), np.log10(scan_max), npts)
    return np.linspace(scan_min, scan_max, npts)


def resolve_xi(mode, scan_value, T0, m_eff0, EF0_eV, Ne0_cm3, xi0, fixed):
    """
    概要:
        現在のスキャン値とモードに基づいて、還元フェルミ準位 (xi)、フェルミ準位 (EF_eV)、
        電子濃度 (Ne_cm3)、温度 (T_K)、および有効質量 (m_eff) を解決します。
    詳細説明:
        この関数は、mode 引数に応じて、scan_value を基に物理量を計算します。
        - mode が "EF" の場合、scan_value は EF_eV として扱われ、Tとm_effは T0 と m_eff0 を使用します。
        - mode が "Ne" の場合、scan_value は Ne_cm3 として扱われ、Tとm_effは T0 と m_eff0 を使用します。
        - mode が "T" の場合、scan_value は T_K として扱われ、m_effは m_eff0 を使用し、
          fixed 引数（"EF", "Ne", "xi"）に基づいて EF_eV または Ne_cm3 を固定して xi を計算します。
        - mode が "m_eff" の場合、scan_value は m_eff として扱われ、Tは T0 を使用し、
          fixed 引数（"EF", "Ne", "xi"）に基づいて EF_eV または Ne_cm3 を固定して xi を計算します。
    引数:
        :param mode: 現在のスキャンモード ("EF", "Ne", "T", "m_eff")。
        :type mode: str
        :param scan_value: 現在のスキャン対象の値。
        :type scan_value: float
        :param T0: 基準温度 [K]。
        :type T0: float
        :param m_eff0: 基準有効質量比 m*/m0。
        :type m_eff0: float
        :param EF0_eV: 基準フェルミ準位 [eV]。mode が "T" または "m_eff" で fixed が "EF" の場合に固定値として使用。
        :type EF0_eV: float
        :param Ne0_cm3: 基準電子濃度 [cm^-3]。mode が "T" または "m_eff" で fixed が "Ne" の場合に固定値として使用。
        :type Ne0_cm3: float
        :param xi0: 基準還元フェルミ準位。mode が "T" または "m_eff" で fixed が "xi" の場合に固定値として使用。
        :type xi0: float
        :param fixed: mode が "T" または "m_eff" の場合に固定する変数 ("EF", "Ne", "xi")。
        :type fixed: str
    戻り値:
        :returns: (xi, EF_eV, Ne_cm3, T_K, m_eff) のタプル。
        :rtype: tuple[float, float, float, float, float]
    例外:
        :raises ValueError: 未知のスキャンモードまたは固定条件が指定された場合。
    """
    if mode == "EF":
        T = T0
        m_eff = m_eff0
        EF_eV = scan_value
        xi = EF_eV / (kB_eV * T)
    elif mode == "Ne":
        T = T0
        m_eff = m_eff0
        Ne_cm3 = scan_value
        xi = xi_from_electron_density(Ne_cm3 * 1e6, T=T, m_eff=m_eff)
        EF_eV = xi * kB_eV * T
    elif mode == "T":
        T = scan_value
        m_eff = m_eff0
        if fixed == "EF":
            EF_eV = EF0_eV
            xi = EF_eV / (kB_eV * T)
        elif fixed == "Ne":
            xi = xi_from_electron_density(Ne0_cm3 * 1e6, T=T, m_eff=m_eff)
            EF_eV = xi * kB_eV * T
        elif fixed == "xi":
            xi = xi0
            EF_eV = xi * kB_eV * T
        else:
            raise ValueError(f"Unknown fixed condition: {fixed}")
    elif mode == "m_eff":
        T = T0
        m_eff = scan_value
        if fixed == "EF":
            EF_eV = EF0_eV
            xi = EF_eV / (kB_eV * T)
        elif fixed == "Ne":
            xi = xi_from_electron_density(Ne0_cm3 * 1e6, T=T, m_eff=m_eff)
            EF_eV = xi * kB_eV * T
        elif fixed == "xi":
            xi = xi0
            EF_eV = xi * kB_eV * T
        else:
            raise ValueError(f"Unknown fixed condition: {fixed}")
    else:
        raise ValueError(f"Unknown mode: {mode}")

    Ne_cm3 = electron_density_from_xi(xi, T=T, m_eff=m_eff) / 1e6
    return xi, EF_eV, Ne_cm3, T, m_eff


def x_label(mode):
    """
    概要:
        指定されたスキャンモードに対応するX軸ラベル文字列を返します。
    詳細説明:
        プロットのX軸に表示する適切なラベルを返します。
        返されるラベルはLaTeX形式の数式を含む場合があります。
    引数:
        :param mode: スキャンモード ("EF", "Ne", "T", "m_eff")。
        :type mode: str
    戻り値:
        :returns: X軸のラベル文字列。
        :rtype: str
    """
    return {
        "EF": r"$E_F-E_c$ (eV)",
        "Ne": r"$N_e$ (cm$^{-3}$)",
        "T": r"$T$ (K)",
        "m_eff": r"$m^*/m_0$",
    }[mode]


def make_title(mode, carrier, fixed, T0, m_eff0, EF0_eV, Ne0_cm3, xi0):
    """
    概要:
        プロットのタイトル文字列を生成します。
    詳細説明:
        スキャンモード、キャリアタイプ、固定条件などのパラメータに基づいて、グラフのタイトルを作成します。
        タイトルには、スキャンモード、キャリアタイプ、および固定されている物理量とその値が含まれます。
    引数:
        :param mode: スキャンモード ("EF", "Ne", "T", "m_eff")。
        :type mode: str
        :param carrier: キャリアタイプ ("electron" または "hole")。
        :type carrier: str
        :param fixed: 固定された変数 ("EF", "Ne", "xi")。mode が "T" または "m_eff" の場合にのみ関連します。
        :type fixed: str
        :param T0: 基準温度 [K]。
        :type T0: float
        :param m_eff0: 基準有効質量比 m*/m0。
        :type m_eff0: float
        :param EF0_eV: 基準フェルミ準位 [eV]。
        :type EF0_eV: float
        :param Ne0_cm3: 基準電子濃度 [cm^-3]。
        :type Ne0_cm3: float
        :param xi0: 基準還元フェルミ準位。
        :type xi0: float
    戻り値:
        :returns: 生成されたタイトル文字列。
        :rtype: str
    """
    base = f"Seebeck coefficient scan ({mode}, {carrier})"
    if mode in ["EF", "Ne"]:
        return base + f"\nT={T0:g} K, m*/m0={m_eff0:g}"
    if fixed == "EF":
        fixed_text = f"fixed EF={EF0_eV:g} eV"
    elif fixed == "Ne":
        fixed_text = f"fixed Ne={Ne0_cm3:.3e} cm^-3"
    else: # fixed == "xi"
        fixed_text = f"fixed xi={xi0:g}"
    return base + f"\n{fixed_text}, T={T0:g} K, m*/m0={m_eff0:g}"


def save_csv(outfile, rows):
    """
    概要:
        計算結果をCSVファイルに保存します。
    詳細説明:
        計算されたゼーベック係数と関連する物理量を、指定されたファイル名でCSV形式で保存します。
        outfile が空文字列の場合、ファイル保存はスキップされます。
        CSVヘッダーは固定されており、以下のフィールドが含まれます:
        "mode", "scan_value", "r", "S_V_per_K", "S_uV_per_K",
        "xi", "EF_eV", "Ne_cm3", "T_K", "m_eff_m0"
    引数:
        :param outfile: 出力CSVファイルのパス。空文字列の場合、ファイルは保存されません。
        :type outfile: str
        :param rows: 各行が辞書形式で表現されたデータのリスト。
        :type rows: list[dict]
    戻り値:
        :returns: None
        :rtype: None
    """
    if not outfile:
        return
    fieldnames = [
        "mode", "scan_value", "r", "S_V_per_K", "S_uV_per_K",
        "xi", "EF_eV", "Ne_cm3", "T_K", "m_eff_m0"
    ]
    with open(outfile, "w", newline="", encoding="utf-8-sig") as f:
        writer = csv.DictWriter(f, fieldnames=fieldnames)
        writer.writeheader()
        writer.writerows(rows)


def main():
    """
    概要:
        ゼーベック係数スキャンプログラムのメインエントリポイント。
    詳細説明:
        コマンドライン引数をパースし、指定されたモードとパラメータに基づいてゼーベック係数を計算し、結果をプロットまたはCSVファイルに保存します。
        argparse を使用して、スキャンモード、固定条件、キャリアタイプ、基準温度、有効質量、フェルミ準位、電子濃度、還元フェルミ準位、散乱因子のリスト、スキャン範囲、点数、対数スキャンオプション、X軸スケール、CSV出力ファイル名、プロット表示の有無などの引数を定義します。
        デフォルトのスキャン設定を適用し、スキャン値の配列を生成します。
        指定された散乱因子のリスト r_list の各値に対して、スキャン値 (xvals) ごとにゼーベック係数を計算します。
        各計算結果は rows リストに追加され、後でCSVとして保存するために使用されます。
        計算結果は matplotlib を使用してプロットされ、適切なラベル、タイトル、凡例が設定されます。
        --csv オプションが指定された場合、save_csv 関数を呼び出して結果を保存します。
        --no-show オプションが指定されない限り、プロットを表示します。
    引数:
        None
    戻り値:
        None
    """
    parser = argparse.ArgumentParser(
        description="Unified Seebeck scan script for EF / Ne / T / m_eff"
    )
    parser.add_argument("--mode", choices=["EF", "Ne", "T", "m_eff"], default="EF",
                        help="scan target")
    parser.add_argument("--fixed", choices=["EF", "Ne", "xi"], default="EF",
                        help="for mode=T or m_eff, what to keep fixed")
    parser.add_argument("--carrier", choices=["electron", "hole"], default="electron")
    parser.add_argument("--T", type=float, default=300.0,
                        help="reference temperature [K]")
    parser.add_argument("--m_eff", type=float, default=1.0,
                        help="reference effective mass m*/m0")
    parser.add_argument("--EF", type=float, default=0.10,
                        help="reference EF-Ec [eV] (or corresponding band-edge reference)")
    parser.add_argument("--Ne", type=float, default=1e18,
                        help="reference electron density [cm^-3]")
    parser.add_argument("--xi", type=float, default=0.0,
                        help="reference reduced Fermi energy xi")
    parser.add_argument("--r-list", default="2,1.5,1,0.5,0,-0.5",
                        help="comma-separated r list")
    parser.add_argument("--scan-min", type=float, default=None)
    parser.add_argument("--scan-max", type=float, default=None)
    parser.add_argument("--n", type=int, default=None,
                        help="number of scan points")
    parser.add_argument("--logscan", action="store_true",
                        help="use log spacing for scan values")
    parser.add_argument("--xscale", choices=["auto", "linear", "log"], default="auto")
    parser.add_argument("--csv", default="",
                        help="optional CSV output filename")
    parser.add_argument("--no-show", action="store_true",
                        help="do not display plot")
    args = parser.parse_args()

    r_list = parse_r_list(args.r_list)

    dmin, dmax, dn, dlog = default_scan_settings(args.mode)
    scan_min = dmin if args.scan_min is None else args.scan_min
    scan_max = dmax if args.scan_max is None else args.scan_max
    npts = dn if args.n is None else args.n
    logscan = args.logscan or dlog

    if npts < 2:
        raise ValueError("n must be >= 2")
    if scan_max <= scan_min:
        raise ValueError("scan-max must be larger than scan-min")
    if args.T <= 0.0:
        raise ValueError("T must be positive")
    if args.m_eff <= 0.0:
        raise ValueError("m_eff must be positive")
    if args.Ne <= 0.0:
        raise ValueError("Ne must be positive")

    xvals = make_scan_values(args.mode, scan_min, scan_max, npts, logscan)

    rows = []
    plt.figure(figsize=(8, 6))

    for r in r_list:
        yvals_uV = []
        for x in xvals:
            xi, EF_eV, Ne_cm3, T_K, m_eff = resolve_xi(
                mode=args.mode,
                scan_value=float(x),
                T0=args.T,
                m_eff0=args.m_eff,
                EF0_eV=args.EF,
                Ne0_cm3=args.Ne,
                xi0=args.xi,
                fixed=args.fixed,
            )
            S = seebeck_from_xi(xi, r=r, carrier=args.carrier)
            yvals_uV.append(S * 1e6)
            rows.append({
                "mode": args.mode,
                "scan_value": float(x),
                "r": float(r),
                "S_V_per_K": float(S),
                "S_uV_per_K": float(S * 1e6),
                "xi": float(xi),
                "EF_eV": float(EF_eV),
                "Ne_cm3": float(Ne_cm3),
                "T_K": float(T_K),
                "m_eff_m0": float(m_eff),
            })

        plt.plot(xvals, yvals_uV, label=f"r = {r:g}")

    save_csv(args.csv, rows)

    plt.axhline(0.0, linewidth=1.0, alpha=0.4)
    plt.xlabel(x_label(args.mode))
    plt.ylabel(r"$S$ ($\mu$V/K)")
    plt.title(make_title(args.mode, args.carrier, args.fixed, args.T, args.m_eff, args.EF, args.Ne, args.xi))

    xscale = args.xscale
    if xscale == "auto":
        xscale = "log" if (args.mode == "Ne" or logscan) else "linear"
    plt.xscale(xscale)

    plt.grid(True, which="both", alpha=0.25)
    plt.legend()
    plt.tight_layout()

    if not args.no_show:
        plt.show()


if __name__ == "__main__":
    main()