import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
from scipy.constants import h, c, k


# -----------------------------
# 1. 数据预处理
# -----------------------------
# PDMS 可见光折射率数据
wl_data = np.array([0.3500,0.3535,0.3570,0.3605,0.3640,0.3675,0.3710,0.3745,0.3780,0.3815,
                    0.3850,0.3885,0.3920,0.3955,0.3990,0.4025,0.4060,0.4095,0.4130,0.4165,
                    0.4200,0.4235,0.4270,0.4305,0.4340,0.4375,0.4410,0.4445,0.4480,0.4515,
                    0.4550,0.4585,0.4620,0.4655,0.4690,0.4725,0.4760,0.4795,0.4830,0.4865,
                    0.4900,0.4935,0.4970,0.5005,0.5040,0.5075,0.5110,0.5145,0.5180,0.5215,
                    0.5250,0.5285,0.5320,0.5355,0.5390,0.5425,0.5460,0.5495,0.5530,0.5565,
                    0.5600,0.5635,0.5670,0.5705,0.5740,0.5775,0.5810,0.5845,0.5880,0.5915,
                    0.5950,0.5985,0.6020,0.6055,0.6090,0.6125,0.6160,0.6195,0.6230,0.6265,
                    0.6300,0.6335,0.6370,0.6405,0.6440,0.6475,0.6510,0.6545,0.6580,0.6615,
                    0.6650,0.6685,0.6720,0.6755,0.6790,0.6825,0.6860,0.6895,0.6930,0.6965,
                    0.7000])
n_data = np.array([1.4585,1.4576,1.4567,1.4559,1.4550,1.4542,1.4535,1.4527,1.45197,1.45126,
                   1.45057,1.44990,1.44926,1.44863,1.44802,1.44742,1.44685,1.44628,1.44574,1.44521,
                   1.44470,1.44420,1.44371,1.44324,1.44277,1.44232,1.44188,1.44145,1.44104,1.44063,
                   1.44024,1.43985,1.43947,1.43911,1.43875,1.43840,1.43805,1.43772,1.43739,1.43707,
                   1.43676,1.43645,1.43616,1.43587,1.43558,1.43531,1.43503,1.43477,1.43450,1.43425,
                   1.43399,1.43375,1.43351,1.43328,1.43305,1.43283,1.43260,1.43238,1.43217,1.43197,
                   1.43177,1.43157,1.43137,1.43118,1.43098,1.43080,1.43062,1.43044,1.43027,1.43009,
                   1.42993,1.42976,1.42960,1.42944,1.42929,1.42913,1.42898,1.42883,1.42869,1.42855,
                   1.42841,1.42827,1.42813,1.42799,1.42787,1.42773,1.42761,1.42749,1.42736,1.42724,
                   1.42712,1.42701,1.42689,1.42677,1.42666,1.42656,1.42645,1.42634,1.42624,1.42613,
                   1.42604])

# 三次样条插值
cs_n = CubicSpline(wl_data, n_data)

# 定义厚度序列（μm）
thicknesses = [0.5, 1.0, 1.5, 2.0]

# -----------------------------
# 2. 发射率计算（小问1）
# -----------------------------
def fresnel_reflectance(n1, n2):
    return ((n1 - n2)/(n1 + n2))**2

def thin_film_reflectance(n_film, d, wl):
    R12 = fresnel_reflectance(1.0, n_film)
    R23 = fresnel_reflectance(n_film, 1.0)
    delta = 2 * np.pi * n_film * d / wl
    R = (R12 + R23 + 2*np.sqrt(R12*R23)*np.cos(2*delta)) / (1 + R12*R23 + 2*np.sqrt(R12*R23)*np.cos(2*delta))
    return R

# 波长范围 0.35-0.7 μm，步长 0.001
wl_fine = np.linspace(0.35, 0.7, 500)
plt.figure(figsize=(8,5))

emission_dict = {}
for d in thicknesses:
    R = thin_film_reflectance(cs_n(wl_fine), d, wl_fine)
    epsilon = 1 - R
    emission_dict[d] = epsilon
    plt.plot(wl_fine, epsilon, label=f"d={d} μm")

plt.xlabel("Wavelength (μm)")
plt.ylabel("Emissivity ε(λ)")
plt.title("PDMS Thin Film Spectral Emissivity")
plt.legend()
plt.grid(True)
plt.show()

# -----------------------------
# 问题2
# -----------------------------
T_film = 420  # 薄膜温度 K
T_sky = 280   # 天空温度 K
alpha = 0.1   # 吸收可见光比例
I_sun_total = 1000  # W/m²，太阳总辐射
I_sun = np.ones_like(wl_fine) * I_sun_total / len(wl_fine)  # 均匀分布
tau_atm = 0.9  # 大气透射率假设

# -----------------------------
# 4. Planck 黑体辐射函数
# -----------------------------
def planck_spectrum(wl, T):
    wl_m = wl * 1e-6  # μm → m
    return (2*h*c**2 / wl_m**5) / (np.exp(h*c/(wl_m*k*T)) - 1)

# -----------------------------
# 5. 计算净辐射功率
# -----------------------------
P_net_list = []
for d in thicknesses:
    epsilon = emission_dict[d]
    P_emit = np.trapz(epsilon * planck_spectrum(wl_fine, T_film), wl_fine)
    P_sun = np.trapz(alpha * I_sun, wl_fine)
    P_atm = np.trapz(epsilon * planck_spectrum(wl_fine, T_sky) * tau_atm, wl_fine)
    P_net = P_emit - P_sun - P_atm
    P_net_list.append(P_net)

# -----------------------------
# 6. 绘图：净辐射功率 vs 厚度
# -----------------------------
plt.figure(figsize=(8,5))
plt.plot(thicknesses, P_net_list, marker='o', linestyle='-', color='blue', label='Net Radiative Cooling')
for x, y in zip(thicknesses, P_net_list):
    plt.text(x, y, f"{y:.1f}", ha='center', va='bottom')

plt.xlabel("Film Thickness (μm)")
plt.ylabel("Net Radiative Cooling Power (W/m²)")
plt.title("PDMS Thin Film Radiative Cooling Performance")
plt.grid(True)
plt.legend()
plt.show()