"""
Fermi-Dirac積分を高速かつ堅牢に計算するためのライブラリ。

このモジュールは、フェルミ・ディラック積分 `F_r(eta)` の様々な次数 `r` に対する効率的な計算を提供します。
半導体物理学や統計物理学におけるキャリア密度や関連する熱力学的特性の計算に利用できます。
特に、整数次数 `r=-1, 0, 1` に対しては解析解を用いた高速化が施されており、
その他の次数では漸近展開や数値積分を組み合わせて高い精度と性能を実現しています。

:doc:`tkfermi_integral_usage`
"""

from math import exp, gamma, isclose
from functools import lru_cache
import numpy as np
from scipy import integrate
from scipy.special import expit, log_expit
from scipy.special import zeta
from scipy.special import spence  # 追加：Li2(x)の計算に使用


xlim_exp = 40.0
#xlim_exp = 300.0

kF12 = 2.0 / np.sqrt(np.pi)  # 2/sqrt(pi)


# ------------------------------------------------------------
# Sommerfeld expansion coefficients for degenerate limit
# z_k = (1 - 2^(-(2k+1))) * zeta(2k+2)
# where zeta(s, 1) = Riemann zeta(s)
# for validation: _deg_coeffs = (0.822467, 0.947033, 0.985551, 0.996233, 0.999040, 0.999758, 0.999939, 0.999985)
# ------------------------------------------------------------
n_deg_coeffs = 10

"""
_deg_coeffs = [
    (1.0 - 2.0 ** (-(2 * k + 1))) * float(zeta(2 * k + 2, 1.0))
    for k in range(n_deg_coeffs)
]
"""

def generate_sommerfeld_coeffs(n_terms):
    """
    ゾンマーフェルト展開に使用される係数を生成します。

    この関数は、フェルミ・ディラック積分の縮退近似（ゾンマーフェルト展開）
    において使用される定数 `z_k = (1 - 2^(1-2k)) * zeta(2k)` を計算します。
    `zeta(s, q)` はフルヴィッツのゼータ関数です。

    :param n_terms: 生成する係数の数。
    :type n_terms: int
    :returns: 生成されたゾンマーフェルト展開係数のタプル。
    :rtype: tuple[float]
    """
    coeffs = []
    for k in range(1, n_terms + 1):
        # z_k = (1 - 2^(1-2k)) * zeta(2k)
        val = (1 - 2.0**(1 - 2*k)) * zeta(2*k, 1)
        coeffs.append(float(val))
    return tuple(coeffs)

_deg_coeffs = generate_sommerfeld_coeffs(n_deg_coeffs)

def _round_key(x, ndigit=10):
    """
    数値を指定された小数点以下の桁数で丸めます。

    この関数は、`lru_cache` のキーとして使用される数値を標準化するために設計されています。
    浮動小数点数の比較における微細な違いを吸収し、キャッシュヒット率を向上させます。

    :param x: 丸める対象の数値。
    :type x: float
    :param ndigit: 丸める小数点以下の桁数。デフォルトは10です。
    :type ndigit: int
    :returns: 指定された桁数で丸められた浮動小数点数。
    :rtype: float
    """
    return round(float(x), ndigit)


# ------------------------------------------------------------
# Fermi–Dirac occupations (x = E/kT, eta = Ef/kT)
# ------------------------------------------------------------
def fermi_dirac_e_npvector(x):
    """
    与えられたエネルギー `x = (E - E_f) / kT` のNumpy配列に対するフェルミ・ディラック分布関数を計算します。

    この関数は、数値的な安定性を確保するために、`x` の値を `[-xlim_exp, xlim_exp]` の範囲にクリップします。
    また、`x` の符号に応じて計算方法を切り替えることで、大きな `x` の値に対してもオーバーフローやアンダーフローを防ぎます。

    :param x: 無次元エネルギー `(E - E_f) / kT` を表すNumpy配列。
    :type x: numpy.ndarray
    :returns: 各エネルギー `x` に対応するフェルミ・ディラック分布関数の値。
    :rtype: numpy.ndarray
    """
    x = np.clip(x, -xlim_exp, xlim_exp)
    pos = x > 0
    out = np.empty_like(x)
    exp_positive = np.exp(-x[pos])
    out[pos] = exp_positive / (1 + exp_positive)
    out[~pos] = 1 / (1 + np.exp(x[~pos]))
    return out

def fermi_dirac_e_if(x, eta, g=1.0):
    """
    フェルミ・ディラック分布関数 `f = 1 / (1 + g * exp(x - eta))` を計算します。

    この関数は、`x - eta` の値に基づいて条件分岐を使用し、
    指数関数 `exp()` の引数が非常に大きくなったり小さくなったりする場合の
    数値的な安定性（オーバーフローやアンダーフロー）を確保します。
    `xlim_exp` の範囲外では、結果を `0.0` または `1.0` にクリップします。

    :param x: 無次元エネルギー `E / kT`。
    :type x: float
    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float
    :param g: スピン縮重因子。デフォルトは1.0です。
    :type g: float
    :returns: フェルミ・ディラック分布関数の値。
    :rtype: float
    """
    xe = x - eta
    if xe < -xlim_exp: return 1.0
    if xe > xlim_exp: return 0.0
    return 1.0 / (1 + g * exp(xe))

def fermi_dirac_e(x, eta, g=1.0):
    """
    フェルミ・ディラック分布関数 `f = 1 / (1 + g * exp(x - eta))` を計算します。

    この関数は、`scipy.special.expit` (ロジスティックシグモイド関数) を使用して、
    `exp(x - eta)` が非常に大きいまたは小さい場合に発生する可能性のある
    数値的な不安定性（オーバーフローやアンダーフロー）を回避し、
    よりロバストな計算を提供します。

    :param x: 無次元エネルギー `E / kT`。
    :type x: float or numpy.ndarray
    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float or numpy.ndarray
    :param g: スピン縮重因子。デフォルトは1.0です。
    :type g: float
    :returns: フェルミ・ディラック分布関数の値。`x` または `eta` がNumpy配列の場合、結果もNumpy配列となります。
    :rtype: float or numpy.ndarray
    """
    return expit((eta - x) - np.log(g))

def log_fermi_dirac_e(x, eta, g=1.0):
    """
    フェルミ・ディラック分布関数 `f` の自然対数 `log(f)` を安定的に計算します。

    この関数は `f = 1 / (1 + g * exp(x - eta))` の対数を計算しますが、
    `scipy.special.log_expit` を使用することで、`f` が非常に小さい（0に近い）場合に
    `log(f)` を直接計算する際のアンダーフロー問題を回避します。

    :param x: 無次元エネルギー `E / kT`。
    :type x: float or numpy.ndarray
    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float or numpy.ndarray
    :param g: スピン縮重因子。デフォルトは1.0です。
    :type g: float
    :returns: フェルミ・ディラック分布関数の自然対数値。`x` または `eta` がNumpy配列の場合、結果もNumpy配列となります。
    :rtype: float or numpy.ndarray
    """
    return log_expit((eta - x) - np.log(g))

def fermi_dirac_h(x, eta, g=1.0):
    """
    正孔（ホール）の占有確率 `1 - f` を計算します。

    この関数は、`1 - f = 1 / (1 + g * exp(eta - x))` という形で定義される正孔の確率を、
    `scipy.special.expit` を使用して数値的に安定な方法で計算します。
    フェルミ・ディラック分布関数 `f` と同様に、`g` はスピン縮重因子です。

    :param x: 無次元エネルギー `E / kT`。
    :type x: float or numpy.ndarray
    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float or numpy.ndarray
    :param g: スピン縮重因子。デフォルトは1.0です。
    :type g: float
    :returns: 正孔の占有確率 `1 - f` の値。`x` または `eta` がNumpy配列の場合、結果もNumpy配列となります。
    :rtype: float or numpy.ndarray
    """
    return expit((x - eta) - np.log(g))

def log_fermi_dirac_h(x, eta, g=1.0):
    """
    正孔（ホール）の占有確率 `1 - f` の自然対数 `log(1 - f)` を安定的に計算します。

    この関数は `1 - f = 1 / (1 + g * exp(eta - x))` の対数を計算しますが、
    `scipy.special.log_expit` を使用することで、`1 - f` が非常に小さい（0に近い）場合に
    `log(1 - f)` を直接計算する際のアンダーフロー問題を回避します。

    :param x: 無次元エネルギー `E / kT`。
    :type x: float or numpy.ndarray
    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float or numpy.ndarray
    :param g: スピン縮重因子。デフォルトは1.0です。
    :type g: float
    :returns: 正孔の占有確率の自然対数値。`x` または `eta` がNumpy配列の場合、結果もNumpy配列となります。
    :rtype: float or numpy.ndarray
    """
    return log_expit((x - eta) - np.log(g))

def ionized_acceptor_frac(eta, xA, g=1.0):
    """
    イオン化されたアクセプターの割合を計算します。

    この関数は、アクセプター準位が電子で占有されている確率として定義される、
    イオン化されたアクセプターの割合 `N_A^- / N_A` を計算します。
    基本的なフェルミ・ディラック分布関数 `fermi_dirac_e` のラッパーです。

    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float or numpy.ndarray
    :param xA: 無次元アクセプター準位エネルギー `E_A / kT`。
    :type xA: float or numpy.ndarray
    :param g: スピン縮重因子。デフォルトは1.0です。
    :type g: float
    :returns: イオン化されたアクセプターの割合。`eta` または `xA` がNumpy配列の場合、結果もNumpy配列となります。
    :rtype: float or numpy.ndarray
    """
    return fermi_dirac_e(xA, eta, g=g)

def ionized_donor_frac(eta, xD, g=1.0):
    """
    イオン化されたドナーの割合を計算します。

    この関数は、ドナー準位が正孔で占有されている（つまり電子を放出した）確率として定義される、
    イオン化されたドナーの割合 `N_D^+ / N_D` を計算します。
    基本的な正孔の占有確率関数 `fermi_dirac_h` のラッパーです。

    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float or numpy.ndarray
    :param xD: 無次元ドナー準位エネルギー `E_D / kT`。
    :type xD: float or numpy.ndarray
    :param g: スピン縮重因子。デフォルトは1.0です。
    :type g: float
    :returns: イオン化されたドナーの割合。`eta` または `xD` がNumpy配列の場合、結果もNumpy配列となります。
    :rtype: float or numpy.ndarray
    """
    return fermi_dirac_h(xD, eta, g=g)

# ============================================================
# internal: degenerate limit (Sommerfeld-like)
# ============================================================

def _DegenerateFermiDirac(eta, r):
    """
    縮退領域におけるフェルミ・ディラック積分 `F_r(eta)` を計算します。

    この関数は、化学ポテンシャル `eta` が十分に大きい場合（縮退領域）に適用される
    ゾンマーフェルト展開に基づいています。`_deg_coeffs` は、
    展開の係数 `z_k = (1 - 2^(1-2k)) * zeta(2k)` を事前に計算したものです。
    これにより、数値積分の代わりに解析的な近似を用いることで高速な計算が可能です。

    :param eta: 無次元化学ポテンシャル `E_f / kT`。`eta` は正で十分に大きいと仮定されます。
    :type eta: float
    :param r: フェルミ・ディラック積分の次数。
    :type r: float
    :returns: 縮退領域におけるフェルミ・ディラック積分の近似値。
    :rtype: float
    """
    gi = eta ** (r + 1.0) / (r + 1.0)
    t = r
    gg = t / (eta * eta)
    gs = gg * _deg_coeffs[0]
    t -= 1.0

    for k in range(1, n_deg_coeffs):
        gg *= t * (t - 1.0) / (eta * eta)
        gs += gg * _deg_coeffs[k]
        t -= 2.0

    return gi * (1.0 + 2.0 * (r + 1.0) * gs)


# ============================================================
# internal: non-degenerate limit (alternating series)
# ============================================================

def _NonDegenerateFermiDirac(eta, r, reltol=1e-10, max_terms=200000):
    """
    非縮退領域におけるフェルミ・ディラック積分 `F_r(eta)` を計算します。

    この関数は、化学ポテンシャル `eta` が小さい場合（非縮退領域）に適用されます。
    `eta` が非常に小さい場合は `Gamma(r+1) * exp(eta)` の近似を使用し、
    それ以外の非縮退領域では、交代級数展開を用いて積分を計算します。
    級数の収束は `reltol` で制御されます。

    :param eta: 無次元化学ポテンシャル `E_f / kT`。`eta` は負で小さいと仮定されます。
    :type eta: float
    :param r: フェルミ・ディラック積分の次数。
    :type r: float
    :param reltol: 級数展開の収束を判定するための相対許容誤差。デフォルトは1e-10です。
    :type reltol: float
    :param max_terms: 級数展開の最大項数。デフォルトは200000です。
    :type max_terms: int
    :returns: 非縮退領域におけるフェルミ・ディラック積分の近似値。
    :rtype: float
    """
    if eta < -20.0:
        return gamma(r + 1.0) * exp(eta)

    s = 1.0
    total = 0.0
    scale = 0.0

    for _ in range(max_terms // 2):
        t1 = exp(s * eta) / s**(r + 1.0)
        s += 1.0
        t2 = exp(s * eta) / s**(r + 1.0)
        s += 1.0

        subtotal = t1 - t2
        total += subtotal

        scale = max(scale, abs(t1), abs(t2))

        if abs(subtotal) < reltol * max(abs(total), 1.0):
            break

    return total * gamma(r+1.0)


# ============================================================
# cache-safe wrapper
# ============================================================

def _func_Fr_x2(x, eta, r):
    """
    フェルミ・ディラック積分 `F_r(eta)` の数値積分で使用される補助関数。

    この関数は、フェルミ・ディラック積分の定義
    `F_r(eta) = ∫_0^∞ ε^r / (1 + exp(ε - eta)) dε`
    において、変数変換 `ε = x^2` を行った後の被積分関数の一部を計算します。
    具体的には、`2.0 * x^(2.0 * r + 1.0) * f(x^2, eta)` に相当します。
    `expit` 関数を用いて数値安定性を確保しています。

    :param x: 積分変数 `x` (元のエネルギー `ε` の平方根に相当)。
    :type x: float
    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float
    :param r: フェルミ・ディラック積分の次数。
    :type r: float
    :returns: 被積分関数の値。
    :rtype: float
    """
    return 2.0 * x ** (2.0 * r + 1.0) * expit(eta - x*x)

def _FermiIntegral_core(eta, r, epsabs, epsrel, limit):
    """
    フェルミ・ディラック積分 `F_r(eta)` の中心的な計算ロジック。

    この関数は、`eta` の値に基づいて、以下の3つの主要な計算方法を切り替えます。

    1.  **縮退領域 (`eta >= 25.0 + 2.0 * r`)**:
        `_DegenerateFermiDirac` を使用したゾンマーフェルト展開による近似。
    2.  **非縮退領域 (`eta < -0.1`)**:
        `_NonDegenerateFermiDirac` を使用した交代級数展開による近似。
    3.  **中間領域**:
        `scipy.integrate.quad` を使用した数値積分。
        数値安定性と効率のために、積分範囲 `[0, 1]`, `[1, 4]`, `[4, lim]` の
        3つの区間に分割して積分を実行します。`lim` は動的に決定されます。

    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float
    :param r: フェルミ・ディラック積分の次数。
    :type r: float
    :param epsabs: 数値積分の絶対許容誤差。
    :type epsabs: float
    :param epsrel: 数値積分の相対許容誤差。
    :type epsrel: float
    :param limit: `scipy.integrate.quad` で使用される最大部分区間数。
    :type limit: int
    :returns: フェルミ・ディラック積分の値。
    :rtype: float
    """
    # degenerate
    if eta >= 25.0 + 2.0 * r:
        return _DegenerateFermiDirac(eta, r)

    # non-degenerate
    if eta < -0.1:
        return _NonDegenerateFermiDirac(eta, r)

    # intermediate: z = x^2 + quad
    lim = 17.0 + 4.5 * r + eta
    lim = max(lim, 4.0)
    f = lambda x: _func_Fr_x2(x, eta, r)

    ret = 0.0
    ret += integrate.quad(f, 0.0, 1.0, epsabs=epsabs, epsrel=epsrel, limit=limit)[0]
    ret += integrate.quad(f, 1.0, 4.0, epsabs=epsabs, epsrel=epsrel, limit=limit)[0]
    ret += integrate.quad(f, 4.0, lim, epsabs=epsabs, epsrel=epsrel, limit=limit)[0]

    return ret

@lru_cache(maxsize=1024)
def _FermiIntegral_cached(eta, r, epsabs, epsrel, limit):
    """
    フェルミ・ディラック積分の計算結果をキャッシュするラッパー関数。

    この関数は `_FermiIntegral_core` を呼び出して実際の計算を実行しますが、
    `functools.lru_cache` デコレータにより、同じ引数セットで複数回呼び出された場合、
    以前の計算結果を返します。これにより、計算コストの高い積分処理のパフォーマンスが向上します。
    キャッシュは最大1024エントリを保持します。

    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float
    :param r: フェルミ・ディラック積分の次数。
    :type r: float
    :param epsabs: 数値積分の絶対許容誤差。
    :type epsabs: float
    :param epsrel: 数値積分の相対許容誤差。
    :type epsrel: float
    :param limit: `scipy.integrate.quad` で使用される最大部分区間数。デフォルトは50です。
    :type limit: int
    :returns: フェルミ・ディラック積分の値。
    :rtype: float
    """
    return _FermiIntegral_core(eta, r, epsabs, epsrel, limit)

# ============================================================
# public API
# ============================================================

def FermiIntegral_fast(
    eta, r,
    *,
    epsabs=1.0e-10,
    epsrel=1.0e-8,
    limit=50,
    use_cache = True
):
    """
    一般的なフェルミ・ディラック積分 `F_r(eta)` を高速かつ高精度に計算します。

    この関数は、以下の優先順位で計算方法を適用します。

    1.  **整数次数 `r = -1, 0, 1` の場合**:
        それぞれ `expit(eta)`, `log(1 + exp(eta))`, および `dilogarithm`
        （`scipy.special.spence`）を用いた解析解または特殊関数によって高速に計算されます。
        `r = 1` の非常に大きな `eta` に対しては、漸近展開も利用されます。
    2.  **その他の次数および `eta` の範囲**:
        `_FermiIntegral_core` を呼び出して、縮退近似、非縮退近似、または数値積分を適用します。
        `use_cache` が `True` の場合、計算結果は `_FermiIntegral_cached` を通じてキャッシュされます。

    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float
    :param r: フェルミ・ディラック積分の次数。
    :type r: float
    :param epsabs: (キーワード引数) 数値積分の絶対許容誤差。デフォルトは1.0e-10です。
    :type epsabs: float
    :param epsrel: (キーワード引数) 数値積分の相対許容誤差。デフォルトは1.0e-8です。
    :type epsrel: float
    :param limit: (キーワード引数) `scipy.integrate.quad` で使用される最大部分区間数。デフォルトは50です。
    :type limit: int
    :param use_cache: (キーワード引数) 計算結果をキャッシュするかどうか。デフォルトはTrueです。
    :type use_cache: bool
    :returns: フェルミ・ディラック積分の値。
    :rtype: float
    """
    # --- 整数次数の解析解 / 特殊関数による高速化 ---
    is_integer_r = isclose(r, round(r), abs_tol=1e-12)
    ri = int(round(r)) if is_integer_r else None
    
    if is_integer_r:
    # r = -1: F_{-1}(eta) = 1 / (1 + exp(-eta))
        if ri == -1:
            return expit(eta)

    # r = 0: F_0(eta) = ln(1 + exp(eta))
        if ri == 0:
            if eta > 0:
                return eta + np.log1p(np.exp(-eta))
            else:
                return np.log1p(np.exp(eta))

    # r = 1: F_1(eta) = -Li2(-exp(eta))
    # scipy.special.spence(z) は Li2(z) そのものではなく Li2(z)に関連した定義
    # Li2(z) = spence(1-z) なので、F_1(eta) = -spence(1 + exp(eta))
        if ri == 1:
            if eta > 100:
            # 漸近展開 (eta^2 / 2 + pi^2 / 6)
                return 0.5 * eta**2 + np.pi**2 / 6.0
            return -float(spence(1.0 + np.exp(eta)))

    # キャッシュキーの作成と呼び出し
    key = (
        _round_key(eta),
        _round_key(r),
        epsabs,
        epsrel,
        limit,
    )
    if use_cache:
        return _FermiIntegral_cached(*key)
    else:
        return _FermiIntegral_core(*key)

def FermiIntegral_half(
    eta,
    *,
    epsabs=1.0e-10,
    epsrel=1.0e-8,
    limit=50
):
    """
    フェルミ・ディラック積分 `F_{1/2}(eta)` を計算します。

    この関数は、`FermiIntegral_fast` 関数のラッパーであり、
    積分の次数 `r` を `0.5` に固定して呼び出すことで、特定の次数を便利に計算できるようにします。
    半導体中の電子濃度計算によく用いられます。

    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float
    :param epsabs: (キーワード引数) 数値積分の絶対許容誤差。デフォルトは1.0e-10です。
    :type epsabs: float
    :param epsrel: (キーワード引数) 数値積分の相対許容誤差。デフォルトは1.0e-8です。
    :type epsrel: float
    :param limit: (キーワード引数) `scipy.integrate.quad` で使用される最大部分区間数。デフォルトは50です。
    :type limit: int
    :returns: フェルミ・ディラック積分 `F_{1/2}(eta)` の値。
    :rtype: float
    """

    return FermiIntegral_fast(
        eta, 0.5,
        epsabs=epsabs,
        epsrel=epsrel,
        limit=limit
    )


def FermiIntegral_3half(
    eta,
    *,
    epsabs=1.0e-10,
    epsrel=1.0e-8,
    limit=50
):
    """
    フェルミ・ディラック積分 `F_{3/2}(eta)` を計算します。

    この関数は、`FermiIntegral_fast` 関数のラッパーであり、
    積分の次数 `r` を `1.5` に固定して呼び出すことで、特定の次数を便利に計算できるようにします。
    非放物型バンドにおけるキャリア密度や、その他の物理量の計算に利用されます。

    :param eta: 無次元化学ポテンシャル `E_f / kT`。
    :type eta: float
    :param epsabs: (キーワード引数) 数値積分の絶対許容誤差。デフォルトは1.0e-10です。
    :type epsabs: float
    :param epsrel: (キーワード引数) 数値積分の相対許容誤差。デフォルトは1.0e-8です。
    :type epsrel: float
    :param limit: (キーワード引数) `scipy.integrate.quad` で使用される最大部分区間数。デフォルトは50です。
    :type limit: int
    :returns: フェルミ・ディラック積分 `F_{3/2}(eta)` の値。
    :rtype: float
    """

    return FermiIntegral_fast(
        eta, 1.5,
        epsabs=epsabs,
        epsrel=epsrel,
        limit=limit
    )


# ============================================================
# carrier density helpers (Gamma-normalization NOT used)
# ============================================================

def electron_density(Nc, eta_c, *, epsabs=1.0e-10, epsrel=1.0e-8, limit=50):
    """
    導電帯における電子濃度を計算します（非放物型効果は無視し、標準的な放物型バンドを仮定）。

    このライブラリで定義される（非正規化）フェルミ・ディラック積分 `I_{1/2}(eta)` は、
    `∫_0^∞ ε^{1/2} / (1 + exp(ε - eta)) dε` であり、ガンマ関数による除算は行われません。
    したがって、電子濃度 `n` は以下の式で与えられます。

    `n = Nc * (2/sqrt(pi)) * I_{1/2}(eta_c)`

    この関数は `eta_c` がスカラーでもNumpy配列でも対応し、それぞれの `eta_c` に対して
    `FermiIntegral_half` を呼び出して電子濃度を計算します。

    :param Nc: 導電帯の実効状態密度 [cm^-3] または [m^-3]。単位は結果に引き継がれます。
    :type Nc: float
    :param eta_c: 電子に対する無次元化学ポテンシャル `eta_c = (Ef - Ec) / (kT)`。
    :type eta_c: float or numpy.ndarray
    :param epsabs: (キーワード引数) 数値積分の絶対許容誤差。デフォルトは1.0e-10です。
    :type epsabs: float
    :param epsrel: (キーワード引数) 数値積分の相対許容誤差。デフォルトは1.0e-8です。
    :type epsrel: float
    :param limit: (キーワード引数) `scipy.integrate.quad` で使用される最大部分区間数。デフォルトは50です。
    :type limit: int
    :returns: `Nc` と同じ単位の電子濃度。`eta_c` がNumpy配列の場合、結果もNumpy配列となります。
    :rtype: float or numpy.ndarray
    """

    eta_arr = np.asarray(eta_c)
    if eta_arr.ndim == 0:
        return Nc * kF12 * FermiIntegral_half(
            float(eta_arr), epsabs=epsabs, epsrel=epsrel, limit=limit
        )
    out = np.empty_like(eta_arr, dtype=float)
    it = np.nditer(eta_arr, flags=["multi_index"])
    for v in it:
        out[it.multi_index] = Nc * kF12 * FermiIntegral_half(
            float(v), epsabs=epsabs, epsrel=epsrel, limit=limit
        )
    return out


def hole_density(Nv, eta_v, *, epsabs=1.0e-10, epsrel=1.0e-8, limit=50):
    """
    価電子帯における正孔濃度を計算します（標準的な放物型バンドを仮定）。

    このライブラリで定義される（非正規化）フェルミ・ディラック積分 `I_{1/2}` を使用して、
    正孔濃度 `p` は以下の式で計算されます。

    `p = Nv * (2/sqrt(pi)) * I_{1/2}( beta(Ev - Ef) )`

    `eta_v = (Ef - Ev) / (kT)` として渡される場合、`beta(Ev - Ef) = -eta_v` となるため、
    `FermiIntegral_half(-eta_v)` を呼び出す必要があります。
    この関数は `eta_v` がスカラーでもNumpy配列でも対応し、それぞれの `eta_v` に対して
    `FermiIntegral_half` を呼び出して正孔濃度を計算します。

    :param Nv: 価電子帯の実効状態密度 [cm^-3] または [m^-3]。単位は結果に引き継がれます。
    :type Nv: float
    :param eta_v: 無次元化学ポテンシャル `eta_v = (Ef - Ev) / (kT)`。
    :type eta_v: float or numpy.ndarray
    :param epsabs: (キーワード引数) 数値積分の絶対許容誤差。デフォルトは1.0e-10です。
    :type epsabs: float
    :param epsrel: (キーワード引数) 数値積分の相対許容誤差。デフォルトは1.0e-8です。
    :type epsrel: float
    :param limit: (キーワード引数) `scipy.integrate.quad` で使用される最大部分区間数。デフォルトは50です。
    :type limit: int
    :returns: `Nv` と同じ単位の正孔濃度。`eta_v` がNumpy配列の場合、結果もNumpy配列となります。
    :rtype: float or numpy.ndarray
    """

    eta_arr = np.asarray(eta_v)
    if eta_arr.ndim == 0:
        return Nv * kF12 * FermiIntegral_half(
            float(-eta_arr), epsabs=epsabs, epsrel=epsrel, limit=limit
        )
    out = np.empty_like(eta_arr, dtype=float)
    it = np.nditer(eta_arr, flags=["multi_index"])
    for v in it:
        out[it.multi_index] = Nv * kF12 * FermiIntegral_half(
            float(-v), epsabs=epsabs, epsrel=epsrel, limit=limit
        )
    return out


# ============================================================
# self-test
# ============================================================

if __name__ == "__main__":
    Nc = 1.0e19
    Nv = 1.0e19
    eta_list = [-10.0, -3.0, -1.0, 0.0, 2.0, 5.0, 10.0]

    print("Fermi–Dirac integral test")
    print("")

    print("r = 1/2")
    for eta in eta_list:
        print(f"  eta={eta:6.2f}  F1/2={FermiIntegral_half(eta):.8e}")

    print("\nr = 3/2")
    for eta in eta_list:
        n = electron_density(Nc, eta)
        p = hole_density(Nv, eta)
        print(f"  eta={eta:6.2f}  F3/2={FermiIntegral_3half(eta):.8e}  n={n:.8e}  p={p:.8e}")

    print("\ncache info:")
    print(_FermiIntegral_cached.cache_info())