#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ :doc:`pvanalyze_usage` PV特性評価ツール mode=alpha 反射率/透過率ファイルと膜厚から吸収スペクトルを計算し、 プロットとExcel保存を行います。 mode=analyze I-Vデータからパラメータ推定、発電特性解析、量子効率計算、プロットを行います。 """ 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" """ DARK_IV_FILE = None PV_IV_FILE = None R_FILE = None T_FILE = None ALPHA_FILE = None def initialize(): """ コマンドライン引数を解析し、プログラムの初期設定を行います。 詳細説明: `argparse`モジュールを使用して、実行モード、ファイルパス、温度、膜厚などの パラメータを定義し、ユーザー入力から値をパースします。 :returns: `args` (argparse.Namespace): 解析された引数オブジェクト。 `parser` (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: (argparse.Namespace): 解析された引数オブジェクト。 :param E_use: (float): 使用されるフォトンエネルギー (eV)。 :param P0_use: (float or None): 使用される入射光パワー密度 (W/cm^2)。 :returns: なし """ 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: (str): 読み込むCSVファイルのパス。 :param xmin: (float, optional): X軸（電圧）の最小値。これより小さい値はスキップされます。デフォルトはNone。 :param xmax: (float, optional): X軸（電圧）の最大値。これより大きい値はスキップされます。デフォルトはNone。 :param ndataskip: (int, optional): データポイントをスキップする間隔。0の場合、スキップしません。デフォルトは0。 :returns: `xs_list` (list[numpy.ndarray]): 各スイープのX軸（電圧）データのリスト。 `ys_list` (list[numpy.ndarray]): 各スイープのY軸（電流）データのリスト。 `inf` (dict): ファイルから抽出されたメタデータ。 """ 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: (str): 読み込むExcelファイルのパス。 :returns: `dict`: 波長、反射率、透過率、吸収率、吸収係数を含む辞書。 - `"wl_nm"` (numpy.ndarray): 波長 (nm)。 - `"R"` (numpy.ndarray): 反射率。 - `"Tr"` (numpy.ndarray): 透過率。 - `"A"` (numpy.ndarray): 吸収率。 - `"alpha_cm^-1"` (numpy.ndarray): 吸収係数 (cm^-1)。 """ 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: (str): 読み込む光学スペクトルファイルのパス。 :returns: `wl` (numpy.ndarray): 波長 (nm)。 `val` (numpy.ndarray): スペクトル値 (RまたはTのパーセンテージ)。 `info` (dict): ファイルから抽出されたメタデータ（ファイル名、ラベル、Y軸単位など）。 """ 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: (list[numpy.ndarray]): 各スイープのX軸（電圧）データのリスト。 :param ys_list: (list[numpy.ndarray]): 各スイープのY軸（電流）データのリスト。 :param sweep_index: (int, optional): 選択するスイープのインデックス。デフォルトは0。 :returns: `x` (numpy.ndarray): 選択されたスイープのX軸データ。 `y` (numpy.ndarray): 選択されたスイープのY軸データ。 """ 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: (numpy.ndarray): X軸データ。 :param y: (numpy.ndarray): Y軸データ。 :returns: `ux` (numpy.ndarray): 重複が解消されたX軸データ。 `uy` (numpy.ndarray): 重複が解消され、平均化されたY軸データ。 """ 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: (numpy.ndarray): 平滑化するY軸データ。 :param window_points: (int, optional): フィッティングに使用するウィンドウ内のデータポイント数。奇数である必要があります。デフォルトは5。 :param poly_order: (int, optional): 多項式フィッティングの次数。デフォルトは3。 :returns: `ys` (numpy.ndarray): 平滑化されたY軸データ。 """ 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: (numpy.ndarray): X軸データ。 :param y: (numpy.ndarray): Y軸データ。 :param x0: (float): Y値を推定するX座標。 :param npts: (int, optional): フィッティングに使用するデータポイント数。デフォルトは7。 :param order: (int, optional): 多項式フィッティングの次数。デフォルトは3。 :returns: `float`: `x0`における推定されたY値。 """ 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: (numpy.ndarray): X軸データ。 :param y: (numpy.ndarray): Y軸データ。 :returns: `float`: Y値がゼロを横切るX座標。 """ 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: (numpy.ndarray): 参照波長データ (nm)。 :param y_ref: (numpy.ndarray): 参照Y軸データ。 :param wl_target: (numpy.ndarray): 補間対象のターゲット波長データ (nm)。 :returns: `numpy.ndarray`: ターゲット波長スケールに補間されたY軸データ。 """ 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: (float): フォトンエネルギー (eV)。0より大きい場合はこれを使用。 :param lambda_nm: (float): フォトン波長 (nm)。`E_eV`が0の場合にこれを使用。 :returns: `float`: 計算されたフォトンエネルギー (eV)。 """ 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: (float or None): 入射フォトンフラックス (cm^-2 s^-1)。Noneの場合はP0を使用。 :param P0: (float or None): 入射フォトンエネルギー密度 (W/cm^2)。Noneの場合はF0から計算。 :param E_eV: (float): フォトンエネルギー (eV)。F0からP0を計算する場合に使用。 :param lambda_nm: (float): フォトン波長 (nm)。F0からP0を計算する場合に使用。 :returns: `float or None`: 計算された入射フォトンエネルギー密度 (W/cm^2)。 """ if F0 is not None: 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: (numpy.ndarray): 電圧データ (V)。 :param I: (numpy.ndarray): 電流データ (A)。 :param S: (float): 電極面積 (cm^2)。 :returns: `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: (numpy.ndarray): X軸データ。 :param y: (numpy.ndarray): Y軸データ。 :param center_idx: (int): フィッティングの中心となるデータポイントのインデックス。 :param npts: (int, optional): フィッティングに使用するデータポイント数。デフォルトは7。 :returns: `a` (float): フィットされた直線の傾き。 `b` (float): フィットされた直線のY切片。 `xx` (numpy.ndarray): フィッティングに使用されたXデータ。 `yy` (numpy.ndarray): フィッティングに使用されたYデータ。 `i0` (int): フィッティングに使用されたデータの開始インデックス。 `i1` (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: (numpy.ndarray): 電圧データ (V)。 :param I: (numpy.ndarray): 電流データ (A)。 :returns: `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: (numpy.ndarray): 電圧データ (V)。 :param I: (numpy.ndarray): 電流データ (A)。 :returns: `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(I)の二階微分が最小となる点を特定し、 その点での一次微分からダイオード因子を計算します。 :param V: (numpy.ndarray): 電圧データ (V)。 :param I: (numpy.ndarray): 電流データ (A)。 :param T: (float, optional): 温度 (K)。デフォルトは300.0。 :param Ish: (float, optional): 並列抵抗に流れる電流の推定値。これを考慮してIを調整します。デフォルトはNone。 :returns: `dict or None`: ダイオード因子推定に関する情報。計算できなかった場合はNone。 - `"v_rep"` (float): 代表点電圧 (V)。 - `"ndiode"` (float): 推定されたダイオード因子。 - `"slope_logI"` (float): 代表点でのlog(I)の傾き (1/V)。 - `"curvature_logI"` (float): 代表点でのlog(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: (numpy.ndarray): 電圧データ (V)。 :param I_meas: (numpy.ndarray): 測定電流データ (A)。 :param T: (float, optional): 温度 (K)。デフォルトは300.0。 :returns: `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) 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: (numpy.ndarray): 反射率スペクトルの波長データ (nm)。 :param R_percent: (numpy.ndarray): 反射率データ (%)。 :param wl_T_nm: (numpy.ndarray): 透過率スペクトルの波長データ (nm)。 :param T_percent: (numpy.ndarray): 透過率データ (%)。 :param d_nm: (float): PV層の膜厚 (nm)。 :returns: `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: (float or None): 入射フォトンフラックス (cm^-2 s^-1)。Noneまたは0の場合はNaNを返します。 :param A_abs: (float or None): 特定のエネルギーにおける吸収率。NoneまたはNaNの場合はIQEをNaNとします。 :param JPV_A_cm2: (float): 光電流密度 (A/cm^2)。 :param Jsc_A_cm2: (float): 短絡電流密度 (A/cm^2)。 :returns: `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: (dict or None): `analyze_optical`または`read_alpha_from_excel`から返された光学データ辞書。 :param E_eV: (float): 取得したいフォトンエネルギー (eV)。 :returns: `float`: 指定されたエネルギーにおける吸収係数 (cm^-1)。データがない場合は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: (float or None): 吸収係数 (cm^-1)。 :param d_nm: (float): PV層の膜厚 (nm)。 :returns: `float`: 計算された吸収率。 """ 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: (dict): `analyze_optical`関数から返された光学データ辞書。 :param outfile_xlsx: (str): 保存するExcelファイルのパス。 :param meta: (dict, optional): 保存する追加のメタデータ。デフォルトはNone。 :returns: なし """ 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: (dict): `analyze_optical`関数から返された光学データ辞書。 :param outfile: (str, optional): プロットを保存するファイルパス。Noneの場合、画面に表示します。デフォルトはNone。 :returns: `matplotlib.figure.Figure`: 作成されたMatplotlib 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(I)の一次微分を計算し、それに基づいて各電圧点におけるダイオード因子を導出します。 :param V: (numpy.ndarray): 電圧データ (V)。 :param I: (numpy.ndarray): 電流データ (A)。 :param T: (float, optional): 温度 (K)。デフォルトは300.0。 :returns: `Vn` (numpy.ndarray): ダイオード因子が計算された電圧データ。 `nn` (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(|I|)-V曲線（暗電流と光電流）。 - 線形J-V曲線（順方向と逆方向）。RsとRshの接線、Ish(0)点を表示。 - ndiode-V曲線と代表点。 - 光起電力出力（-J vs V）。MPP, Voc, Jsc点を表示。 :param dark_sweeps: (list[tuple[numpy.ndarray, numpy.ndarray]]): 暗電流測定のI-Vスイープデータのリスト。 :param light_sweeps: (list[tuple[numpy.ndarray, numpy.ndarray]]): 光照射測定のI-Vスイープデータのリスト。 :param S: (float): デバイスの面積 (cm^2)。 :param light_metrics: (dict): 光照射下のI-Vデータから計算された光起電力性能指標。 :param dark_params: (dict, optional): 暗電流I-Vデータから推定された初期パラメータ。Rs/Rshの接線表示に使用。デフォルトはNone。 :param light_params: (dict, optional): 光照射I-Vデータから推定された初期パラメータ。I0/IPVの表示に使用。デフォルトはNone。 :param outfile: (str, optional): プロットを保存するファイルパス。Noneの場合、画面に表示します。デフォルトはNone。 :param T: (float, optional): 温度 (K)。ndiode曲線の計算に使用。デフォルトは300.0。 :returns: `matplotlib.figure.Figure`: 作成されたMatplotlib 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とT）を読み込み、吸収率と吸収係数を計算し、 結果をプロットしてExcelファイルに保存します。 :param args: (argparse.Namespace): コマンドライン引数。 :returns: なし """ # 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) print() print("=== make_alpha analysis ===") print(f"d_nm : {args.d:.16g} nm") if args.R and args.Tr: 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(f"R_mean : {np.mean(optical['R']):.16g}") print(f"Tr_mean : {np.mean(optical['Tr']):.16g}") else: optical = read_alpha_from_excel(args.alpha) 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: (argparse.Namespace): コマンドライン引数。 :returns: なし """ 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): _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(): """ プログラムのエントリーポイント。 詳細説明: コマンドライン引数を初期化し、指定されたモードに基づいて 対応する実行関数を呼び出します。 :returns: なし """ 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()