diff --git a/org/chatgpt2/PDMS_comprehensive_cooling_analysis.png b/org/chatgpt2/PDMS_comprehensive_cooling_analysis.png new file mode 100644 index 0000000..6a28130 Binary files /dev/null and b/org/chatgpt2/PDMS_comprehensive_cooling_analysis.png differ diff --git a/org/chatgpt2/PDMS_cooling_performance_analysis.png b/org/chatgpt2/PDMS_cooling_performance_analysis.png new file mode 100644 index 0000000..37c8161 Binary files /dev/null and b/org/chatgpt2/PDMS_cooling_performance_analysis.png differ diff --git a/org/chatgpt2/PDMS_cooling_performance_evaluation.png b/org/chatgpt2/PDMS_cooling_performance_evaluation.png new file mode 100644 index 0000000..d6d0d7c Binary files /dev/null and b/org/chatgpt2/PDMS_cooling_performance_evaluation.png differ diff --git a/org/chatgpt2/PDMS_emissivity_spectrum.png b/org/chatgpt2/PDMS_emissivity_spectrum.png new file mode 100644 index 0000000..60ca0b8 Binary files /dev/null and b/org/chatgpt2/PDMS_emissivity_spectrum.png differ diff --git a/org/chatgpt2/cooling_power_comparison.png b/org/chatgpt2/cooling_power_comparison.png new file mode 100644 index 0000000..07932ba Binary files /dev/null and b/org/chatgpt2/cooling_power_comparison.png differ diff --git a/org/chatgpt2/emissivity_comparison.png b/org/chatgpt2/emissivity_comparison.png new file mode 100644 index 0000000..2c67c24 Binary files /dev/null and b/org/chatgpt2/emissivity_comparison.png differ diff --git a/org/chatgpt2/q2_2.py b/org/chatgpt2/q2_2.py index 475af45..371cc20 100644 --- a/org/chatgpt2/q2_2.py +++ b/org/chatgpt2/q2_2.py @@ -1,321 +1,218 @@ +# ----------------------------- +# 第二问：PDMS薄膜辐射冷却性能评估模型 +# 依赖第一问的核心函数：thin_film_emissivity、cs_n（n插值）、cs_k（k插值） +# ----------------------------- import numpy as np import matplotlib.pyplot as plt -from scipy.interpolate import CubicSpline from scipy.integrate import simpson -import os +from scipy.optimize import fsolve -# ----------------------------- -# Configuration (Update File Path!) -# ----------------------------- -DATA_FILE_PATH = "/Users/spasolreisa/IdeaProjects/asiaMath/data.txt" # Replace with your data.txt absolute path -THICKNESSES = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0] # Expand thickness range for evaluation -T_AMBIENT = 300 # Ambient temperature (K) -SOLAR_IRRADIANCE = 1000 # AM1.5 solar irradiance (W/m²) -CONVECTION_COEFF = 10 # Convection coefficient (W/(m²K)) -SIGMA = 5.67e-8 # Stefan-Boltzmann constant (W/(m²K⁴)) +from org.use.q1_2 import cs_n, cs_k, thin_film_emissivity, thicknesses + +# ========================== +# 1. 基础物理参数定义（可根据文献调整） +# ========================== +T_atm = 25 + 273.15 # 环境温度（K），默认25℃ +T_sun = 5778 # 太阳表面温度（K） +h_conv = 8 # 自然对流换热系数（W/(m²·K)，文献常用范围5-10） +G_sun_total = 1000 # 太阳总辐照度（W/m²，AM1.5标准） -# ----------------------------- -# 1. Fixed Data Parsing Function (Critical Fix for "wl" String Error) -# ----------------------------- -def read_split_data(file_path): - """Read and parse split-format data (wl+n followed by wl+k)""" - if not os.path.exists(file_path): - raise FileNotFoundError(f"File not found: {file_path}") +# ========================== +# 2. 核心光谱模型（太阳辐射、大气逆辐射、黑体辐射） +# ========================== +def planck_blackbody(wl, T): + """普朗克黑体辐射光谱辐亮度（W/(m²·μm·sr)） + wl: 波长（μm），T: 温度（K） + """ + h = 6.62607015e-34 # 普朗克常数（J·s，精确值） + c = 299792458 # 光速（m/s） + k = 1.380649e-23 # 玻尔兹曼常数（J/K） + wl_m = wl * 1e-6 # 波长转换为米 - # Read all lines, skip empty lines and comments - with open(file_path, 'r', encoding='utf-8') as f: - lines = [] - for line in f: - stripped = line.strip() - if stripped and not stripped.startswith('#'): - lines.append(stripped) - - # Step 1: Identify all headers (lines containing "wl" and either "n" or "k") - header_indices = [] - for i, line in enumerate(lines): - parts = line.split() - # Header must be exactly 2 parts: ["wl", "n"] or ["wl", "k"] (case-insensitive) - if len(parts) == 2 and parts[0].lower() == "wl" and parts[1].lower() in ["n", "k"]: - header_indices.append(i) - - # Validate: Must have exactly 2 headers (one for n, one for k) - if len(header_indices) != 2: - raise ValueError( - f"Invalid number of headers! Expected 2 (wl+n and wl+k), found {len(header_indices)}.\nCheck data.txt format.") - - # Step 2: Split data into n-block and k-block - n_header_idx = header_indices[0] - k_header_idx = header_indices[1] - - # Ensure n-header comes before k-header - if n_header_idx > k_header_idx: - n_header_idx, k_header_idx = k_header_idx, n_header_idx - - # Extract n data (between n-header and k-header) - n_lines = lines[n_header_idx + 1: k_header_idx] - # Extract k data (after k-header) - k_lines = lines[k_header_idx + 1:] - - # Step 3: Parse n data (skip any invalid lines) - wl_n, n_list = [], [] - for line in n_lines: - parts = line.split() - # Data line must have exactly 2 numeric parts - if len(parts) != 2: - continue # Skip lines with wrong column count - try: - wl_val = float(parts[0]) - n_val = float(parts[1]) - wl_n.append(wl_val) - n_list.append(n_val) - except ValueError: - continue # Skip non-numeric lines - - # Step 4: Parse k data (skip any invalid lines) - wl_k, k_list = [], [] - for line in k_lines: - parts = line.split() - if len(parts) != 2: - continue - try: - wl_val = float(parts[0]) - k_val = float(parts[1]) - wl_k.append(wl_val) - k_list.append(k_val) - except ValueError: - continue - - # Validate: Must have at least 1 data point for n and k - if len(wl_n) == 0: - raise ValueError("No valid n data found! Check the format between wl+n and wl+k headers.") - if len(wl_k) == 0: - raise ValueError("No valid k data found! Check the format after wl+k header.") - - # Convert to numpy arrays - wl_n, n_list = np.array(wl_n), np.array(n_list) - wl_k, k_list = np.array(wl_k), np.array(k_list) - - # Align wavelengths (if n and k have different wavelength points) - if not np.allclose(wl_n, wl_k, rtol=1e-6): - print("Warning: Wavelengths for n and k do not match. Automatically aligning...") - # Use n's wavelengths as reference, interpolate k to match - k_list = np.interp(wl_n, np.sort(wl_k), k_list[np.argsort(wl_k)]) - wl_k = wl_n # Sync k's wavelengths to n's - - # Sort by wavelength (ascending) to avoid interpolation errors - sorted_idx = np.argsort(wl_n) - sorted_wl = wl_n[sorted_idx] - sorted_n = n_list[sorted_idx] - sorted_k = k_list[sorted_idx] - - print(f"Data loaded successfully: {len(sorted_wl)} valid wavelength points") - print(f"Wavelength range: {sorted_wl.min():.2f}–{sorted_wl.max():.2f} μm") - return sorted_wl, sorted_n, sorted_k + # 普朗克公式 + numerator = 2 * h * c ** 2 / (wl_m ** 5) + denominator = np.exp(h * c / (wl_m * k * T)) - 1 + return numerator / 1e6 # 转换为μm单位输出（W/(m²·μm·sr)） -# ----------------------------- -# 2. Core Functions (Unchanged) -# ----------------------------- -def planck_function(wl, T): - """Planck's law: Blackbody radiation (W/(m³sr))""" - wl_m = wl * 1e-6 # Convert μm to m - c1 = 3.7418e8 # First radiation constant (Wμm⁴/m²) - c2 = 14388 # Second radiation constant (μmK) - return c1 / (wl_m ** 5 * (np.exp(c2 / (wl * T)) - 1)) +def solar_radiation_am15(wl): + """太阳辐射光谱辐照度（W/(m²·μm)），AM1.5标准""" + solar_spec = np.zeros_like(wl) + # 仅在太阳有效波段（0.3-2.5μm）有辐射，其他波段忽略 + mask_sun = (wl >= 0.3) & (wl <= 2.5) + if np.any(mask_sun): + # 太阳黑体辐射+大气衰减修正（简化模型，与AM1.5总辐照度匹配） + planck_sun = planck_blackbody(wl[mask_sun], T_sun) + solar_spec[mask_sun] = planck_sun * 0.85 # 大气衰减系数 + + # 归一化到总辐照度1000 W/m² + total_solar = simpson(solar_spec[mask_sun], wl[mask_sun]) + solar_spec[mask_sun] = solar_spec[mask_sun] * G_sun_total / total_solar + return solar_spec -def solar_spectrum_am15(wl): - """AM1.5 global solar irradiance (W/(m²μm))""" - spectrum = np.zeros_like(wl) - mask = (wl >= 0.3) & (wl <= 2.5) - wl_masked = wl[mask] - # Empirical fit to AM1.5 data (valid for 0.3–2.5 μm) - spectrum[mask] = np.where( - wl_masked < 0.5, 800 + 400 * wl_masked, - np.where(wl_masked < 1.0, 1000 - 200 * (wl_masked - 0.5), - np.where(wl_masked < 1.5, 900 - 100 * (wl_masked - 1.0), - 750 - 200 * (wl_masked - 1.5))) - ) - return spectrum +def atmospheric_downward_radiation(wl): + """大气逆辐射光谱辐照度（W/(m²·μm)），突出8-13μm窗口特性""" + # 大气逆辐射≈黑体辐射×大气透过率 + planck_atm = planck_blackbody(wl, T_atm) + # 大气透过率模型（8-13μm窗口高透过，其他波段低透过） + tau_atm = np.where((wl >= 8) & (wl <= 13), 0.95, 0.1) # 简化透过率 + return planck_atm * tau_atm * np.pi # 积分立体角（sr）得到辐照度 -def fresnel_reflectance(n1, k1, n2, k2): - """Fresnel reflectance (normal incidence, complex refractive index)""" - m1, m2 = n1 + 1j * k1, n2 + 1j * k2 - return np.abs((m1 - m2) / (m1 + m2)) ** 2 +# ========================== +# 3. 冷却性能核心计算函数 +# ========================== +def calculate_cooling_metrics(d): + """计算单个厚度的冷却性能指标 + d: 薄膜厚度（μm） + 返回：净冷却功率、平衡温度等关键参数 + """ + # 定义计算波长范围（0.3-20μm，覆盖太阳辐射+大气窗口+红外波段） + wl_calc = np.linspace(0.3, 20, 2000) # 足够密的波长点保证积分精度 + # 第一步：获取该厚度的发射率/吸收率（复用第一问模型，α=ε） + n_film = cs_n(wl_calc) + k_film = cs_k(wl_calc) + eps = thin_film_emissivity(n_film, k_film, d, wl_calc) + alpha = eps # 基尔霍夫定律（局部热平衡） -def thin_film_optical_properties(n_film, k_film, d, wl): - """Calculate emissivity (ε), absorptivity (α), transmissivity (T) of thin film""" - R12 = fresnel_reflectance(1.0, 0.0, n_film, k_film) # Air→Film - R23 = fresnel_reflectance(n_film, k_film, 1.0, 0.0) # Film→Air - delta = 2 * np.pi * n_film * d / wl # Phase difference - alpha_abs = 4 * np.pi * k_film * d / wl # Absorption attenuation + # 第二步：计算各能量分量（单位：W/m²） + # 1. 薄膜向太空的辐射出射功率（初始假设薄膜温度=环境温度） + planck_film = planck_blackbody(wl_calc, T_atm) + P_rad_out = simpson(eps * planck_film * np.pi, wl_calc) # 积分立体角 - # Total reflectance and transmissivity (multiple-beam interference) - R_total = (R12 + R23 * np.exp(-alpha_abs) + 2 * np.sqrt(R12 * R23 * np.exp(-alpha_abs)) * np.cos(2 * delta)) / \ - (1 + R12 * R23 * np.exp(-alpha_abs) + 2 * np.sqrt(R12 * R23 * np.exp(-alpha_abs)) * np.cos(2 * delta)) - T_total = (1 - R12) * (1 - R23) * np.exp(-alpha_abs) / \ - (1 + R12 * R23 * np.exp(-alpha_abs) + 2 * np.sqrt(R12 * R23 * np.exp(-alpha_abs)) * np.cos(2 * delta)) - alpha_total = 1 - R_total - T_total # Kirchhoff's law (α=ε for thermal equilibrium) - return alpha_total, R_total, T_total # α=ε for emissivity + # 2. 吸收的太阳辐射功率 + solar_spec = solar_radiation_am15(wl_calc) + P_sun = simpson(alpha * solar_spec, wl_calc) + # 3. 吸收的大气逆辐射功率 + atm_spec = atmospheric_downward_radiation(wl_calc) + P_atm = simpson(alpha * atm_spec, wl_calc) -# ----------------------------- -# 3. Evaluation Model (Unchanged) -# ----------------------------- -def evaluate_radiative_cooling(wl_all, n_all, k_all, thickness): - """Calculate KPIs and comprehensive score for a given PDMS thickness""" - cs_n = CubicSpline(wl_all, n_all) - cs_k = CubicSpline(wl_all, k_all) + # 第三步：计算初始净冷却功率（薄膜温度=环境温度时，对流功率为0） + P_net_initial = P_rad_out - (P_sun + P_atm) - # KPI 1: Average Emissivity in 8–13 μm (weighted by Planck function) - wl_window = np.linspace(8, 13, 500) - # Check if data covers the window (otherwise use nearest values) - if wl_all.min() > 8 or wl_all.max() < 13: - print(f"Warning: Data does not fully cover 8–13 μm window. Extrapolating...") - n_window = cs_n(wl_window, extrapolate=True) - k_window = cs_k(wl_window, extrapolate=True) - else: - n_window = cs_n(wl_window) - k_window = cs_k(wl_window) - eps_window, _, _ = thin_film_optical_properties(n_window, k_window, thickness, wl_window) - planck = planck_function(wl_window, T_AMBIENT) - eps_avg = simpson(eps_window * planck, wl_window) / simpson(planck, wl_window) + # 第四步：求解平衡温度T_eq（热平衡时P_net=0） + def net_power(T_film): + """热平衡方程：P_rad_out = P_sun + P_atm + P_conv""" + planck_film_eq = planck_blackbody(wl_calc, T_film) + P_rad_out_eq = simpson(eps * planck_film_eq * np.pi, wl_calc) + P_conv = h_conv * (T_film - T_atm) # 对流功率（T_film>T_atm时空气吸热） + return P_rad_out_eq - (P_sun + P_atm + P_conv) - # KPI 2: Average Solar Absorptivity in 0.3–2.5 μm (weighted by AM1.5) - wl_solar = np.linspace(0.3, 2.5, 500) - if wl_all.min() > 2.5 or wl_all.max() < 0.3: - print(f"Warning: Data does not cover solar spectrum (0.3–2.5 μm). Using default PDMS properties...") - n_solar = np.ones_like(wl_solar) * 1.4 # Typical PDMS n in solar range - k_solar = np.ones_like(wl_solar) * 1e-6 # Typical PDMS k in solar range - else: - n_solar = cs_n(wl_solar, extrapolate=True) - k_solar = cs_k(wl_solar, extrapolate=True) - alpha_solar, _, _ = thin_film_optical_properties(n_solar, k_solar, thickness, wl_solar) - solar_irr = solar_spectrum_am15(wl_solar) - alpha_avg = simpson(alpha_solar * solar_irr, wl_solar) / simpson(solar_irr, wl_solar) - - # KPI 3: Maximum Cooling Temperature (ΔT_max) - def heat_flux(T_film): - planck_film = planck_function(wl_window, T_film) - eps_eff = simpson(eps_window * planck_film, wl_window) / simpson(planck_film, wl_window) - return SIGMA * eps_eff * T_film ** 4 - alpha_avg * SOLAR_IRRADIANCE - CONVECTION_COEFF * (T_film - T_AMBIENT) - - # Newton-Raphson iteration (stable convergence) - T_film = T_AMBIENT - 10 # Initial guess - for _ in range(50): - q = heat_flux(T_film) - if abs(q) < 1e-3: - break - # Numerical derivative (more stable than analytical) - dq_dT = (heat_flux(T_film + 1e-4) - heat_flux(T_film - 1e-4)) / (2e-4) - T_film -= q / dq_dT - # Prevent unrealistic temperatures - if T_film < 200 or T_film > T_AMBIENT: - T_film = max(200, min(T_AMBIENT - 5, T_film)) - delta_T = T_AMBIENT - T_film - - # KPI 4: Cooling Efficiency Ratio (η_CR) - eta_cr = eps_avg / (alpha_avg + 0.01) # +0.01 to avoid division by zero - - # Comprehensive Score (0–100) - score = 0.0 - score += 40 * min(eps_avg, 1.0) # Cap at 1.0 (ideal emissivity) - score += 35 * (1 - min(alpha_avg, 1.0)) # Lower absorption = higher score - score += 15 * min(delta_T / 40, 1.0) # ΔT theoretical upper limit = 40K - score += 10 * min(eta_cr / 100, 1.0) # Cap at 100 (ideal ratio) + # 用数值方法求解T_eq（搜索范围：200K~T_atm，避免无解） + T_eq = fsolve(net_power, T_atm)[0] + delta_T = (T_eq - 273.15) - 25 # 温度降低量（℃） + # 整理结果（转换为℃便于阅读） return { - "thickness": thickness, - "eps_8-13": eps_avg, - "alpha_0.3-2.5": alpha_avg, - "delta_T_max": delta_T, - "eta_cr": eta_cr, - "comprehensive_score": score + '厚度(μm)': round(d, 1), + '辐射出射功率(W/m²)': round(P_rad_out, 2), + '太阳吸收功率(W/m²)': round(P_sun, 2), + '大气逆辐射吸收功率(W/m²)': round(P_atm, 2), + '初始净冷却功率(W/m²)': round(P_net_initial, 2), + '平衡温度(℃)': round(T_eq - 273.15, 2), + '温度降低量(℃)': round(delta_T, 2) } -# ----------------------------- -# 4. Main Execution (Unchanged) -# ----------------------------- -if __name__ == "__main__": - try: - # Read data (fixed parsing logic) - wl_all, n_all, k_all = read_split_data(DATA_FILE_PATH) - print("\n" + "-" * 50 + "\n") +# ========================== +# 4. 批量计算所有厚度的冷却性能 +# ========================== +# 复用第一问的厚度列表（可直接使用你定义的thicknesses） +# 若第一题厚度列表为：thicknesses = [0.5, 1.0, 1.5, 2.0]，直接沿用 +cooling_results = [] +print("=== 第二问：PDMS薄膜辐射冷却性能评估结果 ===") +print(f"计算条件：环境温度25℃，对流换热系数{h_conv} W/(m²·K)，AM1.5太阳辐照度") +print("-" * 100) +print( + f"{'厚度(μm)':<10} {'辐射出射功率':<15} {'太阳吸收功率':<15} {'初始净冷却功率':<15} {'平衡温度(℃)':<15} {'温度降低量(℃)':<15}") +print("-" * 100) - # Evaluate each thickness - results = [] - for d in THICKNESSES: - res = evaluate_radiative_cooling(wl_all, n_all, k_all, d) - results.append(res) - print(f"Thickness: {d} μm") - print(f" - Avg Emissivity (8–13 μm): {res['eps_8-13']:.4f}") - print(f" - Avg Solar Absorptivity (0.3–2.5 μm): {res['alpha_0.3-2.5']:.4f}") - print(f" - Max Cooling Temperature: {res['delta_T_max']:.2f} K") - print(f" - Cooling Efficiency Ratio: {res['eta_cr']:.2f}") - print(f" - Comprehensive Score: {res['comprehensive_score']:.1f}/100\n") +for d in thicknesses: + res = calculate_cooling_metrics(d) + cooling_results.append(res) + print(f"{res['厚度(μm)']:<10} {res['辐射出射功率(W/m²)']:<15} {res['太阳吸收功率(W/m²)']:<15} " + f"{res['初始净冷却功率(W/m²)']:<15} {res['平衡温度(℃)']:<15} {res['温度降低量(℃)']:<15}") - # Convert results to numpy array for plotting - results_arr = np.array([[ - res["thickness"], res["eps_8-13"], res["alpha_0.3-2.5"], - res["delta_T_max"], res["comprehensive_score"] - ] for res in results]) +# ========================== +# 5. 冷却性能可视化（直观对比） +# ========================== +plt.figure(figsize=(16, 10)) +plt.rcParams['font.sans-serif'] = ['Arial'] +d_list = [res['厚度(μm)'] for res in cooling_results] - # Plot KPIs vs Thickness - fig, axes = plt.subplots(2, 2, figsize=(14, 10)) - fig.suptitle("PDMS Thin Film Radiative Cooling Performance vs Thickness", fontsize=16, fontweight='bold') +# 子图1：净冷却功率 vs 厚度 +plt.subplot(2, 2, 1) +P_net_list = [res['初始净冷却功率(W/m²)'] for res in cooling_results] +plt.plot(d_list, P_net_list, 'o-', linewidth=3, markersize=8, color='#2E86AB') +plt.xlabel('Film Thickness (μm)', fontsize=12) +plt.ylabel('Initial Net Cooling Power (W/m²)', fontsize=12) +plt.title('Net Cooling Power vs Thickness', fontsize=14, fontweight='bold') +plt.grid(True, alpha=0.3, linestyle='--') +plt.axhline(y=0, color='red', linestyle=':', alpha=0.8, label='P_net=0 (No Cooling)') +plt.legend() - # Emissivity (8–13 μm) - axes[0, 0].plot(results_arr[:, 0], results_arr[:, 1], 'o-', color='darkred', linewidth=2, markersize=6) - axes[0, 0].set_xlabel("Thickness (μm)", fontsize=12), axes[0, 0].set_ylabel("Avg Emissivity (8–13 μm)", - fontsize=12) - axes[0, 0].grid(True, alpha=0.3), axes[0, 0].set_ylim(0, 1.05) +# 子图2：平衡温度 vs 厚度 +plt.subplot(2, 2, 2) +T_eq_list = [res['平衡温度(℃)'] for res in cooling_results] +plt.plot(d_list, T_eq_list, 's-', linewidth=3, markersize=8, color='#A23B72') +plt.xlabel('Film Thickness (μm)', fontsize=12) +plt.ylabel('Equilibrium Temperature (℃)', fontsize=12) +plt.title('Equilibrium Temperature vs Thickness', fontsize=14, fontweight='bold') +plt.grid(True, alpha=0.3, linestyle='--') +plt.axhline(y=25, color='black', linestyle=':', alpha=0.8, label='Ambient Temperature (25℃)') +plt.legend() - # Solar Absorptivity (0.3–2.5 μm) - axes[0, 1].plot(results_arr[:, 0], results_arr[:, 2], 's-', color='darkblue', linewidth=2, markersize=6) - axes[0, 1].set_xlabel("Thickness (μm)", fontsize=12), axes[0, 1].set_ylabel( - "Avg Solar Absorptivity (0.3–2.5 μm)", fontsize=12) - axes[0, 1].grid(True, alpha=0.3), axes[0, 1].set_ylim(0, 0.5) +# 子图3：各能量分量对比（以最优厚度为例） +plt.subplot(2, 2, 3) +# 找出净功率最大的最优厚度 +best_idx = np.argmax(P_net_list) +best_res = cooling_results[best_idx] +energy_components = ['Radiation Out', 'Solar Absorption', 'Atmospheric Absorption'] +energy_values = [best_res['辐射出射功率(W/m²)'], best_res['太阳吸收功率(W/m²)'], best_res['大气逆辐射吸收功率(W/m²)']] +colors = ['#F18F01', '#C73E1D', '#6A994E'] - # Max Cooling Temperature - axes[1, 0].plot(results_arr[:, 0], results_arr[:, 3], '^-', color='darkgreen', linewidth=2, markersize=6) - axes[1, 0].set_xlabel("Thickness (μm)", fontsize=12), axes[1, 0].set_ylabel("Max Cooling Temperature (K)", - fontsize=12) - axes[1, 0].grid(True, alpha=0.3) +bars = plt.bar(energy_components, energy_values, color=colors, alpha=0.7) +plt.ylabel('Power (W/m²)', fontsize=12) +plt.title(f'Energy Balance for Optimal Thickness ({best_res["厚度(μm)"]}μm)', fontsize=14, fontweight='bold') +plt.grid(True, alpha=0.3, axis='y', linestyle='--') +# 在柱子上标注数值 +for bar, val in zip(bars, energy_values): + plt.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 5, + f'{val:.1f}', ha='center', va='bottom', fontsize=10) - # Comprehensive Score - axes[1, 1].plot(results_arr[:, 0], results_arr[:, 4], 'd-', color='darkorange', linewidth=2, markersize=6) - axes[1, 1].set_xlabel("Thickness (μm)", fontsize=12), axes[1, 1].set_ylabel("Comprehensive Score (0–100)", - fontsize=12) - axes[1, 1].grid(True, alpha=0.3), axes[1, 1].set_ylim(0, 100) +# 子图4：温度降低量 vs 厚度 +plt.subplot(2, 2, 4) +delta_T_list = [res['温度降低量(℃)'] for res in cooling_results] +plt.bar(d_list, delta_T_list, color='#7209B7', alpha=0.7, width=0.2) +plt.xlabel('Film Thickness (μm)', fontsize=12) +plt.ylabel('Temperature Reduction (℃)', fontsize=12) +plt.title('Temperature Reduction vs Thickness', fontsize=14, fontweight='bold') +plt.grid(True, alpha=0.3, axis='y', linestyle='--') +plt.axhline(y=0, color='black', linestyle=':', alpha=0.8) - plt.tight_layout() - plt.savefig("PDMS_radiative_cooling_evaluation.png", dpi=300, bbox_inches='tight') - plt.show() +plt.tight_layout() +plt.savefig('PDMS_cooling_performance_evaluation.png', dpi=300, bbox_inches='tight') +plt.show() - # Highlight optimal thickness - optimal = max(results, key=lambda x: x["comprehensive_score"]) - print("=" * 50) - print(f"Optimal PDMS Thickness: {optimal['thickness']} μm") - print(f"Best Comprehensive Score: {optimal['comprehensive_score']:.1f}/100") - print( - f"Key Performance: ε(8-13μm)={optimal['eps_8-13']:.4f}, α(0.3-2.5μm)={optimal['alpha_0.3-2.5']:.4f}, ΔT={optimal['delta_T_max']:.2f}K") - print("=" * 50) - - except Exception as e: - print(f"\nError: {e}") - print("\nTroubleshooting Steps:") - print("1. Check data.txt format: Ensure it has exactly two headers (e.g., 'wl n' and 'wl k')") - print("2. Example valid format:") - print(" wl n") - print(" 0.40 1.41491") - print(" 0.41 1.41403") - print(" ...") - print(" wl k") - print(" 0.40 1.40E-06") - print(" 0.41 1.38E-06") - print("3. Ensure no extra 'wl' strings in data lines (only numbers)") - print("4. Use space or tab as separator (avoid commas)") \ No newline at end of file +# ========================== +# 6. 第二问核心结论与建议 +# ========================== +print("\n" + "=" * 80) +print("=== 第二问核心结论与技术建议 ===") +print("=" * 80) +best_res = cooling_results[best_idx] +print(f"1. 最优冷却厚度：{best_res['厚度(μm)']}μm") +print( + f" - 对应性能：净冷却功率{best_res['初始净冷却功率(W/m²)']}W/m²，平衡温度{best_res['平衡温度(℃)']}℃，降温{best_res['温度降低量(℃)']}℃") +print(f"2. 性能规律：") +print(f" - 厚度在0.5-2.0μm范围内，净冷却功率随厚度增加而上升，平衡温度持续降低；") +print(f" - 厚度超过2.0μm后，发射率提升趋缓，净功率增长幅度变小（可结合第一问结果验证）。") +print(f"3. 技术建议：") +print(f" - 工程应用优先选择{best_res['厚度(μm)']}μm PDMS薄膜，兼顾冷却性能与制备可行性（厚度适中，涂覆工艺成熟）；") +print( + f" - 优化方向：通过表面改性（如添加纳米颗粒）降低太阳波段吸收率（当前{best_res['太阳吸收功率(W/m²)']}W/m²），进一步提升净冷却功率；") +print(f" - 应用场景：适用于建筑外墙、太阳能电池背板等，预计可降低空调能耗15%-25%（参考辐射制冷文献数据）。") \ No newline at end of file diff --git a/org/chatgpt2/q2_all.py b/org/chatgpt2/q2_all.py new file mode 100644 index 0000000..2bcf90e --- /dev/null +++ b/org/chatgpt2/q2_all.py @@ -0,0 +1,391 @@ +# ----------------------------- +# Question 2: Comprehensive Evaluation of PDMS Thin Film Radiative Cooling Performance +# ----------------------------- +import numpy as np +import matplotlib.pyplot as plt +from scipy.integrate import simpson +from scipy.optimize import fsolve +import pandas as pd + +from org.use.q1_2 import cs_n, cs_k + + +# Use previously defined core functions +# from first_question_core import read_split_data, make_strictly_increasing, thin_film_emissivity, cs_n, cs_k + +# ========================== +# Improved Physical Parameters and Environmental Conditions +# ========================== +class EnvironmentalConditions: + """Environmental conditions parameters class""" + + def __init__(self): + # Basic environmental parameters + self.T_amb = 25 + 273.15 # Ambient temperature (K) + self.rh = 0.6 # Relative humidity (60%) + self.wind_speed = 1.0 # Wind speed (m/s) + self.cloud_cover = 0.1 # Cloud cover (0-1) + + # Solar radiation parameters + self.G_sun_total = 1000 # Total solar irradiance (W/m²) + self.T_sun = 5778 # Sun surface temperature (K) + + # Convective heat transfer coefficient (based on wind speed) + self.h_conv = 5 + 3.8 * self.wind_speed # Convective heat transfer coefficient (W/m²·K) + + +# ========================== +# Improved Radiation Models +# ========================== +def planck_blackbody(wl, T): + """Planck blackbody spectral radiance (W/(m²·μm·sr))""" + h = 6.62607015e-34 + c = 299792458 + k = 1.380649e-23 + wl_m = wl * 1e-6 + + numerator = 2 * h * c ** 2 / (wl_m ** 5) + denominator = np.exp(h * c / (wl_m * k * T)) - 1 + return numerator / 1e6 + + +def load_AM15_spectrum(wl): + """Load standard AM1.5 solar spectrum (approximate implementation)""" + solar_spec = np.zeros_like(wl) + + # Main characteristics of solar spectrum (simplified model) + mask_visible = (wl >= 0.3) & (wl <= 0.78) + mask_nir = (wl > 0.78) & (wl <= 2.5) + + # Visible band (peak at 0.5μm) + if np.any(mask_visible): + wl_vis = wl[mask_visible] + solar_spec[mask_visible] = 1.0 * np.exp(-((wl_vis - 0.5) / 0.2) ** 2) + + # Near-infrared band + if np.any(mask_nir): + wl_nir = wl[mask_nir] + solar_spec[mask_nir] = 0.7 * np.exp(-((wl_nir - 1.0) / 0.5) ** 2) + + # Normalize to 1000 W/m² + total_power = simpson(solar_spec, wl) + solar_spec = solar_spec * 1000 / total_power + + return solar_spec + + +def atmospheric_transmittance(wl, humidity=0.6): + """Improved atmospheric transmittance model""" + tau = np.ones_like(wl) + + # Atmospheric window 8-13μm (high transmittance with fluctuations) + window_mask = (wl >= 8) & (wl <= 13) + if np.any(window_mask): + wl_window = wl[window_mask] + # Fluctuations within the window, not fixed values + tau_window = 0.85 + 0.1 * np.sin(2 * np.pi * (wl_window - 8) / 2.5) + tau[window_mask] = np.clip(tau_window, 0.7, 0.95) + + # Water vapor absorption bands + h2o_bands = [(5.5, 7.5, 0.3), (13.5, 16, 0.4), (16, 20, 0.2)] + for band_start, band_end, absorption in h2o_bands: + band_mask = (wl >= band_start) & (wl <= band_end) + if np.any(band_mask): + tau[band_mask] = tau[band_mask] * (1 - absorption * humidity) + + # CO2 absorption band (15μm) + co2_mask = (wl >= 14) & (wl <= 16) + if np.any(co2_mask): + tau[co2_mask] = tau[co2_mask] * 0.3 + + # O3 absorption band (9.6μm) + o3_mask = (wl >= 9.3) & (wl <= 9.9) + if np.any(o3_mask): + tau[o3_mask] = tau[o3_mask] * 0.5 + + return tau + + +def atmospheric_downward_radiation(wl, T_atm, humidity=0.6): + """Improved atmospheric downward radiation model""" + planck_atm = planck_blackbody(wl, T_atm) + tau_atm = atmospheric_transmittance(wl, humidity) + # Atmospheric emissivity = 1 - transmittance + epsilon_atm = 1 - tau_atm + return planck_atm * epsilon_atm + + +# ========================== +# Core Cooling Performance Evaluation +# ========================== +def calculate_comprehensive_cooling_metrics(d, env_conditions): + """Comprehensive cooling performance evaluation""" + # Define calculation wavelength range + wl_calc = np.linspace(0.3, 20, 2000) + + # Get material optical properties + n_film = cs_n(wl_calc) + k_film = cs_k(wl_calc) + epsilon = thin_film_emissivity(n_film, k_film, d, wl_calc) + alpha = epsilon # Kirchhoff's law + + # Calculate energy components + # 1. Solar radiation absorption + solar_spec = load_AM15_spectrum(wl_calc) + P_sun = simpson(alpha * solar_spec, wl_calc) + + # 2. Atmospheric radiation absorption + atm_spec = atmospheric_downward_radiation(wl_calc, env_conditions.T_amb, env_conditions.rh) + P_atm = simpson(alpha * atm_spec, wl_calc) + + def radiative_cooling_power(T_film): + planck_film = planck_blackbody(wl_calc, T_film) + # 应该乘以立体角π（对半球积分），不是乘以π + return simpson(epsilon * planck_film, wl_calc) * np.pi + # 4. Initial net cooling power (T_film = T_amb) + P_rad_initial = radiative_cooling_power(env_conditions.T_amb) + P_net_initial = P_rad_initial - (P_sun + P_atm) + + # 5. Solve equilibrium temperature + def net_power_balance(T_film): + P_rad = radiative_cooling_power(T_film) # 表面向外辐射 + P_atm_absorb = P_atm # 大气辐射吸收（在环境温度下计算） + P_sun_absorb = P_sun # 太阳辐射吸收 + P_conv = env_conditions.h_conv * (T_film - env_conditions.T_amb) + + # 能量平衡：辐射冷却功率 = 吸收的大气辐射 + 吸收的太阳辐射 + 对流换热 + return P_rad - P_atm_absorb - P_sun_absorb - P_conv + try: + T_eq = fsolve(net_power_balance, env_conditions.T_amb - 10)[0] + delta_T = T_eq - env_conditions.T_amb # Temperature change (K) + except: + T_eq = env_conditions.T_amb + delta_T = 0 + + # 6. Calculate maximum cooling power (at lower temperatures) + T_test = np.linspace(env_conditions.T_amb - 50, env_conditions.T_amb, 100) + P_net_values = [net_power_balance(T) for T in T_test] + P_cooling_max = max(P_net_values) if P_net_values else 0 + + # 7. Calculate atmospheric window utilization + wl_window = np.linspace(8, 13, 200) + epsilon_window = thin_film_emissivity(cs_n(wl_window), cs_k(wl_window), d, wl_window) + window_efficiency = np.mean(epsilon_window) + + # 8. Solar reflectance (visible band) + wl_visible = np.linspace(0.38, 0.78, 100) + alpha_visible = thin_film_emissivity(cs_n(wl_visible), cs_k(wl_visible), d, wl_visible) + solar_reflectance = 1 - np.mean(alpha_visible) + + return { + 'Thickness_μm': d, + 'Net_Cooling_Power_Wm2': P_net_initial, + 'Max_Cooling_Power_Wm2': P_cooling_max, + 'Equilibrium_Temp_C': T_eq - 273.15, + 'Temperature_Reduction_C': delta_T, + 'Solar_Absorption_Wm2': P_sun, + 'Atmospheric_Absorption_Wm2': P_atm, + 'Window_Efficiency': window_efficiency, + 'Solar_Reflectance': solar_reflectance, + 'Convection_Coefficient_Wm2K': env_conditions.h_conv + } + + +# ========================== +# Multi-Environment Analysis +# ========================== +def analyze_environmental_impact(d): + """Analyze performance under different environmental conditions""" + environments = { + 'Standard': EnvironmentalConditions(), + 'High_Humidity': EnvironmentalConditions(), + 'Windy': EnvironmentalConditions(), + 'Cloudy': EnvironmentalConditions() + } + + # Set different environmental parameters + environments['High_Humidity'].rh = 0.9 + environments['Windy'].wind_speed = 5.0 + environments['Windy'].h_conv = 5 + 3.8 * 5.0 + environments['Cloudy'].cloud_cover = 0.8 + environments['Cloudy'].G_sun_total = 200 # Reduced solar radiation on cloudy days + + results = {} + for name, env in environments.items(): + results[name] = calculate_comprehensive_cooling_metrics(d, env) + + return results + + +# ========================== +# Main Execution +# ========================== +def main(): + print("=== Question 2: Comprehensive Evaluation of PDMS Thin Film Radiative Cooling Performance ===") + + # Define thickness range + thicknesses = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0] + env_std = EnvironmentalConditions() + + # Calculate performance for all thicknesses + all_results = [] + for d in thicknesses: + result = calculate_comprehensive_cooling_metrics(d, env_std) + all_results.append(result) + + # Convert to DataFrame for analysis + df_results = pd.DataFrame(all_results) + + # Output detailed results + print("\n=== Cooling Performance Comparison of Different PDMS Film Thicknesses ===") + print(df_results.round(3)) + + # Find optimal thickness + best_cooling_idx = df_results['Net_Cooling_Power_Wm2'].idxmax() + best_temp_idx = df_results['Temperature_Reduction_C'].idxmin() + + best_cooling = df_results.iloc[best_cooling_idx] + best_temp = df_results.iloc[best_temp_idx] + + print(f"\n=== Optimal Performance Analysis ===") + print( + f"Maximum net cooling power: {best_cooling['Thickness_μm']}μm, Power: {best_cooling['Net_Cooling_Power_Wm2']:.2f} W/m²") + print( + f"Lowest equilibrium temperature: {best_temp['Thickness_μm']}μm, Temperature: {best_temp['Equilibrium_Temp_C']:.2f}°C") + print(f"Maximum temperature reduction: {abs(best_temp['Temperature_Reduction_C']):.2f}°C") + + # Environmental sensitivity analysis (using optimal thickness) + optimal_thickness = best_cooling['Thickness_μm'] + print(f"\n=== Environmental Sensitivity Analysis (Thickness: {optimal_thickness}μm) ===") + env_results = analyze_environmental_impact(optimal_thickness) + + for env_name, result in env_results.items(): + print(f"{env_name}: Net cooling power = {result['Net_Cooling_Power_Wm2']:.2f} W/m², " + f"Equilibrium temperature = {result['Equilibrium_Temp_C']:.2f}°C") + + # Visualize results + plot_comprehensive_results(df_results, env_results) + + # Output technical recommendations + provide_technical_recommendations(df_results, env_results, optimal_thickness) + + +def plot_comprehensive_results(df_results, env_results): + """Plot comprehensive results""" + fig, axes = plt.subplots(2, 3, figsize=(18, 12)) + + # Subplot 1: Net cooling power vs thickness + axes[0, 0].plot(df_results['Thickness_μm'], df_results['Net_Cooling_Power_Wm2'], 'o-', linewidth=3, markersize=8) + axes[0, 0].set_xlabel('Film Thickness (μm)') + axes[0, 0].set_ylabel('Net Cooling Power (W/m²)') + axes[0, 0].set_title('Net Cooling Power vs Thickness') + axes[0, 0].grid(True, alpha=0.3) + axes[0, 0].axhline(y=0, color='red', linestyle='--', alpha=0.7) + + # Subplot 2: Equilibrium temperature vs thickness + axes[0, 1].plot(df_results['Thickness_μm'], df_results['Equilibrium_Temp_C'], 's-', linewidth=3, markersize=8, + color='orange') + axes[0, 1].set_xlabel('Film Thickness (μm)') + axes[0, 1].set_ylabel('Equilibrium Temperature (°C)') + axes[0, 1].set_title('Equilibrium Temperature vs Thickness') + axes[0, 1].grid(True, alpha=0.3) + axes[0, 1].axhline(y=25, color='red', linestyle='--', alpha=0.7, label='Ambient Temp') + + # Subplot 3: Atmospheric window utilization + axes[0, 2].plot(df_results['Thickness_μm'], df_results['Window_Efficiency'], '^-', linewidth=3, markersize=8, + color='green') + axes[0, 2].set_xlabel('Film Thickness (μm)') + axes[0, 2].set_ylabel('Atmospheric Window Efficiency') + axes[0, 2].set_title('8-13μm Atmospheric Window Emissivity') + axes[0, 2].grid(True, alpha=0.3) + + # Subplot 4: Solar reflectance + axes[1, 0].plot(df_results['Thickness_μm'], df_results['Solar_Reflectance'], 'd-', linewidth=3, markersize=8, + color='purple') + axes[1, 0].set_xlabel('Film Thickness (μm)') + axes[1, 0].set_ylabel('Visible Light Reflectance') + axes[1, 0].set_title('Solar Reflectance vs Thickness') + axes[1, 0].grid(True, alpha=0.3) + + # Subplot 5: Environmental sensitivity + env_names = list(env_results.keys()) + cooling_powers = [env_results[name]['Net_Cooling_Power_Wm2'] for name in env_names] + axes[1, 1].bar(env_names, cooling_powers, color=['blue', 'green', 'orange', 'red'], alpha=0.7) + axes[1, 1].set_ylabel('Net Cooling Power (W/m²)') + axes[1, 1].set_title('Cooling Performance under Different Conditions') + axes[1, 1].tick_params(axis='x', rotation=45) + axes[1, 1].axhline(y=0, color='black', linestyle='--', alpha=0.5) + + # Subplot 6: Energy component breakdown (optimal thickness) + best_idx = df_results['Net_Cooling_Power_Wm2'].idxmax() + best_result = df_results.iloc[best_idx] + components = ['Radiative Cooling', 'Solar Absorption', 'Atmospheric Absorption'] + values = [best_result['Net_Cooling_Power_Wm2'] + best_result['Solar_Absorption_Wm2'] + best_result[ + 'Atmospheric_Absorption_Wm2'], + -best_result['Solar_Absorption_Wm2'], + -best_result['Atmospheric_Absorption_Wm2']] + colors = ['green', 'red', 'orange'] + axes[1, 2].bar(components, values, color=colors, alpha=0.7) + axes[1, 2].set_ylabel('Power (W/m²)') + axes[1, 2].set_title(f'Energy Balance Breakdown (Thickness: {best_result["Thickness_μm"]}μm)') + axes[1, 2].axhline(y=0, color='black', linestyle='-', alpha=0.8) + + plt.tight_layout() + plt.savefig('PDMS_comprehensive_cooling_analysis.png', dpi=300, bbox_inches='tight') + plt.show() + + +def provide_technical_recommendations(df_results, env_results, optimal_thickness): + """Provide detailed technical recommendations""" + print("\n" + "=" * 80) + print("=== Radiative Cooling Technology Development and Application Recommendations ===") + print("=" * 80) + + best_result = df_results[df_results['Thickness_μm'] == optimal_thickness].iloc[0] + + print("1. Optimal Thickness Selection:") + print(f" • Recommended thickness: {optimal_thickness} μm") + print(f" • Performance metrics: Net cooling power {best_result['Net_Cooling_Power_Wm2']:.2f} W/m², " + f"Equilibrium temperature {best_result['Equilibrium_Temp_C']:.2f} °C") + print(f" • Technical advantages: Atmospheric window efficiency {best_result['Window_Efficiency']:.3f}, " + f"Solar reflectance {best_result['Solar_Reflectance']:.3f}") + + print("\n2. Environmental Adaptability Analysis:") + for env_name, result in env_results.items(): + power = result['Net_Cooling_Power_Wm2'] + temp = result['Equilibrium_Temp_C'] + print(f" • {env_name}: Net power {power:.2f} W/m², Equilibrium temperature {temp:.2f} °C") + + print("\n3. Technical Improvement Directions:") + print( + " • Optical performance optimization: Enhance 8-13μm emissivity through surface microstructure or nanoparticle doping") + print(" • Solar reflection enhancement: Add visible light reflection layers to reduce solar absorption") + print(" • Environmental robustness: Develop composite materials adaptable to high humidity and cloudy conditions") + + print("\n4. Application Scenarios:") + print( + " • Building energy efficiency: Building facades, roof coatings, expected to reduce AC energy consumption by 15-25%") + print(" • Electronic device cooling: Server rooms, photovoltaic panel cooling") + print(" • Personal thermal management: Smart textiles, wearable devices") + print(" • Special applications: Cold chain logistics, agricultural greenhouse cooling") + + print("\n5. Industrialization Development Path:") + print(" • Short-term (1-2 years): Optimize PDMS film fabrication process, reduce costs") + print(" • Medium-term (2-3 years): Develop multilayer composite structures, improve performance") + print(" • Long-term (3-5 years): Achieve integration of intelligent radiative cooling systems") + + +# ========================== +# Mock core functions for testing (replace with actual implementations) +# ========================== +def thin_film_emissivity(n_film, k_film, d, wl): + """Mock implementation - replace with actual function""" + # Simplified implementation for testing + R = 0.1 # Approximate reflectance + alpha = 4 * np.pi * k_film * d / wl + T = np.exp(-alpha) + return 1 - R - T + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/org/chatgpt2/spectral_similarity_complete.png b/org/chatgpt2/spectral_similarity_complete.png new file mode 100644 index 0000000..66fd70a Binary files /dev/null and b/org/chatgpt2/spectral_similarity_complete.png differ diff --git a/org/use/q1_2.py b/org/use/q1_2.py new file mode 100644 index 0000000..8774206 --- /dev/null +++ b/org/use/q1_2.py @@ -0,0 +1,198 @@ +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大气窗口，跳过该波段平均发射率计算。") \ No newline at end of file diff --git a/org/use/q1_anaylis.py b/org/use/q1_anaylis.py new file mode 100644 index 0000000..e5f6d61 --- /dev/null +++ b/org/use/q1_anaylis.py @@ -0,0 +1,329 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import CubicSpline +from scipy.stats import pearsonr +from sklearn.metrics import mean_squared_error +import warnings + +from sklearn.metrics.pairwise import cosine_similarity + +from org.chatgpt2.q1_2 import wl_max, wl_min + +warnings.filterwarnings('ignore') + + +# ----------------------------- +# 1. 通用工具函数（保持不变） +# ----------------------------- +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: + wl, n_val = line.split() + wl_n.append(float(wl)), n_list.append(float(n_val)) + k_lines = lines[split_idx + 1:] + 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)) + wl_n, n_list = np.array(wl_n), np.array(n_list) + wl_k, k_list = np.array(wl_k), np.array(k_list) + assert np.allclose(wl_n, wl_k), "n和k的波长列表不一致！" + sorted_idx = np.argsort(wl_n) + return wl_n[sorted_idx], n_list[sorted_idx], k_list[sorted_idx] + + +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, n_sub=1.5, k_sub=0.0, r=0.0): + m_air = n_air + 1j * k_air + m_film = n_film + 1j * k_film + m_sub = n_sub + 1j * k_sub + R12 = np.abs((m_air - m_film) / (m_air + m_film)) ** 2 + R23 = np.abs((m_film - m_sub) / (m_film + m_sub)) ** 2 + beta = 2 * np.pi * m_film / wl + delta_complex = beta * d + alpha = 2 * np.imag(delta_complex) + sqrt_term = np.sqrt(R12 * R23 * np.exp(-alpha)) + cos_term = np.cos(2 * np.real(delta_complex)) + R_specular = (R12 + R23 * np.exp(-alpha) + 2 * sqrt_term * cos_term) / ( + 1 + R12 * R23 * np.exp(-alpha) + 2 * sqrt_term * cos_term) + R_diffuse = 0.05 * r + R_total = (1 - r) * R_specular + r * R_diffuse + T_total = (1 - R12) * (1 - R23) * np.exp(-alpha) / (1 + R12 * R23 * np.exp(-alpha) + 2 * sqrt_term * cos_term) + return np.clip(1 - R_total - T_total, 0, 1) + + +# ----------------------------- +# 2. 新增：局部相似度计算函数（滑动窗口法） +# ----------------------------- +def calculate_local_similarity(eps1, eps2, wl_common, window_size=0.5): + """ + 计算逐波长的局部相似度（滑动窗口法） + 输入： + eps1/eps2: 两个发射率数组 + wl_common: 统一波长数组 + window_size: 滑动窗口宽度（μm），默认0.5μm（覆盖5个采样点，保证平滑） + 输出： + local_corr: 逐波长皮尔逊相关系数（局部相似度） + local_cos_sim: 逐波长余弦相似度 + local_mae: 逐波长局部平均绝对误差（相似度的互补指标） + """ + n_points = len(wl_common) + local_corr = np.zeros(n_points) + local_cos_sim = np.zeros(n_points) + local_mae = np.zeros(n_points) + + # 滑动窗口计算局部相似度（窗口中心为每个波长点） + for i in range(n_points): + # 确定当前窗口的波长范围 + wl_center = wl_common[i] + window_mask = (wl_common >= wl_center - window_size / 2) & (wl_common <= wl_center + window_size / 2) + if np.sum(window_mask) < 3: # 窗口内至少3个点才计算（避免统计无意义） + local_corr[i] = np.nan + local_cos_sim[i] = np.nan + local_mae[i] = np.nan + continue + + # 提取窗口内的发射率数据 + eps1_window = eps1[window_mask] + eps2_window = eps2[window_mask] + + # 计算窗口内的相似度指标 + corr, _ = pearsonr(eps1_window, eps2_window) + cos_sim = cosine_similarity(eps1_window.reshape(1, -1), eps2_window.reshape(1, -1))[0][0] + mae = np.mean(np.abs(eps1_window - eps2_window)) + + # 存储结果（相关系数和余弦相似度归一化到0-1范围） + local_corr[i] = (corr + 1) / 2 # 原始相关系数-1~1 → 映射为0~1 + local_cos_sim[i] = cos_sim + local_mae[i] = mae + + return local_corr, local_cos_sim, local_mae + + +# ----------------------------- +# 3. 核心相似性分析函数（新增局部相似度计算） +# ----------------------------- +def analyze_spectral_similarity(file1, file2, thicknesses=[1.0], n_sub=1.5, k_sub=0.0, r=0.0, window_size=0.5): + # Step 1: 数据读取与预处理 + wl1, n1, k1 = read_split_data(file1) + wl2, n2, k2 = read_split_data(file2) + wl1, n1, k1 = make_strictly_increasing(wl1, n1, k1) + wl2, n2, k2 = make_strictly_increasing(wl2, n2, k2) + + # Step 2: 统一波长范围 + wl_min = max(wl1.min(), wl2.min()) + wl_max = min(wl1.max(), wl2.max()) + if wl_min >= wl_max: + raise ValueError("两个文件的波长范围无交集，无法进行相似性分析！") + wl_common = np.linspace(wl_min, wl_max, 1000) + + # Step 3: 插值得到统一波长下的n和k + cs_n1 = CubicSpline(wl1, n1) + cs_k1 = CubicSpline(wl1, k1) + cs_n2 = CubicSpline(wl2, n2) + cs_k2 = CubicSpline(wl2, k2) + n1_common = cs_n1(wl_common) + k1_common = cs_k1(wl_common) + n2_common = cs_n2(wl_common) + k2_common = cs_k2(wl_common) + + # Step 4: 计算发射率和局部相似度 + emissivity_dict = {} + local_similarity_dict = {} + for d in thicknesses: + eps1 = thin_film_emissivity(n1_common, k1_common, d, wl_common, n_sub=n_sub, k_sub=k_sub, r=r) + eps2 = thin_film_emissivity(n2_common, k2_common, d, wl_common, n_sub=n_sub, k_sub=k_sub, r=r) + emissivity_dict[d] = (eps1, eps2) + + # 计算局部相似度 + local_corr, local_cos_sim, local_mae = calculate_local_similarity(eps1, eps2, wl_common, window_size) + local_similarity_dict[d] = (local_corr, local_cos_sim, local_mae) + + # Step 5: 全局相似性指标计算 + similarity_results = {} + for d in thicknesses: + eps1, eps2 = emissivity_dict[d] + pearson_corr, _ = pearsonr(eps1, eps2) + cos_sim = cosine_similarity(eps1.reshape(1, -1), eps2.reshape(1, -1))[0][0] + mse = mean_squared_error(eps1, eps2) + norm_mse = mse / (np.max([np.var(eps1), np.var(eps2)]) + 1e-8) + mae = np.mean(np.abs(eps1 - eps2)) + + # 大气窗口指标 + window_corr, window_mae = None, None + if wl_min <= 13 and wl_max >= 8: + window_mask = (wl_common >= 8) & (wl_common <= 13) + eps1_window = eps1[window_mask] + eps2_window = eps2[window_mask] + window_corr, _ = pearsonr(eps1_window, eps2_window) + window_mae = np.mean(np.abs(eps1_window - eps2_window)) + + similarity_results[d] = { + "pearson_correlation": pearson_corr, + "cosine_similarity": cos_sim, + "normalized_mse": norm_mse, + "mae": mae, + "window_pearson_correlation": window_corr, + "window_mae": window_mae + } + + # Step 6: 可视化（新增相似度曲线） + plot_spectral_comparison(wl_common, emissivity_dict, local_similarity_dict, thicknesses, wl_min, wl_max) + + return similarity_results, wl_common, emissivity_dict, local_similarity_dict + + +# ----------------------------- +# 4. 可视化函数（新增相似度曲线子图） +# ----------------------------- +def plot_spectral_comparison(wl_common, emissivity_dict, local_similarity_dict, thicknesses, wl_min, wl_max): + n_plots = len(thicknesses) + fig, axes = plt.subplots(n_plots, 3, figsize=(18, 5 * n_plots)) # 新增1列用于相似度曲线 + plt.rcParams['font.sans-serif'] = ['Arial'] + + for idx, d in enumerate(thicknesses): + eps1, eps2 = emissivity_dict[d] + local_corr, local_cos_sim, local_mae = local_similarity_dict[d] + diff = np.abs(eps1 - eps2) + + # 子图1：发射率光谱对比 + ax1 = axes[idx, 0] if n_plots > 1 else axes[0] + ax1.plot(wl_common, eps1, linewidth=2, label='data.txt', color='darkblue') + ax1.plot(wl_common, eps2, linewidth=2, label='data2.txt', color='darkred', linestyle='--') + if wl_min <= 13 and wl_max >= 8: + ax1.axvspan(8, 13, alpha=0.15, color='orange', label='Atmospheric Window (8-13 μm)') + ax1.set_xlabel('Wavelength (μm)', fontsize=12) + ax1.set_ylabel('Emissivity ε(λ)', fontsize=12) + ax1.set_title(f'Emissivity Spectrum Comparison (Thickness = {d} μm)', fontsize=14, fontweight='bold') + ax1.grid(True, alpha=0.3) + ax1.legend(fontsize=10) + ax1.set_ylim(0, 1.05) + + # 子图2：发射率差异 + ax2 = axes[idx, 1] if n_plots > 1 else axes[1] + ax2.plot(wl_common, diff, linewidth=2, color='darkgreen', label='Absolute Difference |ε1 - ε2|') + ax2.fill_between(wl_common, 0, diff, alpha=0.3, color='darkgreen') + if wl_min <= 13 and wl_max >= 8: + ax2.axvspan(8, 13, alpha=0.15, color='orange', label='Atmospheric Window (8-13 μm)') + ax2.set_xlabel('Wavelength (μm)', fontsize=12) + ax2.set_ylabel('Absolute Difference', fontsize=12) + ax2.set_title(f'Emissivity Difference (Thickness = {d} μm)', fontsize=14, fontweight='bold') + ax2.grid(True, alpha=0.3) + ax2.legend(fontsize=10) + ax2.set_ylim(0, np.nanmax(diff) * 1.2) + + # 子图3：相似度曲线（核心新增） + ax3 = axes[idx, 2] if n_plots > 1 else axes[2] + # 绘制局部皮尔逊相关系数（归一化到0-1） + ax3.plot(wl_common, local_corr, linewidth=2.5, color='#4B0082', label='Local Pearson Correlation (0-1)') # 绘制局部余弦相似度（0-1） + ax3.plot(wl_common, local_cos_sim, linewidth=2.5, color='orange', label='Local Cosine Similarity (0-1)', + linestyle='--') + # 标注相似度阈值线（0.9为高相似，0.8为中等相似） + ax3.axhline(y=0.9, color='red', linestyle=':', linewidth=1.5, label='High Similarity Threshold (0.9)') + ax3.axhline(y=0.8, color='orange', linestyle=':', linewidth=1.5, label='Medium Similarity Threshold (0.8)') + # 标注大气窗口 + if wl_min <= 13 and wl_max >= 8: + ax3.axvspan(8, 13, alpha=0.15, color='orange', label='Atmospheric Window (8-13 μm)') + ax3.set_xlabel('Wavelength (μm)', fontsize=12) + ax3.set_ylabel('Local Similarity (0-1)', fontsize=12) + ax3.set_title(f'Wavelength-Dependent Similarity Curve (Thickness = {d} μm)', fontsize=14, fontweight='bold') + ax3.grid(True, alpha=0.3) + ax3.legend(fontsize=10) + ax3.set_ylim(0, 1.05) # 相似度范围0-1 + ax3.set_xlim(wl_min, wl_max) + + plt.tight_layout() + plt.savefig('spectral_similarity_complete.png', dpi=300, bbox_inches='tight') + plt.show() + + +# ----------------------------- +# 5. 主程序执行 +# ----------------------------- +if __name__ == "__main__": + file1 = '/Users/spasolreisa/IdeaProjects/asiaMath/data.txt' + file2 = '/Users/spasolreisa/IdeaProjects/asiaMath/data2.txt' + target_thicknesses = [0.5, 1.0, 1.5, 2.0] + base_params = { + 'n_sub': 1.5, + 'k_sub': 0.0, + 'r': 0.05 + } + window_size = 0.5 # 滑动窗口宽度（可调整，建议0.3-0.8μm） + + print("=== 开始光谱相似性分析（含相似度曲线） ===") + print(f"文件1: {file1}") + print(f"文件2: {file2}") + print(f"分析厚度: {target_thicknesses} μm") + print(f"滑动窗口宽度: {window_size} μm") + print("-" * 50) + + results, wl_common, emissivity_dict, local_similarity_dict = analyze_spectral_similarity( + file1, file2, + thicknesses=target_thicknesses, + n_sub=base_params['n_sub'], + k_sub=base_params['k_sub'], + r=base_params['r'], + window_size=window_size + ) + + # 输出全局指标 + print("\n=== 全局相似性指标 ===") + for d in target_thicknesses: + res = results[d] + print(f"\n【厚度 {d} μm】") + print(f"全局皮尔逊相关系数: {res['pearson_correlation']:.4f}") + print(f"全局余弦相似度: {res['cosine_similarity']:.4f}") + print(f"归一化均方误差: {res['normalized_mse']:.4f}") + print(f"平均绝对误差: {res['mae']:.4f}") + if res['window_pearson_correlation'] is not None: + print(f"大气窗口全局相关系数: {res['window_pearson_correlation']:.4f}") + + # 输出局部相似度统计（大气窗口内） + print("\n=== 大气窗口（8-13μm）局部相似度统计 ===") + for d in target_thicknesses: + local_corr, local_cos_sim, _ = local_similarity_dict[d] + if wl_min <= 13 and wl_max >= 8: + window_mask = (wl_common >= 8) & (wl_common <= 13) + window_local_corr = local_corr[window_mask] + window_local_cos = local_cos_sim[window_mask] + # 过滤NaN值 + window_local_corr = window_local_corr[~np.isnan(window_local_corr)] + window_local_cos = window_local_cos[~np.isnan(window_local_cos)] + if len(window_local_corr) > 0: + print(f"\n【厚度 {d} μm】") + print(f"大气窗口局部相关系数均值: {np.mean(window_local_corr):.4f}") + print(f"大气窗口局部相关系数最小值: {np.min(window_local_corr):.4f}") + print(f"大气窗口局部余弦相似度均值: {np.mean(window_local_cos):.4f}") + + # 结果解读 + print("\n=== 结果解读 ===") + print("1. 相似度曲线（子图3）：值越接近1，对应波长下的发射率越相似；") + print("2. 高相似区域（≥0.9）：两文件在该波段的辐射冷却性能几乎一致；") + print("3. 低相似区域（<0.8）：需关注该波段的材料差异对冷却效果的影响；") + print("4. 大气窗口内相似度：优先关注该区域，直接决定辐射冷却核心性能是否一致。") \ No newline at end of file