import importlib.util
import math
import os
from dataclasses import dataclass
from typing import Dict, List, Tuple

import numpy as np
from matplotlib import pyplot as plt
from scipy.interpolate import CubicSpline

plt.rcParams["font.sans-serif"] = ["DejaVu Sans", "Arial", "Helvetica"]
plt.rcParams["axes.unicode_minus"] = False

SIGMA = 5.670374419e-8
T_AMB = 300.0  # K (≈27 ℃)
T_SKY = 280.0  # Clear dry sky equivalent radiation temperature
SOLAR_IRR = 900.0  # W/m^2, clear sky noon
H_CONV = 8.0  # W/m^2/K, natural convection + light wind


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


def read_split_data(file_path):
    """读取分块格式的光学数据"""
    with open(file_path, 'r', encoding='utf-8') as f:
        lines = [line.strip() for line in f if line.strip() and not line.startswith('#')]

    split_idx = None
    for i, line in enumerate(lines):
        if line == 'wl	k':
            split_idx = i
            break

    n_lines = lines[1:split_idx]
    wl_n, n_list = [], []
    for line in n_lines:
        parts = line.split()
        if len(parts) >= 2:
            wl_n.append(float(parts[0]))
            n_list.append(float(parts[1]))

    k_lines = lines[split_idx + 1:]
    wl_k, k_list = [], []
    for line in k_lines:
        parts = line.split()
        if len(parts) >= 2:
            wl_k.append(float(parts[0]))
            k_list.append(float(parts[1]))

    wl_n = np.array(wl_n)
    n_list = np.array(n_list)
    wl_k = np.array(wl_k)
    k_list = np.array(k_list)

    # 检查波长范围是否一致
    if not np.array_equal(wl_n, wl_k):
        print(f"警告: n和k的波长范围不完全一致")
        print(f"n波长范围: {wl_n.min():.2f} - {wl_n.max():.2f} μm")
        print(f"k波长范围: {wl_k.min():.2f} - {wl_k.max():.2f} μm")

        # 找到共同的波长范围
        common_wl = np.intersect1d(wl_n, wl_k)
        if len(common_wl) == 0:
            raise ValueError("n和k没有共同的波长点！")

        # 插值到共同的波长点
        n_common = np.interp(common_wl, wl_n, n_list)
        k_common = np.interp(common_wl, wl_k, k_list)

        wl_n, n_list, wl_k, k_list = common_wl, n_common, common_wl, k_common

    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


def fresnel_reflectance(n1, k1, n2, k2):
    """菲涅尔反射率计算"""
    m1 = n1 + 1j * k1
    m2 = n2 + 1j * k2
    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):
    """精确的薄膜发射率计算"""
    R12 = fresnel_reflectance(n_air, k_air, n_film, k_film)
    R23 = fresnel_reflectance(n_film, k_film, n_air, k_air)

    delta = 2 * np.pi * n_film * d / wl
    alpha = 4 * np.pi * k_film * d / wl

    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

    T_total = (1 - R12) * (1 - R23) * np.exp(-alpha) / denominator

    epsilon = 1 - R_total - T_total
    return np.clip(epsilon, 0, 1)


def solar_absorptance(thickness_um: float, cs_n, cs_k, wl_range) -> float:
    """
    基于精确光学模型计算太阳吸收率
    """
    # 太阳光谱范围 - 确保在数据范围内
    wl_solar_min = max(0.3, wl_range.min())
    wl_solar_max = min(2.5, wl_range.max())

    if wl_solar_min >= wl_solar_max:
        print(f"警告: 数据波长范围 {wl_range.min():.2f}-{wl_range.max():.2f}μm 不覆盖太阳光谱")
        return min(0.05 + 0.001 * thickness_um, 0.15)

    wl_solar = np.linspace(wl_solar_min, wl_solar_max, 200)

    try:
        n_solar = cs_n(wl_solar)
        k_solar = cs_k(wl_solar)
        alpha_spectral = thin_film_emissivity(n_solar, k_solar, thickness_um, wl_solar)
        avg_alpha = np.mean(alpha_spectral)
        return float(np.clip(avg_alpha, 0, 1))

    except Exception as e:
        print(f"太阳吸收率计算失败 {thickness_um}μm: {e}")
        return min(0.05 + 0.001 * thickness_um, 0.15)


@dataclass
class CoolingResult:
    thickness_um: float
    eps_window: float
    alpha_solar: float
    net_power_at_amb: float
    eq_temp_K: float

    @property
    def eq_temp_C(self) -> float:
        return self.eq_temp_K - 273.15

    @property
    def delta_T(self) -> float:
        return self.eq_temp_K - T_AMB


def net_cooling_power(temp_K: float, emissivity: float, alpha_s: float) -> float:
    """计算净冷却功率"""
    radiative = emissivity * SIGMA * (temp_K ** 4 - T_SKY ** 4)
    solar_gain = alpha_s * SOLAR_IRR
    convective = H_CONV * (temp_K - T_AMB)
    return -solar_gain - convective + radiative


def solve_equilibrium(emissivity: float, alpha_s: float) -> float:
    """求解平衡温度"""
    low, high = 250.0, 330.0
    f_low = net_cooling_power(low, emissivity, alpha_s)
    f_high = net_cooling_power(high, emissivity, alpha_s)

    if f_low * f_high > 0:
        return low if abs(f_low) < abs(f_high) else high

    for _ in range(80):
        mid = 0.5 * (low + high)
        f_mid = net_cooling_power(mid, emissivity, alpha_s)
        if abs(f_mid) < 1e-4:
            return mid
        if f_low * f_mid <= 0:
            high, f_high = mid, f_mid
        else:
            low, f_low = mid, f_mid
    return 0.5 * (low + high)


def load_optical_data():
    """加载光学数据并创建插值函数"""
    data_path = '/Users/spasolreisa/IdeaProjects/asiaMath/data.txt'
    print(f"正在加载光学数据: {data_path}")
    wl_all, n_all, k_all = read_split_data(data_path)
    wl_all, n_all, k_all = make_strictly_increasing(wl_all, n_all, k_all)

    print(f"数据波长范围: {wl_all.min():.2f} - {wl_all.max():.2f} μm")
    print(f"数据点数: {len(wl_all)}")
    print(f"折射率n范围: {n_all.min():.3f} - {n_all.max():.3f}")
    print(f"消光系数k范围: {k_all.min():.3f} - {k_all.max():.3f}")

    cs_n = CubicSpline(wl_all, n_all)
    cs_k = CubicSpline(wl_all, k_all)

    return wl_all, cs_n, cs_k


def compute_emissivity_for_thickness(thicknesses_um, cs_n, cs_k, wl_range):
    """计算各厚度在大气窗口的平均发射率"""
    # 大气窗口波长范围 - 确保在数据范围内
    wl_window_min = max(8.0, wl_range.min())
    wl_window_max = min(13.0, wl_range.max())

    if wl_window_min >= wl_window_max:
        print(f"警告: 数据波长范围 {wl_range.min():.2f}-{wl_range.max():.2f}μm 不覆盖大气窗口")
        # 返回默认值
        return {}, {d: 0.8 for d in thicknesses_um}

    wl_window = np.linspace(wl_window_min, wl_window_max, 300)

    emissivity_map = {}
    window_avg = {}

    for d in thicknesses_um:
        try:
            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)

            emissivity_map[d] = epsilon_window
            window_avg[d] = float(avg_epsilon)
            print(f"厚度 {d}μm: 大气窗口发射率 = {avg_epsilon:.4f}")
        except Exception as e:
            print(f"发射率计算失败 {d}μm: {e}")
            emissivity_map[d] = np.full_like(wl_window, 0.8)
            window_avg[d] = 0.8

    return emissivity_map, window_avg


def evaluate() -> Tuple[List[CoolingResult], str]:
    """主评估函数"""
    # 加载光学数据
    try:
        wl_all, cs_n, cs_k = load_optical_data()
    except Exception as e:
        print(f"加载光学数据失败: {e}")
        return [], ""

    # 定义研究的厚度 - 使用与原始代码相同的厚度
    THICKNESSES = [1, 5, 10, 25, 50, 100]

    # 计算发射率
    print("\n计算大气窗口发射率...")
    emissivity_map, window_avg = compute_emissivity_for_thickness(THICKNESSES, cs_n, cs_k, wl_all)

    results: List[CoolingResult] = []
    for d in THICKNESSES:
        eps_window = window_avg[d]
        alpha_s = solar_absorptance(d, cs_n, cs_k, wl_all)
        net_amb = net_cooling_power(T_AMB, eps_window, alpha_s)
        eq_temp = solve_equilibrium(eps_window, alpha_s)

        results.append(
            CoolingResult(
                thickness_um=d,
                eps_window=eps_window,
                alpha_solar=alpha_s,
                net_power_at_amb=net_amb,
                eq_temp_K=eq_temp,
            )
        )

    # 保存结果
    outdir = os.path.join(os.path.dirname(__file__), "outputs")
    os.makedirs(outdir, exist_ok=True)

    csv_path = os.path.join(outdir, "question2_cooling_summary_advanced.csv")
    with open(csv_path, "w", encoding="utf-8") as f:
        f.write("thickness_um,eps_8_13,alpha_solar,net_power_amb_Wm2,eq_temp_C,delta_T_C\n")
        for res in results:
            f.write(
                f"{res.thickness_um},{res.eps_window:.4f},{res.alpha_solar:.3f},"
                f"{res.net_power_at_amb:.2f},{res.eq_temp_C:.2f},{res.delta_T:.2f}\n"
            )

    fig_path = os.path.join(outdir, "question2_cooling_results_advanced.png")
    plot_results(results, fig_path)

    # 绘制发射率光谱
    plot_emissivity_spectrum(THICKNESSES, cs_n, cs_k, wl_all, outdir)

    return results, csv_path


def plot_results(results: List[CoolingResult], fig_path: str) -> None:
    """绘制冷却性能结果 - 使用与原始代码相同的格式"""
    thickness = [r.thickness_um for r in results]
    net_power = [r.net_power_at_amb for r in results]
    delta_T = [r.delta_T for r in results]

    fig, ax1 = plt.subplots(figsize=(9, 5))

    # 使用与原始代码相同的条形图格式
    ax1.bar(thickness, net_power, width=4, alpha=0.6, label="Net Cooling Power @T_amb")
    ax1.set_xlabel("PDMS Film Thickness (µm)")
    ax1.set_ylabel("Net Cooling Power (W/m²)")
    ax1.axhline(0, color="black", linewidth=0.8)

    ax2 = ax1.twinx()
    ax2.plot(thickness, delta_T, color="tab:red", marker="o",
             linewidth=2, markersize=6, label="Equilibrium Temperature Difference")
    ax2.set_ylabel("Equilibrium Temperature Difference (K)")

    lines, labels = ax1.get_legend_handles_labels()
    lines2, labels2 = ax2.get_legend_handles_labels()
    ax1.legend(lines + lines2, labels + labels2, loc="upper right")
    ax1.grid(True, alpha=0.3)
    fig.tight_layout()
    plt.savefig(fig_path, dpi=300)
    plt.close(fig)
    print(f"冷却性能图保存至: {fig_path}")


def plot_emissivity_spectrum(thicknesses, cs_n, cs_k, wl_range, outdir):
    """绘制发射率光谱"""
    # 使用原始数据的完整波长范围
    wl_min = wl_range.min()
    wl_max = wl_range.max()
    wl_fine = np.linspace(wl_min, wl_max, 1000)

    plt.figure(figsize=(12, 7))

    for d in thicknesses:
        try:
            n_film = cs_n(wl_fine)
            k_film = cs_k(wl_fine)
            epsilon = thin_film_emissivity(n_film, k_film, d, wl_fine)
            plt.plot(wl_fine, epsilon, linewidth=2, label=f'{d} μm')
        except Exception as e:
            print(f"绘图失败 {d}μm: {e}")
            continue

    # 标注大气窗口（如果数据覆盖）
    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', fontsize=16, fontweight='bold')
    plt.grid(True, alpha=0.3, linestyle='--')
    plt.legend(fontsize=10, loc='best')
    plt.ylim(0, 1.05)
    plt.xlim(wl_min, wl_max)

    fig_path = os.path.join(outdir, "advanced_emissivity_spectrum.png")
    plt.tight_layout()
    plt.savefig(fig_path, dpi=300, bbox_inches='tight')
    plt.close()
    print(f"发射率光谱图保存至: {fig_path}")


def main():
    """主函数"""
    print("开始辐射冷却性能分析...")
    results, csv_path = evaluate()

    if results:
        print(f"\n结果已保存至: {csv_path}")
        print("\n=== 辐射冷却性能汇总 ===")
        print("厚度(μm) | 大气窗口发射率 | 太阳吸收率 | 净冷却功率(W/m²) | 平衡温度(°C) | ΔT(K)")
        print("-" * 85)
        for r in results:
            print(
                f"{r.thickness_um:>7.1f} | {r.eps_window:>13.3f} | {r.alpha_solar:>10.3f} | "
                f"{r.net_power_at_amb:>15.1f} | {r.eq_temp_C:>11.1f} | "
                f"{r.delta_T:>5.1f}"
            )
    else:
        print("分析失败，无结果输出。")


if __name__ == "__main__":
    main()