198 lines
7.0 KiB
Python
198 lines
7.0 KiB
Python
import numpy as np
|
||
import matplotlib.pyplot as plt
|
||
from scipy.interpolate import CubicSpline
|
||
|
||
def make_strictly_increasing(wl, n, k):
|
||
# 去除完全相同的点
|
||
unique_wl, indices = np.unique(wl, return_index=True)
|
||
if len(unique_wl) != len(wl):
|
||
print(f"Removed {len(wl) - len(unique_wl)} duplicate wavelength points")
|
||
wl, n, k = wl[indices], n[indices], k[indices]
|
||
|
||
# 确保严格递增
|
||
is_increasing = np.diff(wl) > 0
|
||
if not all(is_increasing):
|
||
# 移除不满足严格递增的点
|
||
valid_indices = np.concatenate([[True], is_increasing])
|
||
wl, n, k = wl[valid_indices], n[valid_indices], k[valid_indices]
|
||
|
||
return wl, n, k
|
||
# -----------------------------
|
||
# 1. 从data.txt读取分块格式数据（先wl+n，再wl+k）
|
||
# -----------------------------
|
||
def read_split_data(file_path):
|
||
"""
|
||
解析分块数据格式：
|
||
第一部分：wl n（多行数据）
|
||
第二部分：wl k（多行数据）
|
||
返回：sorted_wl, n, k（波长已排序，保证n和k一一对应）
|
||
"""
|
||
with open(file_path, 'r', encoding='utf-8') as f:
|
||
lines = [line.strip() for line in f if line.strip() and not line.startswith('#')]
|
||
|
||
# 分割n数据和k数据（以"wl k"为分界）
|
||
split_idx = None
|
||
for i, line in enumerate(lines):
|
||
if line == 'wl k': # 找到k数据的表头
|
||
split_idx = i
|
||
break
|
||
|
||
# 提取n数据（表头后到split_idx前）
|
||
n_lines = lines[1:split_idx] # 跳过"wl n"表头
|
||
wl_n = []
|
||
n_list = []
|
||
for line in n_lines:
|
||
wl, n_val = line.split()
|
||
wl_n.append(float(wl))
|
||
n_list.append(float(n_val))
|
||
|
||
# 提取k数据（split_idx后）
|
||
k_lines = lines[split_idx + 1:] # 跳过"wl k"表头
|
||
wl_k = []
|
||
k_list = []
|
||
for line in k_lines:
|
||
wl, k_val = line.split()
|
||
wl_k.append(float(wl))
|
||
k_list.append(float(k_val))
|
||
|
||
# 转换为numpy数组
|
||
wl_n = np.array(wl_n)
|
||
n_list = np.array(n_list)
|
||
wl_k = np.array(wl_k)
|
||
k_list = np.array(k_list)
|
||
|
||
# 确保n和k的波长完全一致（否则插值会出错）
|
||
assert np.allclose(wl_n, wl_k), "n和k的波长列表不一致！"
|
||
|
||
# 排序（按波长递增，避免插值异常）
|
||
sorted_idx = np.argsort(wl_n)
|
||
sorted_wl = wl_n[sorted_idx]
|
||
sorted_n = n_list[sorted_idx]
|
||
sorted_k = k_list[sorted_idx]
|
||
|
||
return sorted_wl, sorted_n, sorted_k
|
||
|
||
|
||
# 读取数据
|
||
wl_all, n_all, k_all = read_split_data('/Users/spasolreisa/IdeaProjects/asiaMath/data.txt')
|
||
|
||
wl_all, n_all, k_all = make_strictly_increasing(wl_all, n_all, k_all)
|
||
#
|
||
# 三次样条插值（覆盖全波段，保证计算精度）
|
||
cs_n = CubicSpline(wl_all, n_all) # 折射率n的插值函数
|
||
cs_k = CubicSpline(wl_all, k_all) # 消光系数k的插值函数
|
||
|
||
# 定义待研究的PDMS薄膜厚度（μm），可按需调整
|
||
thicknesses = [0.5, 1.0, 1.5, 2.0]
|
||
|
||
|
||
# -----------------------------
|
||
# 2. 核心物理模型：薄膜发射率计算
|
||
# -----------------------------
|
||
def fresnel_reflectance(n1, k1, n2, k2):
|
||
"""
|
||
菲涅尔反射率（垂直入射近似，考虑消光系数k）
|
||
输入：介质1的n/k，介质2的n/k
|
||
输出：垂直入射时的反射率R（0-1）
|
||
"""
|
||
m1 = n1 + 1j * k1 # 介质1的复折射率
|
||
m2 = n2 + 1j * k2 # 介质2的复折射率
|
||
return np.abs((m1 - m2) / (m1 + m2)) ** 2
|
||
|
||
|
||
def thin_film_emissivity(n_film, k_film, d, wl, n_air=1.0, k_air=0.0):
|
||
"""
|
||
薄膜发射率计算（考虑多光束干涉和吸收）
|
||
输入：
|
||
n_film/k_film: 薄膜的折射率/消光系数
|
||
d: 薄膜厚度（μm）
|
||
wl: 波长（μm）
|
||
n_air/k_air: 空气的折射率/消光系数（默认n=1, k=0）
|
||
输出：
|
||
epsilon: 发射率（0-1）
|
||
"""
|
||
# 1. 计算上下表面的菲涅尔反射率
|
||
R12 = fresnel_reflectance(n_air, k_air, n_film, k_film) # 空气→薄膜
|
||
R23 = fresnel_reflectance(n_film, k_film, n_air, k_air) # 薄膜→空气
|
||
|
||
# 2. 计算相位差和吸收衰减
|
||
delta = 2 * np.pi * n_film * d / wl # 干涉相位差（无吸收时）
|
||
alpha = 4 * np.pi * k_film * d / wl # 吸收导致的振幅衰减系数
|
||
|
||
# 3. 多光束干涉反射率（考虑吸收）
|
||
numerator = R12 + R23 * np.exp(-alpha) + 2 * np.sqrt(R12 * R23 * np.exp(-alpha)) * np.cos(2 * delta)
|
||
denominator = 1 + R12 * R23 * np.exp(-alpha) + 2 * np.sqrt(R12 * R23 * np.exp(-alpha)) * np.cos(2 * delta)
|
||
R_total = numerator / denominator
|
||
|
||
# 4. 多光束干涉透射率（考虑吸收）
|
||
T_total = (1 - R12) * (1 - R23) * np.exp(-alpha) / denominator
|
||
|
||
# 5. 基尔霍夫定律：ε = 1 - R - T（局部热平衡）
|
||
epsilon = 1 - R_total - T_total
|
||
return epsilon
|
||
|
||
|
||
# -----------------------------
|
||
# 3. 计算并绘制不同厚度的发射率光谱
|
||
# -----------------------------
|
||
# 定义计算波长范围（覆盖数据的有效波长区间，步长0.01μm保证平滑）
|
||
wl_min = wl_all.min()
|
||
wl_max = wl_all.max()
|
||
wl_fine = np.linspace(wl_min, wl_max, int((wl_max - wl_min) / 0.01) + 1)
|
||
|
||
# 创建绘图
|
||
plt.figure(figsize=(12, 7))
|
||
plt.rcParams['font.sans-serif'] = ['Arial'] # 统一字体
|
||
|
||
# 存储不同厚度的发射率结果（用于后续分析）
|
||
emission_dict = {}
|
||
|
||
for d in thicknesses:
|
||
# 插值得到当前波长下的n和k
|
||
n_film = cs_n(wl_fine)
|
||
k_film = cs_k(wl_fine)
|
||
|
||
# 计算发射率
|
||
epsilon = thin_film_emissivity(n_film, k_film, d, wl_fine)
|
||
emission_dict[d] = epsilon
|
||
|
||
# 绘制光谱曲线
|
||
plt.plot(wl_fine, epsilon, linewidth=2, label=f'Thickness = {d} μm')
|
||
|
||
# -----------------------------
|
||
# 4. 图表美化与标注（突出辐射制冷关键波段）
|
||
# -----------------------------
|
||
# 标注大气透明窗口（8-13μm，辐射制冷核心波段）
|
||
if wl_min <= 13 and wl_max >= 8: # 仅当数据覆盖该波段时才标注
|
||
plt.axvspan(8, 13, alpha=0.15, color='red', label='Atmospheric Window (8-13 μm)')
|
||
|
||
# 坐标轴与标题
|
||
plt.xlabel('Wavelength (μm)', fontsize=14)
|
||
plt.ylabel('Emissivity ε(λ)', fontsize=14)
|
||
plt.title('PDMS Thin Film Spectral Emissivity for Different Thicknesses', fontsize=16, fontweight='bold')
|
||
|
||
# 网格、图例与范围设置
|
||
plt.grid(True, alpha=0.3, linestyle='--')
|
||
plt.legend(fontsize=12, loc='best')
|
||
plt.ylim(0, 1.05) # 发射率范围0-1.05（留出余量）
|
||
plt.xlim(wl_min, wl_max)
|
||
|
||
# 保存图片（高分辨率）
|
||
plt.tight_layout()
|
||
plt.savefig('PDMS_emissivity_spectrum.png', dpi=300, bbox_inches='tight')
|
||
plt.show()
|
||
|
||
# -----------------------------
|
||
# 5. 输出关键波段（8-13μm）的平均发射率（若数据覆盖该波段）
|
||
# -----------------------------
|
||
if wl_min <= 13 and wl_max >= 8:
|
||
print("=== 8-13 μm 大气窗口平均发射率 ===")
|
||
wl_window = np.linspace(8, 13, 500) # 大气窗口内的波长点
|
||
for d in thicknesses:
|
||
n_window = cs_n(wl_window)
|
||
k_window = cs_k(wl_window)
|
||
epsilon_window = thin_film_emissivity(n_window, k_window, d, wl_window)
|
||
avg_epsilon = np.mean(epsilon_window)
|
||
print(f"厚度 {d} μm: 平均发射率 = {avg_epsilon:.4f}")
|
||
else:
|
||
print("数据未覆盖8-13μm大气窗口，跳过该波段平均发射率计算。") |