pvanalyze.py ダウンロード/コピー
pvanalyze.py
pvanalyze.py
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3
4"""
5PV特性評価ツールスクリプト。
6
7概要:
8 太陽電池の電気的・光学的特性を評価するツールです。
9詳細説明:
10 mode=alpha または mode=make_alpha では、反射率と透過率のスペクトルから
11 吸収スペクトルを計算し、プロットおよびExcelへの保存を行います。
12 mode=analyze では、I-Vデータからパラメータ推定、発電特性解析、量子効率計算、
13 および結果のプロットを行います。
14主な機能:
15 - 反射率/透過率から吸収スペクトルの計算と保存
16 - I-Vデータからの太陽電池パラメータ（I0, ndiode, IPV, Rs, Rsh）推定
17 - 発電特性（Voc, Jsc, FF, Pmaxなど）の解析
18 - 量子効率（EQE, IQE）の計算
19 - 解析結果のプロットとExcelファイルへの保存
20関連リンク:
21 pvanalyze_usage
22"""
23
24import sys
25import argparse
26import builtins
27from pathlib import Path
28import traceback
29import csv
30import math
31
32import matplotlib.pyplot as plt
33import numpy as np
34from openpyxl import Workbook, load_workbook
35from openpyxl.styles import Font, PatternFill
36
37
38KB = 1.380649e-23
39E_CHARGE = 1.602176634e-19
40H = 6.62607015e-34
41C = 2.99792458e8
42
43PARAM_NAMES = ["I0", "ndiode", "IPV", "Rs", "Rsh"]
44EPS_I = 1.0e-15
45
46DPV_NM = 46.837
47AREA_MM2 = 0.5 * 0.5 * math.pi
48AREA_CM2 = AREA_MM2 * 0.01 # 1 mm^2 = 0.01 cm^2
49PHOTON_FLUX = 1.95804e18 # cm^-2 s^-1
50PHOTON_NM = 1363.0 # nm
51PIN = "" # W/cm^2, if None and F0 is not None then compute from hc/lambda
52
53fontsize = 16
54
55DARK_IV_FILE = "I_V SweepTest SMU1 [(1) ; 2026_03_24 14_15_24]-SY251213-1-LD0.0V.csv"
56PV_IV_FILE = "I_V SweepTest SMU1 [(10) ; 2026_03_24 14_18_37]-SY251213-1-LD1.8V.csv"
57R_FILE = "SY251024-1-Bi4O6S21-STO001-225oC-10mJ_R(UDS).txt"
58T_FILE = "SY251024-1-Bi4O6S21-STO001-225oC-10mJ_T(UDS).txt"
59ALPHA_FILE = "alpha_Bi2OS2.xlsx"
60
61
62def initialize():
63 """
64 概要:
65 コマンドライン引数を解析し、プログラムの初期設定を行います。
66 詳細説明:
67 argparseモジュールを使用して、実行モード、ファイルパス、温度、膜厚などの
68 パラメータを定義し、ユーザー入力から値をパースします。
69 引数:
70 なし
71 戻り値:
72 :returns: 解析された引数オブジェクトと引数パーサーオブジェクトのタプル。
73 :rtype: tuple[argparse.Namespace, argparse.ArgumentParser]
74 """
75 parser = argparse.ArgumentParser(description="PV characterization tool with alpha/analyze modes.")
76 parser.add_argument("--mode", default="analyze", choices=["alpha", "make_alpha", "analyze"],
77 help="Execution mode")
78 parser.add_argument("--dark", default=DARK_IV_FILE, help="Dark IV CSV file")
79 parser.add_argument("--light", default=PV_IV_FILE, help="Illuminated IV CSV file")
80 parser.add_argument("--R", default=R_FILE, help="Reflectance spectrum file")
81 # NOTE: --T is already requested for temperature, so transmittance file uses --Tr.
82 parser.add_argument("--Tr", default=T_FILE, help="Transmittance spectrum file")
83 parser.add_argument("--eq", default=None, help="Equivalent spectral irradiance file")
84 parser.add_argument("--alpha", default=ALPHA_FILE, help="Excel file for absorption spectrum / alpha")
85 parser.add_argument("--T", type=float, default=300.0, help="Temperature (K)")
86 parser.add_argument("--d", type=float, default=DPV_NM, help="PV layer thickness (nm)")
87 parser.add_argument("--sweep_dark", type=int, default=0, help="Sweep index for dark analysis")
88 parser.add_argument("--sweep_light", type=int, default=0, help="Sweep index for light analysis")
89 parser.add_argument("--F0", type=float, default=PHOTON_FLUX, help="Incident photon flux (cm^-2 s^-1)")
90 parser.add_argument("--P0", type=str, default=PIN, help="Incident photon energy density (W/cm^2)")
91 parser.add_argument("--E", type=float, default=0.0, help="Photon energy (eV). If 0, use lambda.")
92 parser.add_argument("--lambda_nm", type=float, default=PHOTON_NM,
93 help="Photon wavelength (nm), used when E=0")
94 parser.add_argument("--S", type=float, default=AREA_CM2, help="Electrode area (cm^2)")
95 parser.add_argument("--outprefix", default="pvanalyze", help="Output file prefix")
96 args = parser.parse_args()
97 return args, parser
98
99
100def print_args_and_derived(args, E_use=None, P0_use=None):
101 """
102 概要:
103 解析されたコマンドライン引数と派生値を標準出力に表示します。
104 引数:
105 :param args: 解析された引数オブジェクト。
106 :type args: argparse.Namespace
107 :param E_use: 使用されるフォトンエネルギー (eV)。
108 :type E_use: float or None
109 :param P0_use: 使用される入射光パワー密度 (W/cm^2)。
110 :type P0_use: float or None
111 戻り値:
112 なし
113 例外:
114 :raises ValueError: フォトン波長 lambda_nm が0以下の場合。
115 :raises ValueError: フォトンエネルギー E が0未満の場合。
116 :raises ValueError: 入射フォトンフラックス F0 が0以下の場合。
117 :raises ValueError: 入射光パワー密度 P0 が0以下の場合。
118 :raises ValueError: 温度 T が0以下の場合。
119 :raises ValueError: 電極面積 S が0以下の場合。
120 :raises ValueError: PV層の膜厚 d が0.1nm以下の場合。
121 :raises ValueError: compute_p0 関数から None が返された場合。
122 """
123 print("=== Arguments ===")
124 print(f"mode = {args.mode}")
125 if args.mode == "make_alpha":
126 print(f"R = {args.R}")
127 print(f"Tr = {args.Tr}")
128 else:
129 print(f"alpha = {args.alpha}")
130 print(f"d = {args.d} nm")
131 print(f"outprefix = {args.outprefix}")
132
133 if "alpha" not in args.mode:
134 print(f"dark = {args.dark}")
135 print(f"light = {args.light}")
136 print(f"T = {args.T} K")
137 print(f"sweep_dark = {args.sweep_dark}")
138 print(f"sweep_light = {args.sweep_light}")
139 if args.F0 is not None:
140 print(f"F0 = {args.F0} cm^-2 s^-1")
141 if args.P0 != "":
142 print(f"P0 = {args.P0} W/cm^2")
143# print(f"eq = {args.eq}")
144 print(f"E = {args.E} eV")
145 print(f"lambda_nm = {args.lambda_nm} nm")
146 print(f"S = {args.S} cm^2")
147 print()
148
149 if args.mode == 'analyze':
150 if args.P0 == "":
151 if args.E is not None and args.E == 0.0:
152 if args.lambda_nm is None or args.lambda_nm <= 0.0:
153 raise ValueError(f"lambda_nm must be > 0 nm, got {args.lambda_nm}")
154 elif args.E is not None and args.E < 0.0:
155 raise ValueError(f"E must be >= 0 eV, got {args.E}")
156
157 if args.F0 <= 0.0:
158 raise ValueError(f"F0 must be > 0 cm^-2 s^-1, got {args.F0}")
159
160 args.E = compute_photon_energy_eV(args.E, args.lambda_nm)
161 print(f"E_use = {args.E} eV")
162 args.P0 = compute_p0(args.F0, args.P0, args.E, args.lambda_nm)
163 if args.P0 is None:
164 raise ValueError(f"Got None from compute_p0() using F0={args.F0} and E={args.E}")
165 print(f"P0_use = {args.P0} W/cm^2" if args.P0 is not None else "P0_use = None")
166 else:
167 args.P0 = float(args.P0)
168 if args.P0 <= 0.0:
169 raise ValueError(f"P0 must be > 0 W/cm^2, got {args.P0}")
170
171 if args.T <= 0.0:
172 raise ValueError(f"T must be > 0 K, got {args.T}")
173 if args.S <= 0.0:
174 raise ValueError(f"S must be > 0 cm^2, got {args.S}")
175
176 if args.d <= 0.1:
177 raise ValueError(f"d must be > 0.1 nm, got {args.d}")
178
179
180def read_data(infile, xmin=None, xmax=None, ndataskip=0):
181 """
182 概要:
183 指定されたCSVファイルからI-Vデータを読み込みます。
184 詳細説明:
185 ファイル内のメタデータ（記録時間、データ名）を抽出し、電圧と電流のデータポイントをパースします。
186 複数のスイープを検出し、それぞれをリストに分割します。
187 引数:
188 :param infile: 読み込むCSVファイルのパス。
189 :type infile: str
190 :param xmin: X軸（電圧）の最小値。これより小さい値はスキップされます。
191 :type xmin: float or None
192 :param xmax: X軸（電圧）の最大値。これより大きい値はスキップされます。
193 :type xmax: float or None
194 :param ndataskip: データポイントをスキップする間隔。0の場合、スキップしません。
195 :type ndataskip: int
196 戻り値:
197 :returns: 各スイープのX軸（電圧）データのリスト、各スイープのY軸（電流）データのリスト、ファイルから抽出されたメタデータのタプル。
198 :rtype: tuple[list[numpy.ndarray], list[numpy.ndarray], dict]
199 例外:
200 :raises ValueError: ファイルから有効なI-Vデータが見つからない場合。
201 """
202 print(f"[I/O] Read IV data: {infile}")
203 raw_x = []
204 raw_y = []
205 inf = {"FileName": infile, "RecordTime": "Unknown", "DataName": ("V", "I")}
206
207 icount = 0
208 data_started = False
209 with open(infile, "r", encoding="utf-8-sig", newline="") as f:
210 reader = csv.reader(f)
211 for row in reader:
212 if not row:
213 continue
214 cols = [c.strip() for c in row]
215
216 if len(cols) >= 3 and cols[0] == "MetaData" and "RecordTime" in cols[1]:
217 inf["RecordTime"] = cols[2]
218
219 if cols[0] == "DataName":
220 if len(cols) >= 3:
221 inf["DataName"] = (cols[1], cols[2])
222 data_started = True
223 continue
224
225 if data_started and cols[0] == "DataValue":
226 if len(cols) < 3:
227 continue
228 icount += 1
229 if ndataskip > 0 and (icount - 1) % ndataskip != 0:
230 continue
231 try:
232 x = float(cols[1])
233 y = float(cols[2])
234 except ValueError:
235 continue
236 if xmin is not None and x < xmin:
237 continue
238 if xmax is not None and x > xmax:
239 continue
240 raw_x.append(x)
241 raw_y.append(y)
242
243 if not raw_x:
244 raise ValueError(f"No valid IV data found in {infile}")
245
246 xs_list = []
247 ys_list = []
248
249 cur_x = [raw_x[0]]
250 cur_y = [raw_y[0]]
251 prev_sign = 0
252
253 for i in range(1, len(raw_x)):
254 dx = raw_x[i] - raw_x[i - 1]
255 sign = 0 if abs(dx) < 1e-15 else (1 if dx > 0 else -1)
256
257 if prev_sign == 0:
258 prev_sign = sign
259
260 if sign != 0 and prev_sign != 0 and sign != prev_sign:
261 xs_list.append(np.asarray(cur_x, dtype=float))
262 ys_list.append(np.asarray(cur_y, dtype=float))
263 cur_x = [raw_x[i - 1], raw_x[i]]
264 cur_y = [raw_y[i - 1], raw_y[i]]
265 prev_sign = sign
266 else:
267 cur_x.append(raw_x[i])
268 cur_y.append(raw_y[i])
269 if sign != 0:
270 prev_sign = sign
271
272 if cur_x:
273 xs_list.append(np.asarray(cur_x, dtype=float))
274 ys_list.append(np.asarray(cur_y, dtype=float))
275
276 inf["Points"] = sum(len(x) for x in xs_list)
277 inf["NSweeps"] = len(xs_list)
278 return xs_list, ys_list, inf
279
280
281def read_alpha_from_excel(infile):
282 """
283 概要:
284 Excelファイルから吸収スペクトル（alpha）データを読み込みます。
285 詳細説明:
286 save_alpha_to_excel関数によって保存された形式のExcelファイルを想定しています。
287 alpha_spectrumシートから波長、R、Tr、A、alphaのデータを抽出します。
288 引数:
289 :param infile: 読み込むExcelファイルのパス。
290 :type infile: str
291 戻り値:
292 :returns: 波長、反射率、透過率、吸収率、吸収係数を含む辞書。
293 :rtype: dict
294 例外:
295 :raises ValueError: 指定されたExcelファイルに'alpha_spectrum'シートが見つからない場合。
296 :raises ValueError: ファイルから数値の吸収スペクトルデータが見つからない場合。
297 """
298 print(f"[I/O] Read alpha Excel: {infile}")
299 wb = load_workbook(infile, data_only=True)
300 if "alpha_spectrum" not in wb.sheetnames:
301 raise ValueError(f"Sheet 'alpha_spectrum' not found in {infile}")
302 ws = wb["alpha_spectrum"]
303
304 wl = []
305 R = []
306 Tr = []
307 A = []
308 alpha = []
309
310 first = True
311 for row in ws.iter_rows(values_only=True):
312 if first:
313 first = False
314 continue
315 if row is None or len(row) < 6:
316 continue
317 try:
318 # columns: photon_energy_eV, wavelength_nm, R, Tr, A, alpha_cm^-1
319 wl.append(float(row[1]))
320 R.append(float(row[2]))
321 Tr.append(float(row[3]))
322 A.append(float(row[4]))
323 alpha.append(float(row[5]))
324 except (TypeError, ValueError):
325 continue
326
327 if not wl:
328 raise ValueError(f"No numeric alpha data found in {infile}")
329
330 return {
331 "wl_nm": np.asarray(wl, dtype=float),
332 "R": np.asarray(R, dtype=float),
333 "Tr": np.asarray(Tr, dtype=float),
334 "A": np.asarray(A, dtype=float),
335 "alpha_cm^-1": np.asarray(alpha, dtype=float),
336 }
337
338
339def read_optical_spectrum(infile):
340 """
341 概要:
342 指定されたテキストファイルから光学スペクトルデータ（反射率Rまたは透過率T）を読み込みます。
343 詳細説明:
344 ファイル内のタブ区切りまたはスペース区切りのデータをパースし、
345 波長とスペクトル値のペアを抽出します。
346 重複する波長を処理し、波長順にソートします。
347 引数:
348 :param infile: 読み込む光学スペクトルファイルのパス。
349 :type infile: str
350 戻り値:
351 :returns: 波長 (nm)のnumpy.ndarray、スペクトル値 (RまたはTのパーセンテージ)のnumpy.ndarray、
352 ファイルから抽出されたメタデータを含む辞書のタプル。
353 :rtype: tuple[numpy.ndarray, numpy.ndarray, dict]
354 例外:
355 :raises RuntimeError: ファイルの読み込みに失敗した場合。
356 :raises ValueError: ファイルから数値のスペクトルデータが見つからない場合。
357 """
358
359 print(f"[I/O] Read optical spectrum: {infile}")
360 info = {"FileName": infile, "Label": Path(infile).stem, "YUnit": "a.u."}
361 lines = []
362
363 encodings = ["cp932", "utf-8-sig", "utf-8", "latin1"]
364 last_err = None
365 for enc in encodings:
366 try:
367 with open(infile, "r", encoding=enc, errors="replace") as f:
368 lines = f.readlines()
369 break
370 except Exception as e:
371 last_err = e
372 if not lines:
373 raise RuntimeError(f"Failed to read optical spectrum: {infile} ({last_err})")
374
375 wl = []
376 val = []
377
378 for line in lines:
379 s = line.strip()
380 if not s:
381 continue
382
383 if "\t" in s:
384 parts = [p.strip() for p in s.split("\t")]
385 if len(parts) >= 2:
386 if "ｻﾝﾌﾟﾙ" in parts[0] or "サンプル" in parts[0]:
387 info["Label"] = parts[1]
388 if "%T" in parts[1]:
389 info["YUnit"] = "%T"
390 elif "%R" in parts[1]:
391 info["YUnit"] = "%R"
392
393 tokens = s.replace(",", " ").replace("\t", " ").split()
394 nums = []
395 for t in tokens:
396 try:
397 nums.append(float(t))
398 except ValueError:
399 pass
400 if len(nums) >= 2:
401 wl.append(nums[0])
402 val.append(nums[1])
403
404 if not wl:
405 raise ValueError(f"No numeric spectrum data found in {infile}")
406
407 wl = np.asarray(wl, dtype=float)
408 val = np.asarray(val, dtype=float)
409
410 uu, idx = np.unique(wl, return_index=True)
411 wl = wl[np.sort(idx)]
412 val = val[np.sort(idx)]
413
414 order = np.argsort(wl)
415 wl = wl[order]
416 val = val[order]
417 return wl, val, info
418
419
420def choose_sweep(xs_list, ys_list, sweep_index=0):
421 """
422 概要:
423 複数のスイープデータから指定されたインデックスの単一スイープを選択します。
424 引数:
425 :param xs_list: 各スイープのX軸（電圧）データのリスト。
426 :type xs_list: list[numpy.ndarray]
427 :param ys_list: 各スイープのY軸（電流）データのリスト。
428 :type ys_list: list[numpy.ndarray]
429 :param sweep_index: 選択するスイープのインデックス。
430 :type sweep_index: int
431 戻り値:
432 :returns: 選択されたスイープのX軸データとY軸データ。
433 :rtype: tuple[numpy.ndarray, numpy.ndarray]
434 例外:
435 :raises ValueError: スイープデータが利用できない場合。
436 :raises IndexError: sweep_indexが範囲外の場合。
437 """
438 if len(xs_list) == 0:
439 raise ValueError("No sweep data available.")
440 idx = int(sweep_index)
441 if idx < 0 or idx >= len(xs_list):
442 raise IndexError(f"sweep_index={idx} is out of range (0..{len(xs_list)-1})")
443 return np.asarray(xs_list[idx], dtype=float), np.asarray(ys_list[idx], dtype=float)
444
445
446def consolidate_duplicate_x(x, y):
447 """
448 概要:
449 X軸に重複する値がある場合、Y軸の対応する値を平均して重複を解消します。
450 引数:
451 :param x: X軸データ。
452 :type x: numpy.ndarray
453 :param y: Y軸データ。
454 :type y: numpy.ndarray
455 戻り値:
456 :returns: 重複が解消されたX軸データと、重複が解消され平均化されたY軸データ。
457 :rtype: tuple[numpy.ndarray, numpy.ndarray]
458 """
459 x = np.asarray(x, dtype=float)
460 y = np.asarray(y, dtype=float)
461 order = np.argsort(x)
462 xs = x[order]
463 ys = y[order]
464
465 ux = []
466 uy = []
467 i = 0
468 n = len(xs)
469 while i < n:
470 xv = xs[i]
471 vals = [ys[i]]
472 j = i + 1
473 while j < n and abs(xs[j] - xv) < 1e-15:
474 vals.append(ys[j])
475 j += 1
476 ux.append(xv)
477 uy.append(float(np.mean(vals)))
478 i = j
479 return np.asarray(ux), np.asarray(uy)
480
481
482def smooth_polyfit(y, window_points=5, poly_order=3):
483 """
484 概要:
485 多項式フィッティングを用いたSavitzky-Golay風の平滑化をデータに適用します。
486 詳細説明:
487 各データポイントを中心に指定されたウィンドウ内のデータに対して多項式フィッティングを行い、
488 中心点の値を予測することで平滑化を行います。
489 引数:
490 :param y: 平滑化するY軸データ。
491 :type y: numpy.ndarray
492 :param window_points: フィッティングに使用するウィンドウ内のデータポイント数。奇数である必要があります。
493 :type window_points: int
494 :param poly_order: 多項式フィッティングの次数。
495 :type poly_order: int
496 戻り値:
497 :returns: 平滑化されたY軸データ。
498 :rtype: numpy.ndarray
499 """
500 y = np.asarray(y, dtype=float)
501 n = len(y)
502 if n < 3:
503 return y.copy()
504
505 if window_points < 3:
506 window_points = 3
507 if window_points % 2 == 0:
508 window_points += 1
509 if window_points > n:
510 window_points = n if (n % 2 == 1) else n - 1
511 if poly_order >= window_points:
512 poly_order = window_points - 1
513 poly_order = max(poly_order, 1)
514
515 half = window_points // 2
516 ys = np.empty(n, dtype=float)
517
518 for i in range(n):
519 i0 = max(0, i - half)
520 i1 = min(n, i + half + 1)
521 while (i1 - i0) < window_points:
522 if i0 > 0:
523 i0 -= 1
524 elif i1 < n:
525 i1 += 1
526 else:
527 break
528 idx = np.arange(i0, i1, dtype=float)
529 yy = y[i0:i1]
530 xloc = idx - i
531 deg = min(poly_order, len(yy) - 1)
532 coeff = np.polyfit(xloc, yy, deg)
533 ys[i] = np.polyval(coeff, 0.0)
534 return ys
535
536
537def local_poly_value(x, y, x0, npts=7, order=3):
538 """
539 概要:
540 指定されたX座標の周囲のデータポイントを使用して、局所的な多項式フィッティングを行い、
541 x0におけるY値を推定します。
542 引数:
543 :param x: X軸データ。
544 :type x: numpy.ndarray
545 :param y: Y軸データ。
546 :type y: numpy.ndarray
547 :param x0: Y値を推定するX座標。
548 :type x0: float
549 :param npts: フィッティングに使用するデータポイント数。
550 :type npts: int
551 :param order: 多項式フィッティングの次数。
552 :type order: int
553 戻り値:
554 :returns: x0における推定されたY値。
555 :rtype: float
556 """
557 x = np.asarray(x, dtype=float)
558 y = np.asarray(y, dtype=float)
559 idx = np.argsort(np.abs(x - x0))[:max(2, npts)]
560 xs = x[idx]
561 ys = y[idx]
562 sidx = np.argsort(xs)
563 xs = xs[sidx]
564 ys = ys[sidx]
565 deg = min(order, len(xs) - 1)
566 coeff = np.polyfit(xs, ys, deg)
567 return float(np.polyval(coeff, x0))
568
569
570def zero_crossing_x(x, y):
571 """
572 概要:
573 Y値がゼロを横切るX座標を線形補間によって見つけます。
574 詳細説明:
575 2つの連続するデータポイントの間でY値の符号が変わる点、
576 またはY値が最もゼロに近い点のX座標を返します。
577 引数:
578 :param x: X軸データ。
579 :type x: numpy.ndarray
580 :param y: Y軸データ。
581 :type y: numpy.ndarray
582 戻り値:
583 :returns: Y値がゼロを横切るX座標。
584 :rtype: float
585 """
586 x = np.asarray(x, dtype=float)
587 y = np.asarray(y, dtype=float)
588 order = np.argsort(x)
589 x = x[order]
590 y = y[order]
591
592 for i in range(len(x) - 1):
593 y1, y2 = y[i], y[i + 1]
594 if y1 == 0:
595 return float(x[i])
596 if y1 * y2 < 0:
597 return float(x[i] + (0 - y1) * (x[i + 1] - x[i]) / (y2 - y1))
598 return float(x[np.argmin(np.abs(y))])
599
600
601def interpolate_to_common_wavelength(wl_ref, y_ref, wl_target):
602 """
603 概要:
604 参照波長スケール上のデータをターゲット波長スケールに線形補間します。
605 引数:
606 :param wl_ref: 参照波長データ (nm)。
607 :type wl_ref: numpy.ndarray
608 :param y_ref: 参照Y軸データ。
609 :type y_ref: numpy.ndarray
610 :param wl_target: 補間対象のターゲット波長データ (nm)。
611 :type wl_target: numpy.ndarray
612 戻り値:
613 :returns: ターゲット波長スケールに補間されたY軸データ。
614 :rtype: numpy.ndarray
615 """
616 wl_ref = np.asarray(wl_ref, dtype=float)
617 y_ref = np.asarray(y_ref, dtype=float)
618 wl_target = np.asarray(wl_target, dtype=float)
619 return np.interp(wl_target, wl_ref, y_ref, left=np.nan, right=np.nan)
620
621
622def compute_photon_energy_eV(E_eV, lambda_nm):
623 """
624 概要:
625 フォトンエネルギー（eV）を計算します。
626 詳細説明:
627 フォトンエネルギーが直接与えられている場合はそれを使用し、
628 そうでない場合は波長から計算します。
629 引数:
630 :param E_eV: フォトンエネルギー (eV)。0より大きい場合はこれを使用。
631 :type E_eV: float
632 :param lambda_nm: フォトン波長 (nm)。E_eVが0の場合にこれを使用。
633 :type lambda_nm: float
634 戻り値:
635 :returns: 計算されたフォトンエネルギー (eV)。
636 :rtype: float
637 """
638 if E_eV is not None and E_eV > 0:
639 return float(E_eV)
640 return float(1239.841984 / lambda_nm)
641
642
643def compute_p0(F0, P0, E_eV, lambda_nm):
644 """
645 概要:
646 入射フォトンエネルギー密度（P0）を計算します。
647 詳細説明:
648 入射フォトンフラックス（F0）が与えられている場合はそれから計算し、
649 そうでない場合は直接指定されたP0を使用します。
650 引数:
651 :param F0: 入射フォトンフラックス (cm^-2 s^-1)。Noneまたは0の場合はP0を使用します。
652 :type F0: float or None
653 :param P0: 入射フォトンエネルギー密度 (W/cm^2)。Noneの場合はF0から計算します。
654 :type P0: str or None
655 :param E_eV: フォトンエネルギー (eV)。F0からP0を計算する場合に使用。
656 :type E_eV: float
657 :param lambda_nm: フォトン波長 (nm)。F0からP0を計算する場合に使用。
658 :type lambda_nm: float
659 戻り値:
660 :returns: 計算された入射フォトンエネルギー密度 (W/cm^2)。
661 :rtype: float or None
662 """
663 if F0 is not None and F0 > 0:
664 E_use = compute_photon_energy_eV(E_eV, lambda_nm)
665 return float(F0 * E_use * E_CHARGE)
666 return None if P0 is None else float(P0)
667
668
669def pv_metrics_from_iv(V, I, S):
670 """
671 概要:
672 I-Vデータから太陽電池の主要な性能指標（Voc, Jsc, FF, Pmaxなど）を計算します。
673 引数:
674 :param V: 電圧データ (V)。
675 :type V: numpy.ndarray
676 :param I: 電流データ (A)。
677 :type I: numpy.ndarray
678 :param S: 電極面積 (cm^2)。
679 :type S: float
680 戻り値:
681 :returns: 以下の主要な太陽電池性能指標を含む辞書。
682 :rtype: dict
683 - Voc_V (float): 開放電圧 (V)。
684 - Jsc_A_cm2 (float): 短絡電流密度 (A/cm^2)。
685 - Jsc_mA_cm2 (float): 短絡電流密度 (mA/cm^2)。
686 - Vop_V (float): 最大出力動作電圧 (V)。
687 - Jop_A_cm2 (float): 最大出力動作電流密度 (A/cm^2)。
688 - Jop_mA_cm2 (float): 最大出力動作電流密度 (mA/cm^2)。
689 - Pmax_W_cm2 (float): 最大出力電力密度 (W/cm^2)。
690 - Pmax_mW_cm2 (float): 最大出力電力密度 (mW/cm^2)。
691 - FF (float): 曲線因子。
692 """
693 V = np.asarray(V, dtype=float)
694 I = np.asarray(I, dtype=float)
695 J = I / S
696
697 Jsc = local_poly_value(V, J, 0.0, npts=7, order=3)
698 Voc = zero_crossing_x(V, I)
699
700 Pgen = -V * J
701 idx = int(np.argmax(Pgen))
702 Vop = float(V[idx])
703 Jop = float(J[idx])
704 Pmax = float(Pgen[idx])
705
706 denom = abs(Voc * Jsc)
707 FF = float(Pmax / denom) if denom > 1e-30 else float("nan")
708
709 return {
710 "Voc_V": float(Voc),
711 "Jsc_A_cm2": float(Jsc),
712 "Jsc_mA_cm2": float(Jsc * 1e3),
713 "Vop_V": float(Vop),
714 "Jop_A_cm2": float(Jop),
715 "Jop_mA_cm2": float(Jop * 1e3),
716 "Pmax_W_cm2": float(Pmax),
717 "Pmax_mW_cm2": float(Pmax * 1e3),
718 "FF": float(FF),
719 }
720
721
722def fit_local_line(x, y, center_idx, npts=7):
723 """
724 概要:
725 指定された中心インデックスの周囲のデータポイントを用いて局所的な線形フィッティングを行います。
726 引数:
727 :param x: X軸データ。
728 :type x: numpy.ndarray
729 :param y: Y軸データ。
730 :type y: numpy.ndarray
731 :param center_idx: フィッティングの中心となるデータポイントのインデックス。
732 :type center_idx: int
733 :param npts: フィッティングに使用するデータポイント数。
734 :type npts: int
735 戻り値:
736 :returns: フィットされた直線の傾き、Y切片、フィッティングに使用されたXデータ、
737 Yデータ、開始インデックス、終了インデックス（排他的）のタプル。
738 :rtype: tuple[float, float, numpy.ndarray, numpy.ndarray, int, int]
739 """
740 x = np.asarray(x, dtype=float)
741 y = np.asarray(y, dtype=float)
742 n = len(x)
743 npts = max(2, min(int(npts), n))
744 half = npts // 2
745 i0 = max(0, center_idx - half)
746 i1 = min(n, i0 + npts)
747 i0 = max(0, i1 - npts)
748 xx = x[i0:i1]
749 yy = y[i0:i1]
750 a, b = np.polyfit(xx, yy, 1)
751 return float(a), float(b), xx, yy, i0, i1
752
753
754def estimate_rs_tangent_point(V, I):
755 """
756 概要:
757 I-V曲線から直列抵抗（Rs）を推定するための接点と抵抗値を計算します。
758 詳細説明:
759 順方向バイアス領域でdI/dVが最大となる点を特定し、その点での接線からRsを計算します。
760 引数:
761 :param V: 電圧データ (V)。
762 :type V: numpy.ndarray
763 :param I: 電流データ (A)。
764 :type I: numpy.ndarray
765 戻り値:
766 :returns: 直列抵抗推定に関する情報を含む辞書。
767 :rtype: dict
768 - v_rep (float): 接点電圧 (V)。
769 - i_rep (float): 接点電流 (A)。
770 - Is (float): 接点電流 (A) (互換性のため)。
771 - Rs (float): 推定された直列抵抗 (Ω)。
772 - slope (float): 接線の傾き (A/V)。
773 - intercept (float): 接線のY切片 (A)。
774 - xx (numpy.ndarray): 接線フィッティングに使用されたXデータ。
775 - yy (numpy.ndarray): 接線フィッティングに使用されたYデータ。
776 - idx (int): 接点のインデックス。
777 """
778 xu, yu = consolidate_duplicate_x(V, I)
779 y_sm = smooth_polyfit(yu, window_points=5, poly_order=3)
780 dydv = np.gradient(y_sm, xu)
781
782 pos_idx = np.where(xu > 0.0)[0]
783 if len(pos_idx) == 0:
784 idx_pos = int(np.argmax(dydv))
785 else:
786 pad = min(2, max(0, len(pos_idx) // 4))
787 cand = pos_idx[pad:len(pos_idx)-pad] if len(pos_idx) - 2 * pad >= 1 else pos_idx
788 idx_local = int(np.argmax(dydv[cand]))
789 idx_pos = int(cand[idx_local])
790
791 a, b, xx, yy, i0, i1 = fit_local_line(xu, y_sm, idx_pos, npts=7)
792 v_rep = float(xu[idx_pos])
793 i_rep = float(a * v_rep + b)
794 rs = float("inf") if abs(a) < 1e-30 else float(1.0 / a)
795 return {"v_rep": v_rep, "i_rep": i_rep, "Is": i_rep, "Rs": rs,
796 "slope": float(a), "intercept": float(b), "xx": xx, "yy": yy, "idx": idx_pos}
797
798
799def estimate_rsh_tangent_point(V, I):
800 """
801 概要:
802 I-V曲線から並列抵抗（Rsh）を推定するための接点と抵抗値を計算します。
803 詳細説明:
804 逆方向バイアス領域でdI/dVが最小となる点（絶対値）を特定し、
805 その点での接線からRshを計算します。
806 引数:
807 :param V: 電圧データ (V)。
808 :type V: numpy.ndarray
809 :param I: 電流データ (A)。
810 :type I: numpy.ndarray
811 戻り値:
812 :returns: 並列抵抗推定に関する情報を含む辞書。
813 :rtype: dict
814 - v_rep (float): 接点電圧 (V)。
815 - i_rep (float): 接点電流 (A)。
816 - Ish (float): 接点電流 (A) (互換性のため)。
817 - Rsh (float): 推定された並列抵抗 (Ω)。
818 - slope (float): 接線の傾きの絶対値 (A/V)。
819 - intercept (float): 接線のY切片 (A)。
820 - slope_fit_raw (float): フィットされた直線の元の傾き (A/V)。
821 - intercept_fit_raw (float): フィットされた直線の元のY切片 (A)。
822 - xx (numpy.ndarray): 接線フィッティングに使用されたXデータ。
823 - yy (numpy.ndarray): 接線フィッティングに使用されたYデータ。
824 - idx (int): 接点のインデックス。
825 """
826 xu, yu = consolidate_duplicate_x(V, I)
827 y_sm = smooth_polyfit(yu, window_points=5, poly_order=3)
828 dydv = np.gradient(y_sm, xu)
829
830 neg_mask = xu < 0.0
831 if np.any(neg_mask):
832 dneg = np.abs(dydv[neg_mask])
833 idx_neg_local = int(np.argmin(dneg))
834 idx_neg = np.where(neg_mask)[0][idx_neg_local]
835 else:
836 idx_neg = int(np.argmin(np.abs(dydv)))
837
838 a_fit, b_fit, xx, yy, i0, i1 = fit_local_line(xu, y_sm, idx_neg, npts=7)
839 v_rep = float(xu[idx_neg])
840 i_fit_ref = float(a_fit * v_rep + b_fit)
841 a = abs(float(a_fit))
842 b = float(i_fit_ref - a * v_rep)
843 i_rep = float(a * v_rep + b)
844 rsh = float("inf") if abs(a) < 1e-30 else float(1.0 / a)
845 return {"v_rep": v_rep, "i_rep": i_rep, "Ish": i_rep, "Rsh": rsh,
846 "slope": float(a), "intercept": float(b), "slope_fit_raw": float(a_fit),
847 "intercept_fit_raw": float(b_fit), "xx": xx, "yy": yy, "idx": idx_neg}
848
849
850def estimate_ndiode_representative_point(V, I, T=300.0, Ish=None):
851 """
852 概要:
853 I-V曲線からダイオード因子（ndiode）を推定するための代表点と値を計算します。
854 詳細説明:
855 順方向バイアス領域でlog(abs(I))の二階微分が最小となる点を特定し、
856 その点での一次微分からダイオード因子を計算します。
857 引数:
858 :param V: 電圧データ (V)。
859 :type V: numpy.ndarray
860 :param I: 電流データ (A)。
861 :type I: numpy.ndarray
862 :param T: 温度 (K)。
863 :type T: float
864 :param Ish: 並列抵抗に流れる電流の推定値。これを考慮してIを調整します。
865 :type Ish: float or None
866 戻り値:
867 :returns: ダイオード因子推定に関する情報を含む辞書。計算できなかった場合はNone。
868 :rtype: dict or None
869 - v_rep (float): 代表点電圧 (V)。
870 - ndiode (float): 推定されたダイオード因子。
871 - slope_logI (float): 代表点でのlog(abs(I))の傾き (1/V)。
872 - curvature_logI (float): 代表点でのlog(abs(I))の二階微分 (1/V^2)。
873 """
874 xu, yu = consolidate_duplicate_x(V, I)
875 order = np.argsort(xu)
876 xu = xu[order]
877 yu = yu[order]
878
879 ish_abs = abs(Ish) if Ish is not None else 0.0
880 mask = (xu > 0.0) & (np.abs(yu) > ish_abs)
881 if np.count_nonzero(mask) < 5:
882 return None
883
884 x = xu[mask]
885 y = np.abs(yu[mask]) + EPS_I
886
887 logI = np.log(y)
888 logI_sm = smooth_polyfit(logI, window_points=5, poly_order=3)
889 d1 = np.gradient(logI_sm, x)
890 d2 = np.gradient(d1, x)
891
892 good = np.isfinite(d1) & np.isfinite(d2) & (d1 > 0.0)
893 if np.count_nonzero(good) < 3:
894 return None
895
896 xg = x[good]
897 d1g = d1[good]
898 d2g = d2[good]
899
900 idx = int(np.argmin(np.abs(d2g)))
901 v_rep = float(xg[idx])
902 slope = float(d1g[idx])
903 ndiode = float(E_CHARGE / (KB * T * slope)) if abs(slope) > 1e-30 else float("nan")
904 return {"v_rep": v_rep, "ndiode": ndiode, "slope_logI": slope, "curvature_logI": float(d2g[idx])}
905
906
907def estimate_initial_params(V, I_meas, T=300.0):
908 """
909 概要:
910 I-Vデータから単一ダイオードモデルの初期パラメータ（I0, ndiode, IPV, Rs, Rsh）を推定します。
911 詳細説明:
912 estimate_rs_tangent_point、estimate_rsh_tangent_point、
913 estimate_ndiode_representative_pointなどの補助関数を用いて、
914 各パラメータの初期値を経験的に決定します。
915 引数:
916 :param V: 電圧データ (V)。
917 :type V: numpy.ndarray
918 :param I_meas: 測定電流データ (A)。
919 :type I_meas: numpy.ndarray
920 :param T: 温度 (K)。
921 :type T: float
922 戻り値:
923 :returns: 推定された初期パラメータを含む辞書。
924 :rtype: dict
925 - I0 (float): 逆方向飽和電流 (A)。
926 - ndiode (float): ダイオード因子。
927 - IPV (float): 光電流 (A)。
928 - Rs (float): 直列抵抗 (Ω)。
929 - Rsh (float): 並列抵抗 (Ω)。
930 - Ish (float): 並列抵抗に流れる電流 (A)。
931 - Vsh (float): Rsh推定における代表点電圧 (V)。
932 - Vnd (float): ndiode推定における代表点電圧 (V)。
933 """
934 xu, yu = consolidate_duplicate_x(V, I_meas)
935 order = np.argsort(xu)
936 xu = xu[order]
937 yu = yu[order]
938 y_sm = smooth_polyfit(yu, window_points=5, poly_order=3)
939
940 rs_info = estimate_rs_tangent_point(V, I_meas)
941 Rs = abs(float(rs_info["Rs"]))
942
943 rsh_info = estimate_rsh_tangent_point(V, I_meas)
944 Rsh = abs(float(rsh_info["Rsh"]))
945 Ish = float(rsh_info["Ish"])
946
947 I0 = abs(Ish)
948 IPV = -local_poly_value(xu, y_sm, 0.0, npts=7, order=3)
949
950 nd_info = estimate_ndiode_representative_point(V, I_meas, T=T, Ish=Ish)
951 if nd_info is None or (not np.isfinite(nd_info["ndiode"])):
952 print()
953 print("########################################################################")
954 print(" Warning!!!: Could not get valid ndiode estimation")
955 print(" Use ndiode = 1000.0 instead but don't refer to this value")
956 print("########################################################################")
957 ndiode = 1000.0
958 Vnd = float("nan")
959 else:
960 ndiode = float(nd_info["ndiode"])
961 Vnd = float(nd_info["v_rep"])
962
963 return {
964 "I0": I0,
965 "ndiode": ndiode,
966 "IPV": IPV,
967 "Rs": Rs,
968 "Rsh": Rsh,
969 "Ish": Ish,
970 "Vsh": float(rsh_info["v_rep"]),
971 "Vnd": Vnd,
972 }
973
974
975def analyze_optical(wl_R_nm, R_percent, wl_T_nm, T_percent, d_nm):
976 """
977 概要:
978 反射率（R）と透過率（T）のスペクトルデータから吸収率（A）と吸収係数（α）を計算します。
979 詳細説明:
980 異なる波長スケールのRとTを共通の波長スケールに補間し、
981 RとTのデータに基づいてAとαを計算します。
982 引数:
983 :param wl_R_nm: 反射率スペクトルの波長データ (nm)。
984 :type wl_R_nm: numpy.ndarray
985 :param R_percent: 反射率データ (%)。
986 :type R_percent: numpy.ndarray
987 :param wl_T_nm: 透過率スペクトルの波長データ (nm)。
988 :type wl_T_nm: numpy.ndarray
989 :param T_percent: 透過率データ (%)。
990 :type T_percent: numpy.ndarray
991 :param d_nm: PV層の膜厚 (nm)。
992 :type d_nm: float
993 戻り値:
994 :returns: 計算された光学特性値を含む辞書。
995 :rtype: dict
996 - wl_nm (numpy.ndarray): 共通波長データ (nm)。
997 - R (numpy.ndarray): 反射率 (0-1)。
998 - Tr (numpy.ndarray): 透過率 (0-1)。
999 - A (numpy.ndarray): 吸収率 (0-1)。
1000 - alpha_cm^-1 (numpy.ndarray): 吸収係数 (cm^-1)。
1001 """
1002 wl_common = np.unique(np.concatenate([wl_R_nm, wl_T_nm]))
1003 wl_common.sort()
1004
1005 R = interpolate_to_common_wavelength(wl_R_nm, R_percent, wl_common) / 100.0
1006 Tr = interpolate_to_common_wavelength(wl_T_nm, T_percent, wl_common) / 100.0
1007
1008 mask = np.isfinite(R) & np.isfinite(Tr)
1009 wl = wl_common[mask]
1010 R = np.clip(R[mask], 0.0, 1.0)
1011 Tr = np.clip(Tr[mask], 0.0, 1.0)
1012
1013 A = np.clip(1.0 - R - Tr, 0.0, 1.0)
1014
1015 d_cm = d_nm * 1e-7
1016 denom = np.clip(1.0 - R, 1e-12, None)
1017 alpha_cm = -np.log(np.clip(Tr / denom, 1e-12, None)) / max(d_cm, 1e-30)
1018
1019 return {"wl_nm": wl, "R": R, "Tr": Tr, "A": A, "alpha_cm^-1": alpha_cm}
1020
1021
1022def quantum_efficiencies(F0, A_abs, JPV_A_cm2, Jsc_A_cm2):
1023 """
1024 概要:
1025 量子効率（EQEとIQE）を計算します。
1026 詳細説明:
1027 入射フォトンフラックス、吸収率、光電流（JPV）、短絡電流密度（Jsc）に基づいて、
1028 生成効率と収集効率を評価します。
1029 引数:
1030 :param F0: 入射フォトンフラックス (cm^-2 s^-1)。Noneまたは0の場合はNaNを返します。
1031 :type F0: float or None
1032 :param A_abs: 特定のエネルギーにおける吸収率。NoneまたはNaNの場合はIQEをNaNとします。
1033 :type A_abs: float or None
1034 :param JPV_A_cm2: 光電流密度 (A/cm^2)。
1035 :type JPV_A_cm2: float
1036 :param Jsc_A_cm2: 短絡電流密度 (A/cm^2)。
1037 :type Jsc_A_cm2: float
1038 戻り値:
1039 :returns: 計算された量子効率を含む辞書。
1040 :rtype: dict
1041 - EQE_gen (float): 生成外部量子効率。
1042 - IQE_gen (float): 生成内部量子効率。
1043 - EQE (float): 外部量子効率。
1044 - IQE (float): 内部量子効率。
1045 """
1046 out = {"EQE_gen": float("nan"), "IQE_gen": float("nan"), "EQE": float("nan"), "IQE": float("nan")}
1047 if F0 is None or F0 <= 0:
1048 return out
1049
1050 denom_ext = E_CHARGE * F0
1051 out["EQE_gen"] = abs(JPV_A_cm2) / denom_ext
1052 out["EQE"] = abs(Jsc_A_cm2) / denom_ext
1053
1054 if A_abs is not None and np.isfinite(A_abs) and A_abs > 0:
1055 denom_int = E_CHARGE * F0 * A_abs
1056 out["IQE_gen"] = abs(JPV_A_cm2) / denom_int
1057 out["IQE"] = abs(Jsc_A_cm2) / denom_int
1058 return out
1059
1060
1061def alpha_at_energy(optical, E_eV):
1062 """
1063 概要:
1064 指定されたフォトンエネルギーにおける吸収係数（α）を光学スペクトルデータから補間して取得します。
1065 引数:
1066 :param optical: analyze_opticalまたはread_alpha_from_excelから返された光学データ辞書。
1067 :type optical: dict or None
1068 :param E_eV: 取得したいフォトンエネルギー (eV)。
1069 :type E_eV: float
1070 戻り値:
1071 :returns: 指定されたエネルギーにおける吸収係数 (cm^-1)。データがない場合はNaN。
1072 :rtype: float
1073 例外:
1074 :raises ValueError: 補間された吸収係数がNaNの場合。
1075 """
1076 if optical is None:
1077 return float("nan")
1078
1079 wl = np.asarray(optical["wl_nm"], dtype=float)
1080 alpha = np.asarray(optical["alpha_cm^-1"], dtype=float)
1081 energy = 1239.841984 / wl
1082 order = np.argsort(energy)
1083 energy = energy[order]
1084 alpha = alpha[order]
1085
1086 alpha_interp = np.interp(E_eV, energy, alpha, left=np.nan, right=np.nan)
1087 if not np.isfinite(alpha_interp):
1088 raise ValueError(f"got nan for alpha_at_energy()")
1089
1090 return float(alpha_interp)
1091
1092def absorptance_from_alpha(alpha_cm, d_nm):
1093 """
1094 概要:
1095 吸収係数と膜厚から吸収率を計算します。
1096 詳細説明:
1097 A = 1 - exp(-alpha * d) の式を用いて計算します。
1098 引数:
1099 :param alpha_cm: 吸収係数 (cm^-1)。
1100 :type alpha_cm: float or None
1101 :param d_nm: PV層の膜厚 (nm)。
1102 :type d_nm: float
1103 戻り値:
1104 :returns: 計算された吸収率。
1105 :rtype: float
1106 例外:
1107 :raises ValueError: alpha_cmがNoneまたは無効な値の場合。
1108 """
1109 if alpha_cm is None or not np.isfinite(alpha_cm):
1110 raise ValueError(f"got None or invalid value for alpha_cm")
1111 d_cm = d_nm * 1e-7
1112 return float(1.0 - np.exp(-alpha_cm * d_cm))
1113
1114
1115def save_alpha_to_excel(optical, outfile_xlsx, meta=None):
1116 """
1117 概要:
1118 吸収スペクトルデータをExcelファイルに保存します。
1119 詳細説明:
1120 alpha_spectrumというシートにフォトンエネルギー、波長、R、Tr、A、alphaのデータを書き込みます。
1121 summaryシートには膜厚やファイル名、平均値などのメタデータを保存します。
1122 引数:
1123 :param optical: analyze_optical関数から返された光学データ辞書。
1124 :type optical: dict
1125 :param outfile_xlsx: 保存するExcelファイルのパス。
1126 :type outfile_xlsx: str
1127 :param meta: 保存する追加のメタデータ。
1128 :type meta: dict or None
1129 戻り値:
1130 なし
1131 """
1132 print(f"[I/O] Save alpha Excel: {outfile_xlsx}")
1133 wb = Workbook()
1134 ws = wb.active
1135 ws.title = "alpha_spectrum"
1136
1137 header_fill = PatternFill("solid", fgColor="1F4E78")
1138 header_font = Font(color="FFFFFF", bold=True)
1139
1140 headers = ["photon_energy_eV", "wavelength_nm", "R", "Tr", "A", "alpha_cm^-1"]
1141 ws.append(headers)
1142 for c in ws[1]:
1143 c.fill = header_fill
1144 c.font = header_font
1145
1146 hc_eVnm = 1239.841984
1147 energy_eV = hc_eVnm / optical["wl_nm"]
1148 order = np.argsort(energy_eV)
1149
1150 for i in order:
1151 ws.append([
1152 float(energy_eV[i]),
1153 float(optical["wl_nm"][i]),
1154 float(optical["R"][i]),
1155 float(optical["Tr"][i]),
1156 float(optical["A"][i]),
1157 float(optical["alpha_cm^-1"][i]),
1158 ])
1159
1160 for col in ["A", "B", "C", "D", "E", "F"]:
1161 ws.column_dimensions[col].width = 18
1162
1163 ws2 = wb.create_sheet("summary")
1164 ws2["A1"] = "d_nm"
1165 ws2["B1"] = meta.get("d_nm", "") if meta else ""
1166 ws2["A2"] = "R_file"
1167 ws2["B2"] = meta.get("R", "") if meta else ""
1168 ws2["A3"] = "Tr_file"
1169 ws2["B3"] = meta.get("Tr", "") if meta else ""
1170 ws2["A5"] = "R_mean"
1171 ws2["B5"] = float(np.mean(optical["R"]))
1172 ws2["A6"] = "Tr_mean"
1173 ws2["B6"] = float(np.mean(optical["Tr"]))
1174 ws2["A7"] = "A_mean"
1175 ws2["B7"] = float(np.mean(optical["A"]))
1176 ws2["A8"] = "alpha_mean_cm^-1"
1177 ws2["B8"] = float(np.mean(optical["alpha_cm^-1"]))
1178
1179 for col in ["A", "B"]:
1180 ws2.column_dimensions[col].width = 24
1181
1182 wb.save(outfile_xlsx)
1183 print(f"[I/O] Saved alpha Excel: {outfile_xlsx}")
1184
1185
1186def plot_alpha(optical, outfile=None, pause=False):
1187 """
1188 概要:
1189 吸収係数スペクトル（α vs フォトンエネルギー）をプロットします。
1190 引数:
1191 :param optical: analyze_optical関数から返された光学データ辞書。
1192 :type optical: dict
1193 :param outfile: プロットを保存するファイルパス。Noneの場合、画面に表示します。
1194 :type outfile: str or None
1195 :param pause: プロットを閉じるまで一時停止するかどうか。Trueの場合、plt.show()を呼び出します。
1196 :type pause: bool
1197 戻り値:
1198 :returns: 作成されたMatplotlib Figureオブジェクト。
1199 :rtype: matplotlib.figure.Figure
1200 """
1201 fig, ax = plt.subplots(1, 1, figsize=(6.2, 4.8))
1202 hc_eVnm = 1239.841984
1203 energy_eV = hc_eVnm / optical["wl_nm"]
1204
1205 order = np.argsort(energy_eV)
1206 x = energy_eV[order]
1207 y = optical["alpha_cm^-1"][order]
1208
1209 ax.plot(x, y, "-", linewidth=1.5)
1210 ax.tick_params(axis='x', labelsize=fontsize)
1211 ax.tick_params(axis='y', labelsize=fontsize)
1212 ax.set_xlabel("Photon energy / eV", fontsize=fontsize)
1213 ax.set_ylabel(r"$\alpha$ / cm$^{-1}$", fontsize=fontsize)
1214 ax.set_title(r"Absorption coefficient $\alpha(E)$", fontsize=fontsize)
1215 ax.grid(True)
1216
1217 fig.tight_layout()
1218 if outfile:
1219 print(f"[I/O] Save plot: {outfile}")
1220 fig.savefig(outfile, dpi=160)
1221
1222 if pause:
1223 plt.show()
1224 else:
1225 plt.pause(0.01)
1226 return fig
1227
1228
1229def estimate_ndiode_curve(V, I, T=300.0):
1230 """
1231 概要:
1232 I-Vデータから電圧に対するダイオード因子（ndiode）の曲線を推定します。
1233 詳細説明:
1234 log(abs(I))の一次微分を計算し、それに基づいて各電圧点におけるダイオード因子を導出します。
1235 引数:
1236 :param V: 電圧データ (V)。
1237 :type V: numpy.ndarray
1238 :param I: 電流データ (A)。
1239 :type I: numpy.ndarray
1240 :param T: 温度 (K)。
1241 :type T: float
1242 戻り値:
1243 :returns: ダイオード因子が計算された電圧データと、各電圧に対応する推定されたダイオード因子。
1244 :rtype: tuple[numpy.ndarray, numpy.ndarray]
1245 """
1246 V = np.asarray(V, dtype=float)
1247 I = np.asarray(I, dtype=float)
1248 xu, yu = consolidate_duplicate_x(V, I)
1249 order = np.argsort(xu)
1250 xu = xu[order]
1251 yu = yu[order]
1252
1253 mask = (xu > 0.0) & (yu > 0.0)
1254 if np.count_nonzero(mask) < 5:
1255 return np.array([]), np.array([])
1256
1257 x = xu[mask]
1258 y = yu[mask]
1259 logI = np.log(np.abs(y) + EPS_I)
1260 logI_sm = smooth_polyfit(logI, window_points=5, poly_order=3)
1261 dlogIdV = np.gradient(logI_sm, x)
1262
1263 with np.errstate(divide="ignore", invalid="ignore"):
1264 ndiode = E_CHARGE / (KB * T * dlogIdV)
1265 good = np.isfinite(ndiode) & (ndiode > 0)
1266 return x[good], ndiode[good]
1267
1268
1269def plot_iv_comparison(dark_sweeps, light_sweeps, S, light_metrics,
1270 dark_params=None, light_params=None, outfile=None, T=300.0, pause=False):
1271 """
1272 概要:
1273 暗電流I-Vおよび光照射I-Vデータと、関連する推定値や解析結果を比較プロットします。
1274 詳細説明:
1275 以下の4つのサブプロットを作成します。
1276 - log_10(abs(I))-V曲線（暗電流と光電流）。
1277 - 線形J-V曲線（順方向と逆方向）。RsとRshの接線、Ish(0)点を表示。
1278 - ndiode-V曲線と代表点。
1279 - 光起電力出力（-J vs V）。MPP, Voc, Jsc点を表示。
1280 引数:
1281 :param dark_sweeps: 暗電流測定のI-Vスイープデータのリスト。
1282 :type dark_sweeps: list[tuple[numpy.ndarray, numpy.ndarray]]
1283 :param light_sweeps: 光照射測定のI-Vスイープデータのリスト。
1284 :type light_sweeps: list[tuple[numpy.ndarray, numpy.ndarray]]
1285 :param S: デバイスの面積 (cm^2)。
1286 :type S: float
1287 :param light_metrics: 光照射下のI-Vデータから計算された光起電力性能指標。
1288 :type light_metrics: dict
1289 :param dark_params: 暗電流I-Vデータから推定された初期パラメータ。Rs/Rshの接線表示に使用。
1290 :type dark_params: dict or None
1291 :param light_params: 光照射I-Vデータから推定された初期パラメータ。I0/IPVの表示に使用。
1292 :type light_params: dict or None
1293 :param outfile: プロットを保存するファイルパス。Noneの場合、画面に表示します。
1294 :type outfile: str or None
1295 :param T: 温度 (K)。ndiode曲線の計算に使用。
1296 :type T: float
1297 :param pause: プロットを閉じるまで一時停止するかどうか。Trueの場合、plt.show()を呼び出します。
1298 :type pause: bool
1299 戻り値:
1300 :returns: 作成されたMatplotlib Figureオブジェクト。
1301 :rtype: matplotlib.figure.Figure
1302 """
1303 fig = plt.figure(figsize=(13, 8))
1304 gs = fig.add_gridspec(2, 2, height_ratios=[3, 1.5])
1305
1306 ax_log = fig.add_subplot(gs[0, 0])
1307 ax_lin_f = fig.add_subplot(gs[0, 1])
1308 ax_nd = fig.add_subplot(gs[1, 0], sharex=ax_log)
1309 ax_pv = fig.add_subplot(gs[1, 1])
1310
1311 ax_lin_r = ax_lin_f.twinx()
1312
1313 def style_ax(ax, xlabel=None, ylabel=None, title=None):
1314 ax.tick_params(axis='x', labelsize=fontsize)
1315 ax.tick_params(axis='y', labelsize=fontsize)
1316 if xlabel is not None:
1317 ax.set_xlabel(xlabel, fontsize=fontsize)
1318 if ylabel is not None:
1319 ax.set_ylabel(ylabel, fontsize=fontsize)
1320 if title is not None:
1321 ax.set_title(title, fontsize=fontsize)
1322
1323 def plot_split_linear(ax_f, ax_r, sweeps, label_prefix, color=None):
1324 first_f = True
1325 first_r = True
1326 for V, I in sweeps:
1327 V = np.asarray(V, dtype=float)
1328 I = np.asarray(I, dtype=float)
1329
1330 mask_f = V >= 0.0
1331 mask_r = V <= 0.0
1332
1333 if np.any(mask_f):
1334 Jf = I[mask_f] / S * 1e3
1335 ax_f.plot(np.abs(V[mask_f]), Jf, "-", linewidth=1.2, color=color,
1336 label=label_prefix if first_f else None)
1337 first_f = False
1338
1339 if np.any(mask_r):
1340 Jr = I[mask_r] / S * 1e3
1341 ax_r.plot(np.abs(V[mask_r]), Jr, "--", linewidth=1.2, color=color,
1342 label=f"{label_prefix} (rev)" if first_r else None)
1343 first_r = False
1344
1345 plot_split_linear(ax_lin_f, ax_lin_r, dark_sweeps, "dark")
1346 plot_split_linear(ax_lin_f, ax_lin_r, light_sweeps, "light")
1347
1348 if dark_params is not None and len(dark_sweeps) > 0:
1349 Vd, Id = dark_sweeps[0]
1350 rs_info = estimate_rs_tangent_point(Vd, Id)
1351 a_rs = rs_info["slope"] / S * 1e3
1352 b_rs = rs_info["intercept"] / S * 1e3
1353 v_rep_rs = rs_info["v_rep"]
1354 i_rep_rs = rs_info["i_rep"] / S * 1e3
1355
1356 if abs(a_rs) > 1e-30:
1357 v_zero = float(-b_rs / a_rs)
1358 x0 = max(0.0, min(abs(v_rep_rs), abs(v_zero)))
1359 x1 = max(abs(v_rep_rs), abs(v_zero))
1360 x_tan_rs = np.array([x0, x1])
1361 y_tan_rs = a_rs * x_tan_rs + b_rs
1362 ax_lin_f.plot(x_tan_rs, y_tan_rs, "-", linewidth=0.9, label="Rs tangent")
1363 ax_lin_f.plot(abs(v_rep_rs), i_rep_rs, "o", markersize=5, label="Rs point")
1364
1365 if dark_params is not None and len(dark_sweeps) > 0:
1366 Vd, Id = dark_sweeps[0]
1367 tang = estimate_rsh_tangent_point(Vd, Id)
1368 a = tang["slope"] / S * 1e3
1369 b = tang["intercept"] / S * 1e3
1370 v_rep = tang["v_rep"]
1371 i_rep = tang["i_rep"] / S * 1e3
1372
1373 x_rev_all = []
1374 for VV, II in dark_sweeps:
1375 VV = np.asarray(VV, dtype=float)
1376 mask_r_all = VV <= 0.0
1377 if np.any(mask_r_all):
1378 x_rev_all.extend(list(np.abs(VV[mask_r_all])))
1379 if len(x_rev_all) == 0:
1380 x_tan = np.array([0.0, abs(v_rep)])
1381 else:
1382 x_tan = np.array([0.0, float(np.max(x_rev_all))])
1383
1384 y_tan = (-a) * x_tan + b
1385 ax_lin_r.plot(x_tan, y_tan, "-", linewidth=0.9, label="Rsh tangent", zorder=5)
1386 ax_lin_r.plot(abs(v_rep), i_rep, "o", markersize=5, label="Rsh point", zorder=6)
1387 ish0 = float(b)
1388 ax_lin_r.plot(0.0, ish0, "s", markersize=5, label="Ish(0)", zorder=6)
1389
1390 style_ax(ax_lin_f, xlabel="|V| / V", ylabel="J forward / mA cm$^{-2}$, V >= 0",
1391 title="Linear J-V (forward / reverse separated)")
1392 ax_lin_r.tick_params(axis='y', labelsize=fontsize)
1393 ax_lin_r.set_ylabel("J reverse / mA cm$^{-2}$, V <= 0", fontsize=fontsize)
1394 ax_lin_f.grid(True)
1395
1396 h1, l1 = ax_lin_f.get_legend_handles_labels()
1397 h2, l2 = ax_lin_r.get_legend_handles_labels()
1398 ax_lin_f.legend(h1 + h2, l1 + l2, fontsize=8, loc="best")
1399
1400 for i, (V, I) in enumerate(dark_sweeps):
1401 ax_log.plot(V, np.log10(np.abs(I) + EPS_I), "-", linewidth=1.2, label="dark" if i == 0 else None)
1402 for i, (V, I) in enumerate(light_sweeps):
1403 ax_log.plot(V, np.log10(np.abs(I) + EPS_I), "-", linewidth=1.2, label="light" if i == 0 else None)
1404 if light_params is not None:
1405 ax_log.axhline(np.log10(abs(light_params["I0"]) + EPS_I), linestyle=":", linewidth=1.0, label="-I0")
1406 ax_log.axhline(np.log10(abs(light_params["I0"] + light_params["IPV"]) + EPS_I),
1407 linestyle=":", linewidth=1.0, label="-I0-IPV")
1408 style_ax(ax_log, xlabel="V / V", ylabel=r"log$_{10}$(|I|)", title="Dark / Illuminated log10(|I|)-V")
1409 ax_log.grid(True)
1410 ax_log.legend(fontsize=8)
1411
1412 Vdark0, Idark0 = dark_sweeps[0]
1413 Vn, nn = estimate_ndiode_curve(Vdark0, Idark0, T=T)
1414 if len(Vn) > 0:
1415 ax_nd.plot(Vn, nn, "-", linewidth=1.3, label="ndiode(V)")
1416
1417 ish_for_nd = dark_params.get("Ish", None) if dark_params is not None else None
1418 nd_rep = estimate_ndiode_representative_point(Vdark0, Idark0, T=T, Ish=ish_for_nd)
1419 if nd_rep is not None and np.isfinite(nd_rep["ndiode"]):
1420 ax_nd.plot(nd_rep["v_rep"], nd_rep["ndiode"], "o", markersize=6, label="ndiode point")
1421
1422 style_ax(ax_nd, xlabel="V / V", ylabel="ndiode", title="Estimated ndiode vs V")
1423 ax_nd.grid(True)
1424 ax_nd.legend(fontsize=8)
1425
1426 V0, I0 = light_sweeps[0]
1427 J0_mA = I0 / S * 1e3
1428 ypv = -J0_mA
1429 ax_pv.plot(V0, ypv, "-", linewidth=1.5, label="-J (light)")
1430 ax_pv.plot(light_metrics["Vop_V"], -light_metrics["Jop_mA_cm2"], "o", markersize=7, label="MPP")
1431 ax_pv.plot(light_metrics["Voc_V"], 0.0, "s", markersize=7, label="Voc")
1432 ax_pv.plot(0.0, -light_metrics["Jsc_mA_cm2"], "^", markersize=7, label="Jsc")
1433
1434 style_ax(ax_pv, xlabel="V / V", ylabel="-J / mA/cm$^2$", title="Photovoltaic output")
1435 ymax = max(float(np.max(ypv)), 0.0, float(-light_metrics["Jsc_mA_cm2"]))
1436 ax_pv.set_ylim(0.0, ymax if ymax > 0 else 1.0)
1437 ax_pv.grid(True)
1438 ax_pv.legend(fontsize=8)
1439
1440 fig.tight_layout()
1441 if outfile:
1442 print(f"[I/O] Save plot: {outfile}")
1443 fig.savefig(outfile, dpi=160)
1444
1445 if pause:
1446 plt.show()
1447 else:
1448 plt.pause(0.01)
1449 return fig
1450
1451
1452def exec_alpha(args):
1453 """
1454 概要:
1455 alphaまたはmake_alphaモードでプログラムを実行するメイン関数。
1456 詳細説明:
1457 光学スペクトルファイル（RとTr）を読み込み、吸収率と吸収係数を計算し、
1458 結果をプロットしてExcelファイルに保存します。
1459 引数:
1460 :param args: コマンドライン引数。
1461 :type args: argparse.Namespace
1462 戻り値:
1463 なし
1464 """
1465
1466# E_use = compute_photon_energy_eV(args.E, args.lambda_nm)
1467# P0_use = compute_p0(args.F0, args.P0, args.E, args.lambda_nm)
1468# print_args_and_derived(args, E_use, P0_use)
1469 print_args_and_derived(args)
1470
1471 wl_R, R_percent, info_R = read_optical_spectrum(args.R)
1472 wl_T, T_percent, info_T = read_optical_spectrum(args.Tr)
1473
1474 optical = analyze_optical(wl_R, R_percent, wl_T, T_percent, args.d)
1475
1476 print()
1477 print("=== make_alpha analysis ===")
1478 print(f"d_nm : {args.d:.16g} nm")
1479 print(f"R_mean : {np.mean(optical['R']):.16g}")
1480 print(f"Tr_mean : {np.mean(optical['Tr']):.16g}")
1481 print(f"A_mean : {np.mean(optical['A']):.16g}")
1482 print(f"alpha_mean : {np.mean(optical['alpha_cm^-1']):.16g} cm^-1")
1483
1484 print()
1485 outxlsx = f"{args.outprefix}_alpha.xlsx"
1486 if args.mode == 'make_alpha':
1487 save_alpha_to_excel(optical, outxlsx, meta={"d_nm": args.d, "R": args.R, "Tr": args.Tr})
1488 print(f"Saving Excel : {outxlsx}")
1489
1490 outplot = f"{args.outprefix}_alpha_E.png"
1491 print(f"Saving plot : {outplot}")
1492
1493 plot_alpha(optical, outfile=outplot, pause = False)
1494 input("\nPress ENTER to terminate\n")
1495
1496
1497def exec_analyze(args):
1498 """
1499 概要:
1500 analyzeモードでプログラムを実行するメイン関数。
1501 詳細説明:
1502 暗電流および光照射I-Vファイルを読み込み、初期パラメータと光起電力性能指標を推定し、
1503 量子効率を計算します。光学データが利用可能な場合はそれも考慮に入れます。
1504 結果を標準出力に表示し、比較プロットを作成します。
1505 引数:
1506 :param args: コマンドライン引数。
1507 :type args: argparse.Namespace
1508 戻り値:
1509 なし
1510 例外:
1511 :raises ValueError: alpha_at_energy 関数から NaN が返された場合。
1512 """
1513 print()
1514 print("[MODE] analyze")
1515 P0_use = compute_p0(args.F0, args.P0, args.E, args.lambda_nm)
1516 E_use = compute_photon_energy_eV(args.E, args.lambda_nm)
1517 print_args_and_derived(args, E_use, P0_use)
1518
1519 dark_xs, dark_ys, dark_info = read_data(args.dark)
1520 light_xs, light_ys, light_info = read_data(args.light)
1521
1522 dark_sweeps = list(zip(dark_xs, dark_ys))
1523 light_sweeps = list(zip(light_xs, light_ys))
1524
1525 Vd, Id = choose_sweep(dark_xs, dark_ys, args.sweep_dark)
1526 Vl, Il = choose_sweep(light_xs, light_ys, args.sweep_light)
1527
1528 dark_params = estimate_initial_params(Vd, Id, T=args.T)
1529 dark_params["IPV"] = 0.0
1530 light_params = estimate_initial_params(Vl, Il, T=args.T)
1531
1532 dark_metrics = pv_metrics_from_iv(Vd, Id, args.S)
1533 light_metrics = pv_metrics_from_iv(Vl, Il, args.S)
1534
1535 optical = None
1536 alpha_at_E = float("nan")
1537 A_at_E = float("nan")
1538 if args.alpha:
1539 optical = read_alpha_from_excel(args.alpha)
1540
1541 if optical is not None:
1542 alpha_at_E = alpha_at_energy(optical, E_use)
1543 if alpha_at_E is None:
1544 raise ValueError(f"Got nan for alpha_at_E()")
1545
1546 A_at_E = absorptance_from_alpha(alpha_at_E, args.d)
1547
1548 eta = float("nan")
1549 if P0_use is not None and P0_use > 0:
1550 eta = light_metrics["Pmax_W_cm2"] / P0_use
1551
1552 JPV_A_cm2 = light_params["IPV"] / args.S
1553 qe = quantum_efficiencies(args.F0, A_at_E, JPV_A_cm2, light_metrics["Jsc_A_cm2"])
1554
1555 print("=== Device constants ===")
1556 print(f"d = {args.d:.6f} nm")
1557 print(f"S = {args.S:.8f} cm^2")
1558 print(f"T = {args.T:.6f} K")
1559 print()
1560
1561 print("=== Photon input ===")
1562 print(f"lambda = {args.lambda_nm:.16g} nm")
1563 print(f"E = {E_use:.16g} eV")
1564 print(f"F0 = {args.F0:.16g} cm^-2 s^-1" if args.F0 is not None else "F0 = None")
1565 print(f"P0 = {P0_use*1e3:.16g} mW/cm^2" if P0_use is not None else "P0 = None")
1566 print()
1567
1568 print("=== Dark IV estimated parameters ===")
1569 for k in PARAM_NAMES:
1570 print(f"{k:8s} = {dark_params[k]:.16g}")
1571 rsh_info = estimate_rsh_tangent_point(Vd, Id)
1572 rs_info = estimate_rs_tangent_point(Vd, Id)
1573 print(f"{'Vs':8s} = {rs_info['v_rep']:.16g}")
1574 print(f"{'Is':8s} = {rs_info['i_rep']:.16g}")
1575 print(f"{'Vsh':8s} = {rsh_info['v_rep']:.16g}")
1576 print(f"{'Ish':8s} = {rsh_info['Ish']:.16g}")
1577 print(f"{'Ish(0)':8s} = {rsh_info['intercept']:.16g}")
1578 print()
1579
1580 print("=== Illuminated IV estimated parameters ===")
1581 for k in PARAM_NAMES:
1582 print(f"{k:8s} = {light_params[k]:.16g}")
1583 print()
1584
1585 print("=== Photovoltaic metrics ===")
1586 print(f"Voc = {light_metrics['Voc_V']:.16g} V")
1587 print(f"Jsc = {light_metrics['Jsc_mA_cm2']:.16g} mA/cm^2")
1588 print(f"Vop = {light_metrics['Vop_V']:.16g} V")
1589 print(f"Jop = {light_metrics['Jop_mA_cm2']:.16g} mA/cm^2")
1590 print(f"FF = {light_metrics['FF']:.16g}")
1591 print(f"Pmax = {light_metrics['Pmax_mW_cm2']:.16g} mW/cm^2")
1592 print(f"Efficiency = {eta:.16g}" if np.isfinite(eta) else "Efficiency = nan")
1593 print()
1594
1595 if optical is not None:
1596 print("=== Optical metrics ===")
1597 print(f"alpha(E) = {alpha_at_E:.16g} cm^-1")
1598 print(f"absorption length = {args.d} nm")
1599 print(f"A(E) = {A_at_E:.16g}")
1600 print()
1601
1602 print("=== Quantum efficiencies ===")
1603 print(f"JPV = {JPV_A_cm2*1e3:.16g} mA/cm^2")
1604 print(f"EQE_gen = {qe['EQE_gen']:.16g}")
1605 print(f"IQE_gen = {qe['IQE_gen']:.16g}")
1606 print(f"EQE = {qe['EQE']:.16g}")
1607 print(f"IQE = {qe['IQE']:.16g}")
1608
1609 print()
1610 outplot = f"{args.outprefix}_iv.png"
1611 plot_iv_comparison(
1612 dark_sweeps, light_sweeps, args.S, light_metrics,
1613 dark_params=dark_params, light_params=light_params,
1614 outfile=outplot, T=args.T,
1615 pause=False
1616 )
1617 print(f"Saved plot: {outplot}")
1618 input("\nPress ENTER to terminate\n")
1619
1620
1621_original_print = None
1622_redirect_fp = None
1623def dual_print(*args, **kwargs):
1624 """
1625 概要:
1626 標準出力とログファイルの両方にメッセージを出力します。
1627 詳細説明:
1628 通常のprint関数の動作を模倣し、標準出力に加えて、
1629 _redirect_fpに設定されたファイルオブジェクトにも出力します。
1630 引数:
1631 :param args: printに渡される位置引数。
1632 :type args: tuple
1633 :param kwargs: printに渡されるキーワード引数。
1634 :type kwargs: dict
1635 戻り値:
1636 なし
1637 """
1638 _original_print(*args, **kwargs)
1639
1640 file_kwargs = dict(kwargs)
1641 file_kwargs.pop("file", None) # file指定があっても無視してログファイルへ
1642 file_kwargs.setdefault("flush", True)
1643
1644 _original_print(*args, file=_redirect_fp, **file_kwargs)
1645
1646def main():
1647 """
1648 概要:
1649 プログラムのエントリーポイント。
1650 詳細説明:
1651 コマンドライン引数を初期化し、指定されたモードに基づいて
1652 対応する実行関数を呼び出します。
1653 実行ログは自動的にファイルにリダイレクトされます。
1654 引数:
1655 なし
1656 戻り値:
1657 なし
1658 """
1659
1660 global _original_print, _redirect_fp
1661
1662 print()
1663 args, parser = initialize()
1664
1665 _original_print = builtins.print
1666 if "alpha" in args.mode:
1667 logfile = f"{args.outprefix}_alpha.txt"
1668 else:
1669 logfile = f"{args.outprefix}_{args.mode}.txt"
1670 print(f"Open log file [{logfile}]")
1671 _redirect_fp = open(logfile, "w", encoding="utf-8")
1672 builtins.print = dual_print
1673
1674 if args.mode == "alpha" or args.mode == "make_alpha":
1675 exec_alpha(args)
1676 elif args.mode == "analyze":
1677 exec_analyze(args)
1678 else:
1679 parser.error(f"Unsupported mode: {args.mode}")
1680
1681
1682if __name__ == "__main__":
1683 try:
1684 main()
1685 except Exception:
1686 traceback.print_exc()
1687 sys.exit(1)
1688 finally:
1689 if _redirect_fp: _redirect_fp.close()