"""
太陽光スペクトル実測データとPlanckの黒体放射公式の比較モジュール。
概要:
Excelファイルから太陽光スペクトル（AM1.5等）を読み込み、
指定された温度におけるPlanckの法則に基づく理論曲線と比較プロットします。
詳細説明:
1. 'solar.xlsx'から波長とエネルギー密度のデータを読み込みます。
（AM1.5 全天光データをデフォルトとしています）
2. Planckの放射公式を用いて、特定の温度 $T$ における分光放射輝度を計算します。
3. 実測データ（太陽光）と理論データ（黒体）を同一グラフ上に表示し、
太陽光がいかに黒体放射に近いか（あるいは大気吸収でどう変化するか）を視覚化します。
理論式:
Planckの公式による分光放射密度 $u(\lambda, T)$:
$$u(\lambda, T) = \frac{8\pi hc}{\lambda^5} \frac{1}{e^{hc/\lambda k_B T} - 1}$$
関連リンク: :doc:`blackbody_radiation_usage`
"""
import sys
import openpyxl as px
from math import exp
import numpy as np
import matplotlib.pyplot as plt
#=============================
# 物理定数
#=============================
PI = 3.14159265358979323846
H = 6.6260755e-34 # Js
C = 2.99792458e8 # m/s
KB = 1.380658e-23 # J/K
#=============================
# 設定
#=============================
SPECTRUM_FILE_DEF = 'solar.xlsx'
TEMPERATURE_DEF = 6000.0 # K (太陽表面温度の目安)
[ドキュメント]
def calculate_planck(wavelength_nm, temperature):
"""
指定された波長と温度におけるPlanckの公式を計算します。
:param wavelength_nm: float: 波長 [nm]
:param temperature: float: 黒体温度 [K]
:returns: float: エネルギー密度 [W/m^2/nm]
"""
wl_m = wavelength_nm * 1.0e-9
# 指数部分の計算（オーバーフロー対策）
try:
exponent = (H * C) / (KB * temperature * wl_m)
if exponent > 700: return 0.0 # 極端に短い波長ではゼロ
prefactor = (8.0 * PI * H * C) / (wl_m**5)
# 1.5は立体角や幾何学的因子を考慮したスケーリング定数（近似用）
density = prefactor / (exp(exponent) - 1.0) * 1.0e-9 * 1.5
return density
except ZeroDivisionError:
return 0.0
[ドキュメント]
def main():
"""
主要な処理フロー：データの読み込み、理論計算、可視化。
"""
argv = sys.argv
t_input = float(argv[1]) if len(argv) >= 2 else TEMPERATURE_DEF
print(f"\n--- Blackbody Radiation Analysis (T = {t_input} K) ---")
# Excelファイルの読み込み
try:
book = px.load_workbook(SPECTRUM_FILE_DEF, data_only=True)
sheet = book.active
print(f"Reading spectrum data from [{SPECTRUM_FILE_DEF}] (Sheet: {sheet.title})")
except Exception as e:
print(f"Error: Could not open {SPECTRUM_FILE_DEF}: {e}")
return
x_obs = [] # 実測波長 [nm]
y_obs = [] # 実測強度 [W/m2/nm]
# AM1.5全天光（H列:波長[um], J列:強度[W/m2/um]）を想定
# 2行目から読み込み
for row in range(2, 2000):
w_cell = sheet[f"H{row}"].value
i_cell = sheet[f"J{row}"].value
if w_cell is None or i_cell is None:
break
try:
# 単位変換: um -> nm, W/m2/um -> W/m2/nm
x_obs.append(float(w_cell) * 1000.0)
y_obs.append(float(i_cell) * 0.001)
except ValueError:
continue
# 理論値の計算
y_theory = [calculate_planck(wl, t_input) for wl in x_obs]
#=============================
# グラフ表示
#=============================
plt.figure(figsize=(10, 6))
plt.plot(x_obs, y_obs, label=f"Observed Solar Spectrum ({SPECTRUM_FILE_DEF})", alpha=0.8)
plt.plot(x_obs, y_theory, label=f"Planck Equation (T={t_input}K)", color='red', linestyle='--')
plt.xlabel("Wavelength [nm]")
plt.ylabel("Spectral Irradiance [W/m$^2$/nm]")
plt.title("Solar Spectrum vs Blackbody Radiation")
plt.xlim([0, 2500.0])
plt.ylim([0, max(y_obs) * 1.2 if y_obs else 2.5])
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show(block=False)
print("\nPress ENTER to exit>>", end='')
input()
if __name__ == '__main__':
main()