#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
PV特性評価ツールスクリプト。

概要:
    太陽電池の電気的・光学的特性を評価するツールです。
詳細説明:
    mode=alpha または mode=make_alpha では、反射率と透過率のスペクトルから
    吸収スペクトルを計算し、プロットおよびExcelへの保存を行います。
    mode=analyze では、I-Vデータからパラメータ推定、発電特性解析、量子効率計算、
    および結果のプロットを行います。
主な機能:
    - 反射率/透過率から吸収スペクトルの計算と保存
    - I-Vデータからの太陽電池パラメータ（I0, ndiode, IPV, Rs, Rsh）推定
    - 発電特性（Voc, Jsc, FF, Pmaxなど）の解析
    - 量子効率（EQE, IQE）の計算
    - 解析結果のプロットとExcelファイルへの保存
関連リンク:
    pvanalyze_usage
"""

import sys
import argparse
import builtins
from pathlib import Path
import traceback
import csv
import math

import matplotlib.pyplot as plt
import numpy as np
from openpyxl import Workbook, load_workbook
from openpyxl.styles import Font, PatternFill


KB = 1.380649e-23
E_CHARGE = 1.602176634e-19
H = 6.62607015e-34
C = 2.99792458e8

PARAM_NAMES = ["I0", "ndiode", "IPV", "Rs", "Rsh"]
EPS_I = 1.0e-15

DPV_NM = 46.837
AREA_MM2 = 0.5 * 0.5 * math.pi
AREA_CM2 = AREA_MM2 * 0.01  # 1 mm^2 = 0.01 cm^2
PHOTON_FLUX = 1.95804e18  # cm^-2 s^-1
PHOTON_NM = 1363.0        # nm
PIN = ""                  # W/cm^2, if None and F0 is not None then compute from hc/lambda

fontsize = 16

DARK_IV_FILE = "I_V SweepTest SMU1 [(1) ; 2026_03_24 14_15_24]-SY251213-1-LD0.0V.csv"
PV_IV_FILE = "I_V SweepTest SMU1 [(10) ; 2026_03_24 14_18_37]-SY251213-1-LD1.8V.csv"
R_FILE = "SY251024-1-Bi4O6S21-STO001-225oC-10mJ_R(UDS).txt"
T_FILE = "SY251024-1-Bi4O6S21-STO001-225oC-10mJ_T(UDS).txt"
ALPHA_FILE = "alpha_Bi2OS2.xlsx"


def initialize():
    """
    概要:
        コマンドライン引数を解析し、プログラムの初期設定を行います。
    詳細説明:
        argparseモジュールを使用して、実行モード、ファイルパス、温度、膜厚などの
        パラメータを定義し、ユーザー入力から値をパースします。
    引数:
        なし
    戻り値:
        :returns: 解析された引数オブジェクトと引数パーサーオブジェクトのタプル。
        :rtype: tuple[argparse.Namespace, argparse.ArgumentParser]
    """
    parser = argparse.ArgumentParser(description="PV characterization tool with alpha/analyze modes.")
    parser.add_argument("--mode", default="analyze", choices=["alpha", "make_alpha", "analyze"],
                        help="Execution mode")
    parser.add_argument("--dark", default=DARK_IV_FILE, help="Dark IV CSV file")
    parser.add_argument("--light", default=PV_IV_FILE, help="Illuminated IV CSV file")
    parser.add_argument("--R", default=R_FILE, help="Reflectance spectrum file")
    # NOTE: --T is already requested for temperature, so transmittance file uses --Tr.
    parser.add_argument("--Tr", default=T_FILE, help="Transmittance spectrum file")
    parser.add_argument("--eq", default=None, help="Equivalent spectral irradiance file")
    parser.add_argument("--alpha", default=ALPHA_FILE, help="Excel file for absorption spectrum / alpha")
    parser.add_argument("--T", type=float, default=300.0, help="Temperature (K)")
    parser.add_argument("--d", type=float, default=DPV_NM, help="PV layer thickness (nm)")
    parser.add_argument("--sweep_dark", type=int, default=0, help="Sweep index for dark analysis")
    parser.add_argument("--sweep_light", type=int, default=0, help="Sweep index for light analysis")
    parser.add_argument("--F0", type=float, default=PHOTON_FLUX, help="Incident photon flux (cm^-2 s^-1)")
    parser.add_argument("--P0", type=str, default=PIN, help="Incident photon energy density (W/cm^2)")
    parser.add_argument("--E", type=float, default=0.0, help="Photon energy (eV). If 0, use lambda.")
    parser.add_argument("--lambda_nm", type=float, default=PHOTON_NM,
                        help="Photon wavelength (nm), used when E=0")
    parser.add_argument("--S", type=float, default=AREA_CM2, help="Electrode area (cm^2)")
    parser.add_argument("--outprefix", default="pvanalyze", help="Output file prefix")
    args = parser.parse_args()
    return args, parser


def print_args_and_derived(args, E_use=None, P0_use=None):
    """
    概要:
        解析されたコマンドライン引数と派生値を標準出力に表示します。
    引数:
        :param args: 解析された引数オブジェクト。
        :type args: argparse.Namespace
        :param E_use: 使用されるフォトンエネルギー (eV)。
        :type E_use: float or None
        :param P0_use: 使用される入射光パワー密度 (W/cm^2)。
        :type P0_use: float or None
    戻り値:
        なし
    例外:
        :raises ValueError: フォトン波長 lambda_nm が0以下の場合。
        :raises ValueError: フォトンエネルギー E が0未満の場合。
        :raises ValueError: 入射フォトンフラックス F0 が0以下の場合。
        :raises ValueError: 入射光パワー密度 P0 が0以下の場合。
        :raises ValueError: 温度 T が0以下の場合。
        :raises ValueError: 電極面積 S が0以下の場合。
        :raises ValueError: PV層の膜厚 d が0.1nm以下の場合。
        :raises ValueError: compute_p0 関数から None が返された場合。
    """
    print("=== Arguments ===")
    print(f"mode         = {args.mode}")
    if args.mode == "make_alpha":
        print(f"R            = {args.R}")
        print(f"Tr           = {args.Tr}")
    else:
        print(f"alpha        = {args.alpha}")
    print(f"d            = {args.d} nm")
    print(f"outprefix    = {args.outprefix}")

    if "alpha" not in args.mode:
        print(f"dark         = {args.dark}")
        print(f"light        = {args.light}")
        print(f"T            = {args.T} K")
        print(f"sweep_dark   = {args.sweep_dark}")
        print(f"sweep_light  = {args.sweep_light}")
        if args.F0 is not None:
            print(f"F0           = {args.F0} cm^-2 s^-1")
        if args.P0 != "":
            print(f"P0           = {args.P0} W/cm^2")
#        print(f"eq           = {args.eq}")
        print(f"E            = {args.E} eV")
        print(f"lambda_nm    = {args.lambda_nm} nm")
        print(f"S            = {args.S} cm^2")
    print()

    if args.mode == 'analyze':
        if args.P0 == "":
            if args.E is not None and args.E == 0.0:
                if args.lambda_nm is None or args.lambda_nm <= 0.0:
                    raise ValueError(f"lambda_nm must be > 0 nm, got {args.lambda_nm}")
            elif args.E is not None and args.E < 0.0:
                raise ValueError(f"E must be >= 0 eV, got {args.E}")

            if args.F0 <= 0.0:
                raise ValueError(f"F0 must be > 0 cm^-2 s^-1, got {args.F0}")

            args.E = compute_photon_energy_eV(args.E, args.lambda_nm)
            print(f"E_use  = {args.E} eV")
            args.P0 = compute_p0(args.F0, args.P0, args.E, args.lambda_nm)
            if args.P0 is None:
                raise ValueError(f"Got None from compute_p0() using F0={args.F0} and E={args.E}")
            print(f"P0_use = {args.P0} W/cm^2" if args.P0 is not None else "P0_use       = None")
        else:
            args.P0 = float(args.P0)
            if args.P0 <= 0.0:
                raise ValueError(f"P0 must be > 0 W/cm^2, got {args.P0}")

        if args.T <= 0.0:
            raise ValueError(f"T must be > 0 K, got {args.T}")
        if args.S <= 0.0:
            raise ValueError(f"S must be > 0 cm^2, got {args.S}")

        if args.d <= 0.1:
            raise ValueError(f"d must be > 0.1 nm, got {args.d}")


def read_data(infile, xmin=None, xmax=None, ndataskip=0):
    """
    概要:
        指定されたCSVファイルからI-Vデータを読み込みます。
    詳細説明:
        ファイル内のメタデータ（記録時間、データ名）を抽出し、電圧と電流のデータポイントをパースします。
        複数のスイープを検出し、それぞれをリストに分割します。
    引数:
        :param infile: 読み込むCSVファイルのパス。
        :type infile: str
        :param xmin: X軸（電圧）の最小値。これより小さい値はスキップされます。
        :type xmin: float or None
        :param xmax: X軸（電圧）の最大値。これより大きい値はスキップされます。
        :type xmax: float or None
        :param ndataskip: データポイントをスキップする間隔。0の場合、スキップしません。
        :type ndataskip: int
    戻り値:
        :returns: 各スイープのX軸（電圧）データのリスト、各スイープのY軸（電流）データのリスト、ファイルから抽出されたメタデータのタプル。
        :rtype: tuple[list[numpy.ndarray], list[numpy.ndarray], dict]
    例外:
        :raises ValueError: ファイルから有効なI-Vデータが見つからない場合。
    """
    print(f"[I/O] Read IV data: {infile}")
    raw_x = []
    raw_y = []
    inf = {"FileName": infile, "RecordTime": "Unknown", "DataName": ("V", "I")}

    icount = 0
    data_started = False
    with open(infile, "r", encoding="utf-8-sig", newline="") as f:
        reader = csv.reader(f)
        for row in reader:
            if not row:
                continue
            cols = [c.strip() for c in row]

            if len(cols) >= 3 and cols[0] == "MetaData" and "RecordTime" in cols[1]:
                inf["RecordTime"] = cols[2]

            if cols[0] == "DataName":
                if len(cols) >= 3:
                    inf["DataName"] = (cols[1], cols[2])
                data_started = True
                continue

            if data_started and cols[0] == "DataValue":
                if len(cols) < 3:
                    continue
                icount += 1
                if ndataskip > 0 and (icount - 1) % ndataskip != 0:
                    continue
                try:
                    x = float(cols[1])
                    y = float(cols[2])
                except ValueError:
                    continue
                if xmin is not None and x < xmin:
                    continue
                if xmax is not None and x > xmax:
                    continue
                raw_x.append(x)
                raw_y.append(y)

    if not raw_x:
        raise ValueError(f"No valid IV data found in {infile}")

    xs_list = []
    ys_list = []

    cur_x = [raw_x[0]]
    cur_y = [raw_y[0]]
    prev_sign = 0

    for i in range(1, len(raw_x)):
        dx = raw_x[i] - raw_x[i - 1]
        sign = 0 if abs(dx) < 1e-15 else (1 if dx > 0 else -1)

        if prev_sign == 0:
            prev_sign = sign

        if sign != 0 and prev_sign != 0 and sign != prev_sign:
            xs_list.append(np.asarray(cur_x, dtype=float))
            ys_list.append(np.asarray(cur_y, dtype=float))
            cur_x = [raw_x[i - 1], raw_x[i]]
            cur_y = [raw_y[i - 1], raw_y[i]]
            prev_sign = sign
        else:
            cur_x.append(raw_x[i])
            cur_y.append(raw_y[i])
            if sign != 0:
                prev_sign = sign

    if cur_x:
        xs_list.append(np.asarray(cur_x, dtype=float))
        ys_list.append(np.asarray(cur_y, dtype=float))

    inf["Points"] = sum(len(x) for x in xs_list)
    inf["NSweeps"] = len(xs_list)
    return xs_list, ys_list, inf


def read_alpha_from_excel(infile):
    """
    概要:
        Excelファイルから吸収スペクトル（alpha）データを読み込みます。
    詳細説明:
        save_alpha_to_excel関数によって保存された形式のExcelファイルを想定しています。
        alpha_spectrumシートから波長、R、Tr、A、alphaのデータを抽出します。
    引数:
        :param infile: 読み込むExcelファイルのパス。
        :type infile: str
    戻り値:
        :returns: 波長、反射率、透過率、吸収率、吸収係数を含む辞書。
        :rtype: dict
    例外:
        :raises ValueError: 指定されたExcelファイルに'alpha_spectrum'シートが見つからない場合。
        :raises ValueError: ファイルから数値の吸収スペクトルデータが見つからない場合。
    """
    print(f"[I/O] Read alpha Excel: {infile}")
    wb = load_workbook(infile, data_only=True)
    if "alpha_spectrum" not in wb.sheetnames:
        raise ValueError(f"Sheet 'alpha_spectrum' not found in {infile}")
    ws = wb["alpha_spectrum"]

    wl = []
    R = []
    Tr = []
    A = []
    alpha = []

    first = True
    for row in ws.iter_rows(values_only=True):
        if first:
            first = False
            continue
        if row is None or len(row) < 6:
            continue
        try:
            # columns: photon_energy_eV, wavelength_nm, R, Tr, A, alpha_cm^-1
            wl.append(float(row[1]))
            R.append(float(row[2]))
            Tr.append(float(row[3]))
            A.append(float(row[4]))
            alpha.append(float(row[5]))
        except (TypeError, ValueError):
            continue

    if not wl:
        raise ValueError(f"No numeric alpha data found in {infile}")

    return {
        "wl_nm": np.asarray(wl, dtype=float),
        "R": np.asarray(R, dtype=float),
        "Tr": np.asarray(Tr, dtype=float),
        "A": np.asarray(A, dtype=float),
        "alpha_cm^-1": np.asarray(alpha, dtype=float),
    }


def read_optical_spectrum(infile):
    """
    概要:
        指定されたテキストファイルから光学スペクトルデータ（反射率Rまたは透過率T）を読み込みます。
    詳細説明:
        ファイル内のタブ区切りまたはスペース区切りのデータをパースし、
        波長とスペクトル値のペアを抽出します。
        重複する波長を処理し、波長順にソートします。
    引数:
        :param infile: 読み込む光学スペクトルファイルのパス。
        :type infile: str
    戻り値:
        :returns: 波長 (nm)のnumpy.ndarray、スペクトル値 (RまたはTのパーセンテージ)のnumpy.ndarray、
                  ファイルから抽出されたメタデータを含む辞書のタプル。
        :rtype: tuple[numpy.ndarray, numpy.ndarray, dict]
    例外:
        :raises RuntimeError: ファイルの読み込みに失敗した場合。
        :raises ValueError: ファイルから数値のスペクトルデータが見つからない場合。
    """
    
    print(f"[I/O] Read optical spectrum: {infile}")
    info = {"FileName": infile, "Label": Path(infile).stem, "YUnit": "a.u."}
    lines = []

    encodings = ["cp932", "utf-8-sig", "utf-8", "latin1"]
    last_err = None
    for enc in encodings:
        try:
            with open(infile, "r", encoding=enc, errors="replace") as f:
                lines = f.readlines()
            break
        except Exception as e:
            last_err = e
    if not lines:
        raise RuntimeError(f"Failed to read optical spectrum: {infile} ({last_err})")

    wl = []
    val = []

    for line in lines:
        s = line.strip()
        if not s:
            continue

        if "\t" in s:
            parts = [p.strip() for p in s.split("\t")]
            if len(parts) >= 2:
                if "ｻﾝﾌﾟﾙ" in parts[0] or "サンプル" in parts[0]:
                    info["Label"] = parts[1]
                if "%T" in parts[1]:
                    info["YUnit"] = "%T"
                elif "%R" in parts[1]:
                    info["YUnit"] = "%R"

        tokens = s.replace(",", " ").replace("\t", " ").split()
        nums = []
        for t in tokens:
            try:
                nums.append(float(t))
            except ValueError:
                pass
        if len(nums) >= 2:
            wl.append(nums[0])
            val.append(nums[1])

    if not wl:
        raise ValueError(f"No numeric spectrum data found in {infile}")

    wl = np.asarray(wl, dtype=float)
    val = np.asarray(val, dtype=float)

    uu, idx = np.unique(wl, return_index=True)
    wl = wl[np.sort(idx)]
    val = val[np.sort(idx)]

    order = np.argsort(wl)
    wl = wl[order]
    val = val[order]
    return wl, val, info


def choose_sweep(xs_list, ys_list, sweep_index=0):
    """
    概要:
        複数のスイープデータから指定されたインデックスの単一スイープを選択します。
    引数:
        :param xs_list: 各スイープのX軸（電圧）データのリスト。
        :type xs_list: list[numpy.ndarray]
        :param ys_list: 各スイープのY軸（電流）データのリスト。
        :type ys_list: list[numpy.ndarray]
        :param sweep_index: 選択するスイープのインデックス。
        :type sweep_index: int
    戻り値:
        :returns: 選択されたスイープのX軸データとY軸データ。
        :rtype: tuple[numpy.ndarray, numpy.ndarray]
    例外:
        :raises ValueError: スイープデータが利用できない場合。
        :raises IndexError: sweep_indexが範囲外の場合。
    """
    if len(xs_list) == 0:
        raise ValueError("No sweep data available.")
    idx = int(sweep_index)
    if idx < 0 or idx >= len(xs_list):
        raise IndexError(f"sweep_index={idx} is out of range (0..{len(xs_list)-1})")
    return np.asarray(xs_list[idx], dtype=float), np.asarray(ys_list[idx], dtype=float)


def consolidate_duplicate_x(x, y):
    """
    概要:
        X軸に重複する値がある場合、Y軸の対応する値を平均して重複を解消します。
    引数:
        :param x: X軸データ。
        :type x: numpy.ndarray
        :param y: Y軸データ。
        :type y: numpy.ndarray
    戻り値:
        :returns: 重複が解消されたX軸データと、重複が解消され平均化されたY軸データ。
        :rtype: tuple[numpy.ndarray, numpy.ndarray]
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    order = np.argsort(x)
    xs = x[order]
    ys = y[order]

    ux = []
    uy = []
    i = 0
    n = len(xs)
    while i < n:
        xv = xs[i]
        vals = [ys[i]]
        j = i + 1
        while j < n and abs(xs[j] - xv) < 1e-15:
            vals.append(ys[j])
            j += 1
        ux.append(xv)
        uy.append(float(np.mean(vals)))
        i = j
    return np.asarray(ux), np.asarray(uy)


def smooth_polyfit(y, window_points=5, poly_order=3):
    """
    概要:
        多項式フィッティングを用いたSavitzky-Golay風の平滑化をデータに適用します。
    詳細説明:
        各データポイントを中心に指定されたウィンドウ内のデータに対して多項式フィッティングを行い、
        中心点の値を予測することで平滑化を行います。
    引数:
        :param y: 平滑化するY軸データ。
        :type y: numpy.ndarray
        :param window_points: フィッティングに使用するウィンドウ内のデータポイント数。奇数である必要があります。
        :type window_points: int
        :param poly_order: 多項式フィッティングの次数。
        :type poly_order: int
    戻り値:
        :returns: 平滑化されたY軸データ。
        :rtype: numpy.ndarray
    """
    y = np.asarray(y, dtype=float)
    n = len(y)
    if n < 3:
        return y.copy()

    if window_points < 3:
        window_points = 3
    if window_points % 2 == 0:
        window_points += 1
    if window_points > n:
        window_points = n if (n % 2 == 1) else n - 1
    if poly_order >= window_points:
        poly_order = window_points - 1
    poly_order = max(poly_order, 1)

    half = window_points // 2
    ys = np.empty(n, dtype=float)

    for i in range(n):
        i0 = max(0, i - half)
        i1 = min(n, i + half + 1)
        while (i1 - i0) < window_points:
            if i0 > 0:
                i0 -= 1
            elif i1 < n:
                i1 += 1
            else:
                break
        idx = np.arange(i0, i1, dtype=float)
        yy = y[i0:i1]
        xloc = idx - i
        deg = min(poly_order, len(yy) - 1)
        coeff = np.polyfit(xloc, yy, deg)
        ys[i] = np.polyval(coeff, 0.0)
    return ys


def local_poly_value(x, y, x0, npts=7, order=3):
    """
    概要:
        指定されたX座標の周囲のデータポイントを使用して、局所的な多項式フィッティングを行い、
        x0におけるY値を推定します。
    引数:
        :param x: X軸データ。
        :type x: numpy.ndarray
        :param y: Y軸データ。
        :type y: numpy.ndarray
        :param x0: Y値を推定するX座標。
        :type x0: float
        :param npts: フィッティングに使用するデータポイント数。
        :type npts: int
        :param order: 多項式フィッティングの次数。
        :type order: int
    戻り値:
        :returns: x0における推定されたY値。
        :rtype: float
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    idx = np.argsort(np.abs(x - x0))[:max(2, npts)]
    xs = x[idx]
    ys = y[idx]
    sidx = np.argsort(xs)
    xs = xs[sidx]
    ys = ys[sidx]
    deg = min(order, len(xs) - 1)
    coeff = np.polyfit(xs, ys, deg)
    return float(np.polyval(coeff, x0))


def zero_crossing_x(x, y):
    """
    概要:
        Y値がゼロを横切るX座標を線形補間によって見つけます。
    詳細説明:
        2つの連続するデータポイントの間でY値の符号が変わる点、
        またはY値が最もゼロに近い点のX座標を返します。
    引数:
        :param x: X軸データ。
        :type x: numpy.ndarray
        :param y: Y軸データ。
        :type y: numpy.ndarray
    戻り値:
        :returns: Y値がゼロを横切るX座標。
        :rtype: float
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    order = np.argsort(x)
    x = x[order]
    y = y[order]

    for i in range(len(x) - 1):
        y1, y2 = y[i], y[i + 1]
        if y1 == 0:
            return float(x[i])
        if y1 * y2 < 0:
            return float(x[i] + (0 - y1) * (x[i + 1] - x[i]) / (y2 - y1))
    return float(x[np.argmin(np.abs(y))])


def interpolate_to_common_wavelength(wl_ref, y_ref, wl_target):
    """
    概要:
        参照波長スケール上のデータをターゲット波長スケールに線形補間します。
    引数:
        :param wl_ref: 参照波長データ (nm)。
        :type wl_ref: numpy.ndarray
        :param y_ref: 参照Y軸データ。
        :type y_ref: numpy.ndarray
        :param wl_target: 補間対象のターゲット波長データ (nm)。
        :type wl_target: numpy.ndarray
    戻り値:
        :returns: ターゲット波長スケールに補間されたY軸データ。
        :rtype: numpy.ndarray
    """
    wl_ref = np.asarray(wl_ref, dtype=float)
    y_ref = np.asarray(y_ref, dtype=float)
    wl_target = np.asarray(wl_target, dtype=float)
    return np.interp(wl_target, wl_ref, y_ref, left=np.nan, right=np.nan)


def compute_photon_energy_eV(E_eV, lambda_nm):
    """
    概要:
        フォトンエネルギー（eV）を計算します。
    詳細説明:
        フォトンエネルギーが直接与えられている場合はそれを使用し、
        そうでない場合は波長から計算します。
    引数:
        :param E_eV: フォトンエネルギー (eV)。0より大きい場合はこれを使用。
        :type E_eV: float
        :param lambda_nm: フォトン波長 (nm)。E_eVが0の場合にこれを使用。
        :type lambda_nm: float
    戻り値:
        :returns: 計算されたフォトンエネルギー (eV)。
        :rtype: float
    """
    if E_eV is not None and E_eV > 0:
        return float(E_eV)
    return float(1239.841984 / lambda_nm)


def compute_p0(F0, P0, E_eV, lambda_nm):
    """
    概要:
        入射フォトンエネルギー密度（P0）を計算します。
    詳細説明:
        入射フォトンフラックス（F0）が与えられている場合はそれから計算し、
        そうでない場合は直接指定されたP0を使用します。
    引数:
        :param F0: 入射フォトンフラックス (cm^-2 s^-1)。Noneまたは0の場合はP0を使用します。
        :type F0: float or None
        :param P0: 入射フォトンエネルギー密度 (W/cm^2)。Noneの場合はF0から計算します。
        :type P0: str or None
        :param E_eV: フォトンエネルギー (eV)。F0からP0を計算する場合に使用。
        :type E_eV: float
        :param lambda_nm: フォトン波長 (nm)。F0からP0を計算する場合に使用。
        :type lambda_nm: float
    戻り値:
        :returns: 計算された入射フォトンエネルギー密度 (W/cm^2)。
        :rtype: float or None
    """
    if F0 is not None and F0 > 0:
        E_use = compute_photon_energy_eV(E_eV, lambda_nm)
        return float(F0 * E_use * E_CHARGE)
    return None if P0 is None else float(P0)


def pv_metrics_from_iv(V, I, S):
    """
    概要:
        I-Vデータから太陽電池の主要な性能指標（Voc, Jsc, FF, Pmaxなど）を計算します。
    引数:
        :param V: 電圧データ (V)。
        :type V: numpy.ndarray
        :param I: 電流データ (A)。
        :type I: numpy.ndarray
        :param S: 電極面積 (cm^2)。
        :type S: float
    戻り値:
        :returns: 以下の主要な太陽電池性能指標を含む辞書。
        :rtype: dict
            - Voc_V (float): 開放電圧 (V)。
            - Jsc_A_cm2 (float): 短絡電流密度 (A/cm^2)。
            - Jsc_mA_cm2 (float): 短絡電流密度 (mA/cm^2)。
            - Vop_V (float): 最大出力動作電圧 (V)。
            - Jop_A_cm2 (float): 最大出力動作電流密度 (A/cm^2)。
            - Jop_mA_cm2 (float): 最大出力動作電流密度 (mA/cm^2)。
            - Pmax_W_cm2 (float): 最大出力電力密度 (W/cm^2)。
            - Pmax_mW_cm2 (float): 最大出力電力密度 (mW/cm^2)。
            - FF (float): 曲線因子。
    """
    V = np.asarray(V, dtype=float)
    I = np.asarray(I, dtype=float)
    J = I / S

    Jsc = local_poly_value(V, J, 0.0, npts=7, order=3)
    Voc = zero_crossing_x(V, I)

    Pgen = -V * J
    idx = int(np.argmax(Pgen))
    Vop = float(V[idx])
    Jop = float(J[idx])
    Pmax = float(Pgen[idx])

    denom = abs(Voc * Jsc)
    FF = float(Pmax / denom) if denom > 1e-30 else float("nan")

    return {
        "Voc_V": float(Voc),
        "Jsc_A_cm2": float(Jsc),
        "Jsc_mA_cm2": float(Jsc * 1e3),
        "Vop_V": float(Vop),
        "Jop_A_cm2": float(Jop),
        "Jop_mA_cm2": float(Jop * 1e3),
        "Pmax_W_cm2": float(Pmax),
        "Pmax_mW_cm2": float(Pmax * 1e3),
        "FF": float(FF),
    }


def fit_local_line(x, y, center_idx, npts=7):
    """
    概要:
        指定された中心インデックスの周囲のデータポイントを用いて局所的な線形フィッティングを行います。
    引数:
        :param x: X軸データ。
        :type x: numpy.ndarray
        :param y: Y軸データ。
        :type y: numpy.ndarray
        :param center_idx: フィッティングの中心となるデータポイントのインデックス。
        :type center_idx: int
        :param npts: フィッティングに使用するデータポイント数。
        :type npts: int
    戻り値:
        :returns: フィットされた直線の傾き、Y切片、フィッティングに使用されたXデータ、
                  Yデータ、開始インデックス、終了インデックス（排他的）のタプル。
        :rtype: tuple[float, float, numpy.ndarray, numpy.ndarray, int, int]
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    n = len(x)
    npts = max(2, min(int(npts), n))
    half = npts // 2
    i0 = max(0, center_idx - half)
    i1 = min(n, i0 + npts)
    i0 = max(0, i1 - npts)
    xx = x[i0:i1]
    yy = y[i0:i1]
    a, b = np.polyfit(xx, yy, 1)
    return float(a), float(b), xx, yy, i0, i1


def estimate_rs_tangent_point(V, I):
    """
    概要:
        I-V曲線から直列抵抗（Rs）を推定するための接点と抵抗値を計算します。
    詳細説明:
        順方向バイアス領域でdI/dVが最大となる点を特定し、その点での接線からRsを計算します。
    引数:
        :param V: 電圧データ (V)。
        :type V: numpy.ndarray
        :param I: 電流データ (A)。
        :type I: numpy.ndarray
    戻り値:
        :returns: 直列抵抗推定に関する情報を含む辞書。
        :rtype: dict
            - v_rep (float): 接点電圧 (V)。
            - i_rep (float): 接点電流 (A)。
            - Is (float): 接点電流 (A) (互換性のため)。
            - Rs (float): 推定された直列抵抗 (Ω)。
            - slope (float): 接線の傾き (A/V)。
            - intercept (float): 接線のY切片 (A)。
            - xx (numpy.ndarray): 接線フィッティングに使用されたXデータ。
            - yy (numpy.ndarray): 接線フィッティングに使用されたYデータ。
            - idx (int): 接点のインデックス。
    """
    xu, yu = consolidate_duplicate_x(V, I)
    y_sm = smooth_polyfit(yu, window_points=5, poly_order=3)
    dydv = np.gradient(y_sm, xu)

    pos_idx = np.where(xu > 0.0)[0]
    if len(pos_idx) == 0:
        idx_pos = int(np.argmax(dydv))
    else:
        pad = min(2, max(0, len(pos_idx) // 4))
        cand = pos_idx[pad:len(pos_idx)-pad] if len(pos_idx) - 2 * pad >= 1 else pos_idx
        idx_local = int(np.argmax(dydv[cand]))
        idx_pos = int(cand[idx_local])

    a, b, xx, yy, i0, i1 = fit_local_line(xu, y_sm, idx_pos, npts=7)
    v_rep = float(xu[idx_pos])
    i_rep = float(a * v_rep + b)
    rs = float("inf") if abs(a) < 1e-30 else float(1.0 / a)
    return {"v_rep": v_rep, "i_rep": i_rep, "Is": i_rep, "Rs": rs,
            "slope": float(a), "intercept": float(b), "xx": xx, "yy": yy, "idx": idx_pos}


def estimate_rsh_tangent_point(V, I):
    """
    概要:
        I-V曲線から並列抵抗（Rsh）を推定するための接点と抵抗値を計算します。
    詳細説明:
        逆方向バイアス領域でdI/dVが最小となる点（絶対値）を特定し、
        その点での接線からRshを計算します。
    引数:
        :param V: 電圧データ (V)。
        :type V: numpy.ndarray
        :param I: 電流データ (A)。
        :type I: numpy.ndarray
    戻り値:
        :returns: 並列抵抗推定に関する情報を含む辞書。
        :rtype: dict
            - v_rep (float): 接点電圧 (V)。
            - i_rep (float): 接点電流 (A)。
            - Ish (float): 接点電流 (A) (互換性のため)。
            - Rsh (float): 推定された並列抵抗 (Ω)。
            - slope (float): 接線の傾きの絶対値 (A/V)。
            - intercept (float): 接線のY切片 (A)。
            - slope_fit_raw (float): フィットされた直線の元の傾き (A/V)。
            - intercept_fit_raw (float): フィットされた直線の元のY切片 (A)。
            - xx (numpy.ndarray): 接線フィッティングに使用されたXデータ。
            - yy (numpy.ndarray): 接線フィッティングに使用されたYデータ。
            - idx (int): 接点のインデックス。
    """
    xu, yu = consolidate_duplicate_x(V, I)
    y_sm = smooth_polyfit(yu, window_points=5, poly_order=3)
    dydv = np.gradient(y_sm, xu)

    neg_mask = xu < 0.0
    if np.any(neg_mask):
        dneg = np.abs(dydv[neg_mask])
        idx_neg_local = int(np.argmin(dneg))
        idx_neg = np.where(neg_mask)[0][idx_neg_local]
    else:
        idx_neg = int(np.argmin(np.abs(dydv)))

    a_fit, b_fit, xx, yy, i0, i1 = fit_local_line(xu, y_sm, idx_neg, npts=7)
    v_rep = float(xu[idx_neg])
    i_fit_ref = float(a_fit * v_rep + b_fit)
    a = abs(float(a_fit))
    b = float(i_fit_ref - a * v_rep)
    i_rep = float(a * v_rep + b)
    rsh = float("inf") if abs(a) < 1e-30 else float(1.0 / a)
    return {"v_rep": v_rep, "i_rep": i_rep, "Ish": i_rep, "Rsh": rsh,
            "slope": float(a), "intercept": float(b), "slope_fit_raw": float(a_fit),
            "intercept_fit_raw": float(b_fit), "xx": xx, "yy": yy, "idx": idx_neg}


def estimate_ndiode_representative_point(V, I, T=300.0, Ish=None):
    """
    概要:
        I-V曲線からダイオード因子（ndiode）を推定するための代表点と値を計算します。
    詳細説明:
        順方向バイアス領域でlog(abs(I))の二階微分が最小となる点を特定し、
        その点での一次微分からダイオード因子を計算します。
    引数:
        :param V: 電圧データ (V)。
        :type V: numpy.ndarray
        :param I: 電流データ (A)。
        :type I: numpy.ndarray
        :param T: 温度 (K)。
        :type T: float
        :param Ish: 並列抵抗に流れる電流の推定値。これを考慮してIを調整します。
        :type Ish: float or None
    戻り値:
        :returns: ダイオード因子推定に関する情報を含む辞書。計算できなかった場合はNone。
        :rtype: dict or None
            - v_rep (float): 代表点電圧 (V)。
            - ndiode (float): 推定されたダイオード因子。
            - slope_logI (float): 代表点でのlog(abs(I))の傾き (1/V)。
            - curvature_logI (float): 代表点でのlog(abs(I))の二階微分 (1/V^2)。
    """
    xu, yu = consolidate_duplicate_x(V, I)
    order = np.argsort(xu)
    xu = xu[order]
    yu = yu[order]

    ish_abs = abs(Ish) if Ish is not None else 0.0
    mask = (xu > 0.0) & (np.abs(yu) > ish_abs)
    if np.count_nonzero(mask) < 5:
        return None

    x = xu[mask]
    y = np.abs(yu[mask]) + EPS_I

    logI = np.log(y)
    logI_sm = smooth_polyfit(logI, window_points=5, poly_order=3)
    d1 = np.gradient(logI_sm, x)
    d2 = np.gradient(d1, x)

    good = np.isfinite(d1) & np.isfinite(d2) & (d1 > 0.0)
    if np.count_nonzero(good) < 3:
        return None

    xg = x[good]
    d1g = d1[good]
    d2g = d2[good]

    idx = int(np.argmin(np.abs(d2g)))
    v_rep = float(xg[idx])
    slope = float(d1g[idx])
    ndiode = float(E_CHARGE / (KB * T * slope)) if abs(slope) > 1e-30 else float("nan")
    return {"v_rep": v_rep, "ndiode": ndiode, "slope_logI": slope, "curvature_logI": float(d2g[idx])}


def estimate_initial_params(V, I_meas, T=300.0):
    """
    概要:
        I-Vデータから単一ダイオードモデルの初期パラメータ（I0, ndiode, IPV, Rs, Rsh）を推定します。
    詳細説明:
        estimate_rs_tangent_point、estimate_rsh_tangent_point、
        estimate_ndiode_representative_pointなどの補助関数を用いて、
        各パラメータの初期値を経験的に決定します。
    引数:
        :param V: 電圧データ (V)。
        :type V: numpy.ndarray
        :param I_meas: 測定電流データ (A)。
        :type I_meas: numpy.ndarray
        :param T: 温度 (K)。
        :type T: float
    戻り値:
        :returns: 推定された初期パラメータを含む辞書。
        :rtype: dict
            - I0 (float): 逆方向飽和電流 (A)。
            - ndiode (float): ダイオード因子。
            - IPV (float): 光電流 (A)。
            - Rs (float): 直列抵抗 (Ω)。
            - Rsh (float): 並列抵抗 (Ω)。
            - Ish (float): 並列抵抗に流れる電流 (A)。
            - Vsh (float): Rsh推定における代表点電圧 (V)。
            - Vnd (float): ndiode推定における代表点電圧 (V)。
    """
    xu, yu = consolidate_duplicate_x(V, I_meas)
    order = np.argsort(xu)
    xu = xu[order]
    yu = yu[order]
    y_sm = smooth_polyfit(yu, window_points=5, poly_order=3)

    rs_info = estimate_rs_tangent_point(V, I_meas)
    Rs = abs(float(rs_info["Rs"]))

    rsh_info = estimate_rsh_tangent_point(V, I_meas)
    Rsh = abs(float(rsh_info["Rsh"]))
    Ish = float(rsh_info["Ish"])

    I0 = abs(Ish)
    IPV = -local_poly_value(xu, y_sm, 0.0, npts=7, order=3)

    nd_info = estimate_ndiode_representative_point(V, I_meas, T=T, Ish=Ish)
    if nd_info is None or (not np.isfinite(nd_info["ndiode"])):
        print()
        print("########################################################################")
        print("  Warning!!!: Could not get valid ndiode estimation")
        print("     Use ndiode = 1000.0 instead but don't refer to this value")
        print("########################################################################")
        ndiode = 1000.0
        Vnd = float("nan")
    else:
        ndiode = float(nd_info["ndiode"])
        Vnd = float(nd_info["v_rep"])
    
    return {
        "I0": I0,
        "ndiode": ndiode,
        "IPV": IPV,
        "Rs": Rs,
        "Rsh": Rsh,
        "Ish": Ish,
        "Vsh": float(rsh_info["v_rep"]),
        "Vnd": Vnd,
    }


def analyze_optical(wl_R_nm, R_percent, wl_T_nm, T_percent, d_nm):
    """
    概要:
        反射率（R）と透過率（T）のスペクトルデータから吸収率（A）と吸収係数（α）を計算します。
    詳細説明:
        異なる波長スケールのRとTを共通の波長スケールに補間し、
        RとTのデータに基づいてAとαを計算します。
    引数:
        :param wl_R_nm: 反射率スペクトルの波長データ (nm)。
        :type wl_R_nm: numpy.ndarray
        :param R_percent: 反射率データ (%)。
        :type R_percent: numpy.ndarray
        :param wl_T_nm: 透過率スペクトルの波長データ (nm)。
        :type wl_T_nm: numpy.ndarray
        :param T_percent: 透過率データ (%)。
        :type T_percent: numpy.ndarray
        :param d_nm: PV層の膜厚 (nm)。
        :type d_nm: float
    戻り値:
        :returns: 計算された光学特性値を含む辞書。
        :rtype: dict
            - wl_nm (numpy.ndarray): 共通波長データ (nm)。
            - R (numpy.ndarray): 反射率 (0-1)。
            - Tr (numpy.ndarray): 透過率 (0-1)。
            - A (numpy.ndarray): 吸収率 (0-1)。
            - alpha_cm^-1 (numpy.ndarray): 吸収係数 (cm^-1)。
    """
    wl_common = np.unique(np.concatenate([wl_R_nm, wl_T_nm]))
    wl_common.sort()

    R = interpolate_to_common_wavelength(wl_R_nm, R_percent, wl_common) / 100.0
    Tr = interpolate_to_common_wavelength(wl_T_nm, T_percent, wl_common) / 100.0

    mask = np.isfinite(R) & np.isfinite(Tr)
    wl = wl_common[mask]
    R = np.clip(R[mask], 0.0, 1.0)
    Tr = np.clip(Tr[mask], 0.0, 1.0)

    A = np.clip(1.0 - R - Tr, 0.0, 1.0)

    d_cm = d_nm * 1e-7
    denom = np.clip(1.0 - R, 1e-12, None)
    alpha_cm = -np.log(np.clip(Tr / denom, 1e-12, None)) / max(d_cm, 1e-30)

    return {"wl_nm": wl, "R": R, "Tr": Tr, "A": A, "alpha_cm^-1": alpha_cm}


def quantum_efficiencies(F0, A_abs, JPV_A_cm2, Jsc_A_cm2):
    """
    概要:
        量子効率（EQEとIQE）を計算します。
    詳細説明:
        入射フォトンフラックス、吸収率、光電流（JPV）、短絡電流密度（Jsc）に基づいて、
        生成効率と収集効率を評価します。
    引数:
        :param F0: 入射フォトンフラックス (cm^-2 s^-1)。Noneまたは0の場合はNaNを返します。
        :type F0: float or None
        :param A_abs: 特定のエネルギーにおける吸収率。NoneまたはNaNの場合はIQEをNaNとします。
        :type A_abs: float or None
        :param JPV_A_cm2: 光電流密度 (A/cm^2)。
        :type JPV_A_cm2: float
        :param Jsc_A_cm2: 短絡電流密度 (A/cm^2)。
        :type Jsc_A_cm2: float
    戻り値:
        :returns: 計算された量子効率を含む辞書。
        :rtype: dict
            - EQE_gen (float): 生成外部量子効率。
            - IQE_gen (float): 生成内部量子効率。
            - EQE (float): 外部量子効率。
            - IQE (float): 内部量子効率。
    """
    out = {"EQE_gen": float("nan"), "IQE_gen": float("nan"), "EQE": float("nan"), "IQE": float("nan")}
    if F0 is None or F0 <= 0:
        return out

    denom_ext = E_CHARGE * F0
    out["EQE_gen"] = abs(JPV_A_cm2) / denom_ext
    out["EQE"] = abs(Jsc_A_cm2) / denom_ext

    if A_abs is not None and np.isfinite(A_abs) and A_abs > 0:
        denom_int = E_CHARGE * F0 * A_abs
        out["IQE_gen"] = abs(JPV_A_cm2) / denom_int
        out["IQE"] = abs(Jsc_A_cm2) / denom_int
    return out


def alpha_at_energy(optical, E_eV):
    """
    概要:
        指定されたフォトンエネルギーにおける吸収係数（α）を光学スペクトルデータから補間して取得します。
    引数:
        :param optical: analyze_opticalまたはread_alpha_from_excelから返された光学データ辞書。
        :type optical: dict or None
        :param E_eV: 取得したいフォトンエネルギー (eV)。
        :type E_eV: float
    戻り値:
        :returns: 指定されたエネルギーにおける吸収係数 (cm^-1)。データがない場合はNaN。
        :rtype: float
    例外:
        :raises ValueError: 補間された吸収係数がNaNの場合。
    """
    if optical is None:
        return float("nan")

    wl = np.asarray(optical["wl_nm"], dtype=float)
    alpha = np.asarray(optical["alpha_cm^-1"], dtype=float)
    energy = 1239.841984 / wl
    order = np.argsort(energy)
    energy = energy[order]
    alpha = alpha[order]
    
    alpha_interp = np.interp(E_eV, energy, alpha, left=np.nan, right=np.nan)
    if not np.isfinite(alpha_interp):
        raise ValueError(f"got nan for alpha_at_energy()")

    return float(alpha_interp)

def absorptance_from_alpha(alpha_cm, d_nm):
    """
    概要:
        吸収係数と膜厚から吸収率を計算します。
    詳細説明:
        A = 1 - exp(-alpha * d) の式を用いて計算します。
    引数:
        :param alpha_cm: 吸収係数 (cm^-1)。
        :type alpha_cm: float or None
        :param d_nm: PV層の膜厚 (nm)。
        :type d_nm: float
    戻り値:
        :returns: 計算された吸収率。
        :rtype: float
    例外:
        :raises ValueError: alpha_cmがNoneまたは無効な値の場合。
    """
    if alpha_cm is None or not np.isfinite(alpha_cm):
        raise ValueError(f"got None or invalid value for alpha_cm")
    d_cm = d_nm * 1e-7
    return float(1.0 - np.exp(-alpha_cm * d_cm))


def save_alpha_to_excel(optical, outfile_xlsx, meta=None):
    """
    概要:
        吸収スペクトルデータをExcelファイルに保存します。
    詳細説明:
        alpha_spectrumというシートにフォトンエネルギー、波長、R、Tr、A、alphaのデータを書き込みます。
        summaryシートには膜厚やファイル名、平均値などのメタデータを保存します。
    引数:
        :param optical: analyze_optical関数から返された光学データ辞書。
        :type optical: dict
        :param outfile_xlsx: 保存するExcelファイルのパス。
        :type outfile_xlsx: str
        :param meta: 保存する追加のメタデータ。
        :type meta: dict or None
    戻り値:
        なし
    """
    print(f"[I/O] Save alpha Excel: {outfile_xlsx}")
    wb = Workbook()
    ws = wb.active
    ws.title = "alpha_spectrum"

    header_fill = PatternFill("solid", fgColor="1F4E78")
    header_font = Font(color="FFFFFF", bold=True)

    headers = ["photon_energy_eV", "wavelength_nm", "R", "Tr", "A", "alpha_cm^-1"]
    ws.append(headers)
    for c in ws[1]:
        c.fill = header_fill
        c.font = header_font

    hc_eVnm = 1239.841984
    energy_eV = hc_eVnm / optical["wl_nm"]
    order = np.argsort(energy_eV)

    for i in order:
        ws.append([
            float(energy_eV[i]),
            float(optical["wl_nm"][i]),
            float(optical["R"][i]),
            float(optical["Tr"][i]),
            float(optical["A"][i]),
            float(optical["alpha_cm^-1"][i]),
        ])

    for col in ["A", "B", "C", "D", "E", "F"]:
        ws.column_dimensions[col].width = 18

    ws2 = wb.create_sheet("summary")
    ws2["A1"] = "d_nm"
    ws2["B1"] = meta.get("d_nm", "") if meta else ""
    ws2["A2"] = "R_file"
    ws2["B2"] = meta.get("R", "") if meta else ""
    ws2["A3"] = "Tr_file"
    ws2["B3"] = meta.get("Tr", "") if meta else ""
    ws2["A5"] = "R_mean"
    ws2["B5"] = float(np.mean(optical["R"]))
    ws2["A6"] = "Tr_mean"
    ws2["B6"] = float(np.mean(optical["Tr"]))
    ws2["A7"] = "A_mean"
    ws2["B7"] = float(np.mean(optical["A"]))
    ws2["A8"] = "alpha_mean_cm^-1"
    ws2["B8"] = float(np.mean(optical["alpha_cm^-1"]))

    for col in ["A", "B"]:
        ws2.column_dimensions[col].width = 24

    wb.save(outfile_xlsx)
    print(f"[I/O] Saved alpha Excel: {outfile_xlsx}")


def plot_alpha(optical, outfile=None, pause=False):
    """
    概要:
        吸収係数スペクトル（α vs フォトンエネルギー）をプロットします。
    引数:
        :param optical: analyze_optical関数から返された光学データ辞書。
        :type optical: dict
        :param outfile: プロットを保存するファイルパス。Noneの場合、画面に表示します。
        :type outfile: str or None
        :param pause: プロットを閉じるまで一時停止するかどうか。Trueの場合、plt.show()を呼び出します。
        :type pause: bool
    戻り値:
        :returns: 作成されたMatplotlib Figureオブジェクト。
        :rtype: matplotlib.figure.Figure
    """
    fig, ax = plt.subplots(1, 1, figsize=(6.2, 4.8))
    hc_eVnm = 1239.841984
    energy_eV = hc_eVnm / optical["wl_nm"]

    order = np.argsort(energy_eV)
    x = energy_eV[order]
    y = optical["alpha_cm^-1"][order]

    ax.plot(x, y, "-", linewidth=1.5)
    ax.tick_params(axis='x', labelsize=fontsize)
    ax.tick_params(axis='y', labelsize=fontsize)
    ax.set_xlabel("Photon energy / eV", fontsize=fontsize)
    ax.set_ylabel(r"$\alpha$ / cm$^{-1}$", fontsize=fontsize)
    ax.set_title(r"Absorption coefficient $\alpha(E)$", fontsize=fontsize)
    ax.grid(True)

    fig.tight_layout()
    if outfile:
        print(f"[I/O] Save plot: {outfile}")
        fig.savefig(outfile, dpi=160)

    if pause:
        plt.show()
    else:
        plt.pause(0.01)
    return fig


def estimate_ndiode_curve(V, I, T=300.0):
    """
    概要:
        I-Vデータから電圧に対するダイオード因子（ndiode）の曲線を推定します。
    詳細説明:
        log(abs(I))の一次微分を計算し、それに基づいて各電圧点におけるダイオード因子を導出します。
    引数:
        :param V: 電圧データ (V)。
        :type V: numpy.ndarray
        :param I: 電流データ (A)。
        :type I: numpy.ndarray
        :param T: 温度 (K)。
        :type T: float
    戻り値:
        :returns: ダイオード因子が計算された電圧データと、各電圧に対応する推定されたダイオード因子。
        :rtype: tuple[numpy.ndarray, numpy.ndarray]
    """
    V = np.asarray(V, dtype=float)
    I = np.asarray(I, dtype=float)
    xu, yu = consolidate_duplicate_x(V, I)
    order = np.argsort(xu)
    xu = xu[order]
    yu = yu[order]

    mask = (xu > 0.0) & (yu > 0.0)
    if np.count_nonzero(mask) < 5:
        return np.array([]), np.array([])

    x = xu[mask]
    y = yu[mask]
    logI = np.log(np.abs(y) + EPS_I)
    logI_sm = smooth_polyfit(logI, window_points=5, poly_order=3)
    dlogIdV = np.gradient(logI_sm, x)

    with np.errstate(divide="ignore", invalid="ignore"):
        ndiode = E_CHARGE / (KB * T * dlogIdV)
    good = np.isfinite(ndiode) & (ndiode > 0)
    return x[good], ndiode[good]


def plot_iv_comparison(dark_sweeps, light_sweeps, S, light_metrics,
                       dark_params=None, light_params=None, outfile=None, T=300.0, pause=False):
    """
    概要:
        暗電流I-Vおよび光照射I-Vデータと、関連する推定値や解析結果を比較プロットします。
    詳細説明:
        以下の4つのサブプロットを作成します。
        - log_10(abs(I))-V曲線（暗電流と光電流）。
        - 線形J-V曲線（順方向と逆方向）。RsとRshの接線、Ish(0)点を表示。
        - ndiode-V曲線と代表点。
        - 光起電力出力（-J vs V）。MPP, Voc, Jsc点を表示。
    引数:
        :param dark_sweeps: 暗電流測定のI-Vスイープデータのリスト。
        :type dark_sweeps: list[tuple[numpy.ndarray, numpy.ndarray]]
        :param light_sweeps: 光照射測定のI-Vスイープデータのリスト。
        :type light_sweeps: list[tuple[numpy.ndarray, numpy.ndarray]]
        :param S: デバイスの面積 (cm^2)。
        :type S: float
        :param light_metrics: 光照射下のI-Vデータから計算された光起電力性能指標。
        :type light_metrics: dict
        :param dark_params: 暗電流I-Vデータから推定された初期パラメータ。Rs/Rshの接線表示に使用。
        :type dark_params: dict or None
        :param light_params: 光照射I-Vデータから推定された初期パラメータ。I0/IPVの表示に使用。
        :type light_params: dict or None
        :param outfile: プロットを保存するファイルパス。Noneの場合、画面に表示します。
        :type outfile: str or None
        :param T: 温度 (K)。ndiode曲線の計算に使用。
        :type T: float
        :param pause: プロットを閉じるまで一時停止するかどうか。Trueの場合、plt.show()を呼び出します。
        :type pause: bool
    戻り値:
        :returns: 作成されたMatplotlib Figureオブジェクト。
        :rtype: matplotlib.figure.Figure
    """
    fig = plt.figure(figsize=(13, 8))
    gs = fig.add_gridspec(2, 2, height_ratios=[3, 1.5])

    ax_log = fig.add_subplot(gs[0, 0])
    ax_lin_f = fig.add_subplot(gs[0, 1])
    ax_nd = fig.add_subplot(gs[1, 0], sharex=ax_log)
    ax_pv = fig.add_subplot(gs[1, 1])

    ax_lin_r = ax_lin_f.twinx()

    def style_ax(ax, xlabel=None, ylabel=None, title=None):
        ax.tick_params(axis='x', labelsize=fontsize)
        ax.tick_params(axis='y', labelsize=fontsize)
        if xlabel is not None:
            ax.set_xlabel(xlabel, fontsize=fontsize)
        if ylabel is not None:
            ax.set_ylabel(ylabel, fontsize=fontsize)
        if title is not None:
            ax.set_title(title, fontsize=fontsize)

    def plot_split_linear(ax_f, ax_r, sweeps, label_prefix, color=None):
        first_f = True
        first_r = True
        for V, I in sweeps:
            V = np.asarray(V, dtype=float)
            I = np.asarray(I, dtype=float)

            mask_f = V >= 0.0
            mask_r = V <= 0.0

            if np.any(mask_f):
                Jf = I[mask_f] / S * 1e3
                ax_f.plot(np.abs(V[mask_f]), Jf, "-", linewidth=1.2, color=color,
                          label=label_prefix if first_f else None)
                first_f = False

            if np.any(mask_r):
                Jr = I[mask_r] / S * 1e3
                ax_r.plot(np.abs(V[mask_r]), Jr, "--", linewidth=1.2, color=color,
                          label=f"{label_prefix} (rev)" if first_r else None)
                first_r = False

    plot_split_linear(ax_lin_f, ax_lin_r, dark_sweeps, "dark")
    plot_split_linear(ax_lin_f, ax_lin_r, light_sweeps, "light")

    if dark_params is not None and len(dark_sweeps) > 0:
        Vd, Id = dark_sweeps[0]
        rs_info = estimate_rs_tangent_point(Vd, Id)
        a_rs = rs_info["slope"] / S * 1e3
        b_rs = rs_info["intercept"] / S * 1e3
        v_rep_rs = rs_info["v_rep"]
        i_rep_rs = rs_info["i_rep"] / S * 1e3

        if abs(a_rs) > 1e-30:
            v_zero = float(-b_rs / a_rs)
            x0 = max(0.0, min(abs(v_rep_rs), abs(v_zero)))
            x1 = max(abs(v_rep_rs), abs(v_zero))
            x_tan_rs = np.array([x0, x1])
            y_tan_rs = a_rs * x_tan_rs + b_rs
            ax_lin_f.plot(x_tan_rs, y_tan_rs, "-", linewidth=0.9, label="Rs tangent")
            ax_lin_f.plot(abs(v_rep_rs), i_rep_rs, "o", markersize=5, label="Rs point")

    if dark_params is not None and len(dark_sweeps) > 0:
        Vd, Id = dark_sweeps[0]
        tang = estimate_rsh_tangent_point(Vd, Id)
        a = tang["slope"] / S * 1e3
        b = tang["intercept"] / S * 1e3
        v_rep = tang["v_rep"]
        i_rep = tang["i_rep"] / S * 1e3

        x_rev_all = []
        for VV, II in dark_sweeps:
            VV = np.asarray(VV, dtype=float)
            mask_r_all = VV <= 0.0
            if np.any(mask_r_all):
                x_rev_all.extend(list(np.abs(VV[mask_r_all])))
        if len(x_rev_all) == 0:
            x_tan = np.array([0.0, abs(v_rep)])
        else:
            x_tan = np.array([0.0, float(np.max(x_rev_all))])

        y_tan = (-a) * x_tan + b
        ax_lin_r.plot(x_tan, y_tan, "-", linewidth=0.9, label="Rsh tangent", zorder=5)
        ax_lin_r.plot(abs(v_rep), i_rep, "o", markersize=5, label="Rsh point", zorder=6)
        ish0 = float(b)
        ax_lin_r.plot(0.0, ish0, "s", markersize=5, label="Ish(0)", zorder=6)

    style_ax(ax_lin_f, xlabel="|V| / V", ylabel="J forward / mA cm$^{-2}$, V >= 0",
             title="Linear J-V (forward / reverse separated)")
    ax_lin_r.tick_params(axis='y', labelsize=fontsize)
    ax_lin_r.set_ylabel("J reverse / mA cm$^{-2}$, V <= 0", fontsize=fontsize)
    ax_lin_f.grid(True)

    h1, l1 = ax_lin_f.get_legend_handles_labels()
    h2, l2 = ax_lin_r.get_legend_handles_labels()
    ax_lin_f.legend(h1 + h2, l1 + l2, fontsize=8, loc="best")

    for i, (V, I) in enumerate(dark_sweeps):
        ax_log.plot(V, np.log10(np.abs(I) + EPS_I), "-", linewidth=1.2, label="dark" if i == 0 else None)
    for i, (V, I) in enumerate(light_sweeps):
        ax_log.plot(V, np.log10(np.abs(I) + EPS_I), "-", linewidth=1.2, label="light" if i == 0 else None)
    if light_params is not None:
        ax_log.axhline(np.log10(abs(light_params["I0"]) + EPS_I), linestyle=":", linewidth=1.0, label="-I0")
        ax_log.axhline(np.log10(abs(light_params["I0"] + light_params["IPV"]) + EPS_I),
                       linestyle=":", linewidth=1.0, label="-I0-IPV")
    style_ax(ax_log, xlabel="V / V", ylabel=r"log$_{10}$(|I|)", title="Dark / Illuminated log10(|I|)-V")
    ax_log.grid(True)
    ax_log.legend(fontsize=8)

    Vdark0, Idark0 = dark_sweeps[0]
    Vn, nn = estimate_ndiode_curve(Vdark0, Idark0, T=T)
    if len(Vn) > 0:
        ax_nd.plot(Vn, nn, "-", linewidth=1.3, label="ndiode(V)")

    ish_for_nd = dark_params.get("Ish", None) if dark_params is not None else None
    nd_rep = estimate_ndiode_representative_point(Vdark0, Idark0, T=T, Ish=ish_for_nd)
    if nd_rep is not None and np.isfinite(nd_rep["ndiode"]):
        ax_nd.plot(nd_rep["v_rep"], nd_rep["ndiode"], "o", markersize=6, label="ndiode point")

    style_ax(ax_nd, xlabel="V / V", ylabel="ndiode", title="Estimated ndiode vs V")
    ax_nd.grid(True)
    ax_nd.legend(fontsize=8)

    V0, I0 = light_sweeps[0]
    J0_mA = I0 / S * 1e3
    ypv = -J0_mA
    ax_pv.plot(V0, ypv, "-", linewidth=1.5, label="-J (light)")
    ax_pv.plot(light_metrics["Vop_V"], -light_metrics["Jop_mA_cm2"], "o", markersize=7, label="MPP")
    ax_pv.plot(light_metrics["Voc_V"], 0.0, "s", markersize=7, label="Voc")
    ax_pv.plot(0.0, -light_metrics["Jsc_mA_cm2"], "^", markersize=7, label="Jsc")

    style_ax(ax_pv, xlabel="V / V", ylabel="-J / mA/cm$^2$", title="Photovoltaic output")
    ymax = max(float(np.max(ypv)), 0.0, float(-light_metrics["Jsc_mA_cm2"]))
    ax_pv.set_ylim(0.0, ymax if ymax > 0 else 1.0)
    ax_pv.grid(True)
    ax_pv.legend(fontsize=8)

    fig.tight_layout()
    if outfile:
        print(f"[I/O] Save plot: {outfile}")
        fig.savefig(outfile, dpi=160)

    if pause:
        plt.show()
    else:
        plt.pause(0.01)
    return fig


def exec_alpha(args):
    """
    概要:
        alphaまたはmake_alphaモードでプログラムを実行するメイン関数。
    詳細説明:
        光学スペクトルファイル（RとTr）を読み込み、吸収率と吸収係数を計算し、
        結果をプロットしてExcelファイルに保存します。
    引数:
        :param args: コマンドライン引数。
        :type args: argparse.Namespace
    戻り値:
        なし
    """

#    E_use = compute_photon_energy_eV(args.E, args.lambda_nm)
#    P0_use = compute_p0(args.F0, args.P0, args.E, args.lambda_nm)
#    print_args_and_derived(args, E_use, P0_use)
    print_args_and_derived(args)
    
    wl_R, R_percent, info_R = read_optical_spectrum(args.R)
    wl_T, T_percent, info_T = read_optical_spectrum(args.Tr)

    optical = analyze_optical(wl_R, R_percent, wl_T, T_percent, args.d)

    print()
    print("=== make_alpha analysis ===")
    print(f"d_nm          : {args.d:.16g} nm")
    print(f"R_mean        : {np.mean(optical['R']):.16g}")
    print(f"Tr_mean       : {np.mean(optical['Tr']):.16g}")
    print(f"A_mean        : {np.mean(optical['A']):.16g}")
    print(f"alpha_mean    : {np.mean(optical['alpha_cm^-1']):.16g} cm^-1")

    print()
    outxlsx = f"{args.outprefix}_alpha.xlsx"
    if args.mode == 'make_alpha':
        save_alpha_to_excel(optical, outxlsx, meta={"d_nm": args.d, "R": args.R, "Tr": args.Tr})
        print(f"Saving Excel  : {outxlsx}")

    outplot = f"{args.outprefix}_alpha_E.png"
    print(f"Saving plot   : {outplot}")

    plot_alpha(optical, outfile=outplot, pause = False)
    input("\nPress ENTER to terminate\n")


def exec_analyze(args):
    """
    概要:
        analyzeモードでプログラムを実行するメイン関数。
    詳細説明:
        暗電流および光照射I-Vファイルを読み込み、初期パラメータと光起電力性能指標を推定し、
        量子効率を計算します。光学データが利用可能な場合はそれも考慮に入れます。
        結果を標準出力に表示し、比較プロットを作成します。
    引数:
        :param args: コマンドライン引数。
        :type args: argparse.Namespace
    戻り値:
        なし
    例外:
        :raises ValueError: alpha_at_energy 関数から NaN が返された場合。
    """
    print()
    print("[MODE] analyze")
    P0_use = compute_p0(args.F0, args.P0, args.E, args.lambda_nm)
    E_use = compute_photon_energy_eV(args.E, args.lambda_nm)
    print_args_and_derived(args, E_use, P0_use)

    dark_xs, dark_ys, dark_info = read_data(args.dark)
    light_xs, light_ys, light_info = read_data(args.light)

    dark_sweeps = list(zip(dark_xs, dark_ys))
    light_sweeps = list(zip(light_xs, light_ys))

    Vd, Id = choose_sweep(dark_xs, dark_ys, args.sweep_dark)
    Vl, Il = choose_sweep(light_xs, light_ys, args.sweep_light)

    dark_params = estimate_initial_params(Vd, Id, T=args.T)
    dark_params["IPV"] = 0.0
    light_params = estimate_initial_params(Vl, Il, T=args.T)

    dark_metrics = pv_metrics_from_iv(Vd, Id, args.S)
    light_metrics = pv_metrics_from_iv(Vl, Il, args.S)

    optical = None
    alpha_at_E = float("nan")
    A_at_E = float("nan")
    if args.alpha:
        optical = read_alpha_from_excel(args.alpha)

    if optical is not None:
        alpha_at_E = alpha_at_energy(optical, E_use)
        if alpha_at_E is None:
            raise ValueError(f"Got nan for alpha_at_E()")

        A_at_E = absorptance_from_alpha(alpha_at_E, args.d)

    eta = float("nan")
    if P0_use is not None and P0_use > 0:
        eta = light_metrics["Pmax_W_cm2"] / P0_use

    JPV_A_cm2 = light_params["IPV"] / args.S
    qe = quantum_efficiencies(args.F0, A_at_E, JPV_A_cm2, light_metrics["Jsc_A_cm2"])

    print("=== Device constants ===")
    print(f"d            = {args.d:.6f} nm")
    print(f"S            = {args.S:.8f} cm^2")
    print(f"T            = {args.T:.6f} K")
    print()

    print("=== Photon input ===")
    print(f"lambda       = {args.lambda_nm:.16g} nm")
    print(f"E            = {E_use:.16g} eV")
    print(f"F0           = {args.F0:.16g} cm^-2 s^-1" if args.F0 is not None else "F0           = None")
    print(f"P0           = {P0_use*1e3:.16g} mW/cm^2" if P0_use is not None else "P0           = None")
    print()

    print("=== Dark IV estimated parameters ===")
    for k in PARAM_NAMES:
        print(f"{k:8s} = {dark_params[k]:.16g}")
    rsh_info = estimate_rsh_tangent_point(Vd, Id)
    rs_info = estimate_rs_tangent_point(Vd, Id)
    print(f"{'Vs':8s} = {rs_info['v_rep']:.16g}")
    print(f"{'Is':8s} = {rs_info['i_rep']:.16g}")
    print(f"{'Vsh':8s} = {rsh_info['v_rep']:.16g}")
    print(f"{'Ish':8s} = {rsh_info['Ish']:.16g}")
    print(f"{'Ish(0)':8s} = {rsh_info['intercept']:.16g}")
    print()

    print("=== Illuminated IV estimated parameters ===")
    for k in PARAM_NAMES:
        print(f"{k:8s} = {light_params[k]:.16g}")
    print()

    print("=== Photovoltaic metrics ===")
    print(f"Voc        = {light_metrics['Voc_V']:.16g} V")
    print(f"Jsc        = {light_metrics['Jsc_mA_cm2']:.16g} mA/cm^2")
    print(f"Vop        = {light_metrics['Vop_V']:.16g} V")
    print(f"Jop        = {light_metrics['Jop_mA_cm2']:.16g} mA/cm^2")
    print(f"FF         = {light_metrics['FF']:.16g}")
    print(f"Pmax       = {light_metrics['Pmax_mW_cm2']:.16g} mW/cm^2")
    print(f"Efficiency = {eta:.16g}" if np.isfinite(eta) else "Efficiency = nan")
    print()

    if optical is not None:
        print("=== Optical metrics ===")
        print(f"alpha(E)          = {alpha_at_E:.16g} cm^-1")
        print(f"absorption length = {args.d} nm")
        print(f"A(E)              = {A_at_E:.16g}")
        print()

    print("=== Quantum efficiencies ===")
    print(f"JPV         = {JPV_A_cm2*1e3:.16g} mA/cm^2")
    print(f"EQE_gen     = {qe['EQE_gen']:.16g}")
    print(f"IQE_gen     = {qe['IQE_gen']:.16g}")
    print(f"EQE         = {qe['EQE']:.16g}")
    print(f"IQE         = {qe['IQE']:.16g}")

    print()
    outplot = f"{args.outprefix}_iv.png"
    plot_iv_comparison(
        dark_sweeps, light_sweeps, args.S, light_metrics,
        dark_params=dark_params, light_params=light_params,
        outfile=outplot, T=args.T,
        pause=False
    )
    print(f"Saved plot: {outplot}")
    input("\nPress ENTER to terminate\n")


_original_print = None
_redirect_fp = None
def dual_print(*args, **kwargs):
    """
    概要:
        標準出力とログファイルの両方にメッセージを出力します。
    詳細説明:
        通常のprint関数の動作を模倣し、標準出力に加えて、
        _redirect_fpに設定されたファイルオブジェクトにも出力します。
    引数:
        :param args: printに渡される位置引数。
        :type args: tuple
        :param kwargs: printに渡されるキーワード引数。
        :type kwargs: dict
    戻り値:
        なし
    """
    _original_print(*args, **kwargs)

    file_kwargs = dict(kwargs)
    file_kwargs.pop("file", None)   # file指定があっても無視してログファイルへ
    file_kwargs.setdefault("flush", True)

    _original_print(*args, file=_redirect_fp, **file_kwargs)

def main():
    """
    概要:
        プログラムのエントリーポイント。
    詳細説明:
        コマンドライン引数を初期化し、指定されたモードに基づいて
        対応する実行関数を呼び出します。
        実行ログは自動的にファイルにリダイレクトされます。
    引数:
        なし
    戻り値:
        なし
    """

    global _original_print, _redirect_fp

    print()
    args, parser = initialize()

    _original_print = builtins.print
    if "alpha" in args.mode:
        logfile = f"{args.outprefix}_alpha.txt"
    else:
        logfile = f"{args.outprefix}_{args.mode}.txt"
    print(f"Open log file [{logfile}]")
    _redirect_fp = open(logfile, "w", encoding="utf-8")
    builtins.print = dual_print

    if args.mode == "alpha" or args.mode == "make_alpha":
        exec_alpha(args)
    elif args.mode == "analyze":
        exec_analyze(args)
    else:
        parser.error(f"Unsupported mode: {args.mode}")


if __name__ == "__main__":
    try:
        main()
    except Exception:
        traceback.print_exc()
        sys.exit(1)
    finally:
        if _redirect_fp: _redirect_fp.close()