diff --git a/org/other/__init__.py b/org/other/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/org/other/question1_pdms_emissivity.py b/org/other/question1_pdms_emissivity.py new file mode 100644 index 0000000..b56c969 --- /dev/null +++ b/org/other/question1_pdms_emissivity.py @@ -0,0 +1,151 @@ +import os +from dataclasses import dataclass +from typing import Dict, List + +import numpy as np +from matplotlib import pyplot as plt + +AIR_INDEX = 1.0 # normal incidence, non-absorbing ambient +WAVELENGTH_MIN = 0.4 # µm +WAVELENGTH_MAX = 25.0 # µm + + +def cauchy_index(wavelength_um: np.ndarray) -> np.ndarray: + """ + Dispersion model for PDMS based on Cauchy equation fitted to literature data: + n(λ) = A + B / λ^2 + C / λ^4 + The coefficients (A=1.385, B=0.0125, C=0.00045) match n≈1.43 @ 0.55 µm + and n≈1.40 beyond 2 µm. + """ + A = 1.385 + B = 0.0125 + C = 0.00045 + lam2 = wavelength_um ** 2 + return A + B / lam2 + C / (lam2 * wavelength_um ** 2) + + +def extinction_coeff(wavelength_um: np.ndarray) -> np.ndarray: + """ + Construct an approximate extinction coefficient profile by superposing + Gaussians centered at the known vibrational bands of PDMS in the mid-IR. + Peaks taken from FTIR measurements (Si-O-Si at ~9.2 µm, Si-CH3 rocking at + 12.7 µm, CH stretches at 3.4 µm, etc.). + """ + peaks = [ + # (center µm, amplitude, width µm) + (3.4, 0.06, 0.25), + (8.0, 0.30, 0.50), + (9.2, 0.65, 0.60), + (10.6, 0.35, 0.45), + (12.7, 0.45, 0.35), + (14.5, 0.25, 0.60), + ] + k = np.full_like(wavelength_um, 5e-4) + for center, amp, width in peaks: + k += amp * np.exp(-0.5 * ((wavelength_um - center) / width) ** 2) + # enforce plausible upper bound + return np.clip(k, 0, 2.0) + + +def thin_film_emissivity( + wavelength_um: np.ndarray, thickness_um: float, n_complex: np.ndarray +) -> np.ndarray: + """ + Normal-incidence emissivity for a PDMS slab in air computed using + the transfer-matrix solution for a single absorbing layer. + """ + n0 = AIR_INDEX + ns = AIR_INDEX + thickness_m = thickness_um * 1e-6 + # vacuum wavenumber + beta = 2 * np.pi / (wavelength_um * 1e-6) + delta = beta * n_complex * thickness_m + + r01 = (n0 - n_complex) / (n0 + n_complex) + r12 = (n_complex - ns) / (n_complex + ns) + exp_term = np.exp(2j * delta) + numerator_r = r01 + r12 * exp_term + denominator = 1 + r01 * r12 * exp_term + r = numerator_r / denominator + + t01 = 2 * n0 / (n0 + n_complex) + t12 = 2 * n_complex / (n_complex + ns) + exp_half = np.exp(1j * delta) + t = (t01 * t12 * exp_half) / denominator + + R = np.abs(r) ** 2 + T = (ns.real / n0) * np.abs(t) ** 2 + A = 1 - R - T + return np.clip(A.real, 0, 1) + + +@dataclass +class EmissivityResult: + wavelengths: np.ndarray + emissivity_map: Dict[float, np.ndarray] + window_avg: Dict[float, float] + + +def compute_emissivity(thicknesses_um: List[float]) -> EmissivityResult: + wavelengths = np.linspace(WAVELENGTH_MIN, WAVELENGTH_MAX, 1200) + n = cauchy_index(wavelengths) + k = extinction_coeff(wavelengths) + n_complex = n + 1j * k + + emissivity_map = {} + window_avg = {} + window_mask = (wavelengths >= 8) & (wavelengths <= 13) + + for d in thicknesses_um: + eps = thin_film_emissivity(wavelengths, d, n_complex) + emissivity_map[d] = eps + window_avg[d] = float(np.trapz(eps[window_mask], wavelengths[window_mask]) / + np.trapz(np.ones_like(wavelengths[window_mask]), + wavelengths[window_mask])) + + return EmissivityResult(wavelengths, emissivity_map, window_avg) + + +def plot_emissivity(result: EmissivityResult, outdir: str) -> str: + plt.figure(figsize=(9, 5)) + for d, eps in result.emissivity_map.items(): + label = f"{d:.0f} µm (ε̄₈₋₁₃={result.window_avg[d]:.2f})" + plt.plot(result.wavelengths, eps, label=label) + + plt.axvspan(8, 13, color="tab:gray", alpha=0.15, label="Atmospheric window") + plt.xlabel("Wavelength (µm)") + plt.ylabel("Spectral emissivity") + plt.title("PDMS Thin-Film Emissivity vs. Wavelength") + plt.ylim(0, 1.05) + plt.xlim(WAVELENGTH_MIN, WAVELENGTH_MAX) + plt.legend() + plt.grid(alpha=0.3) + + os.makedirs(outdir, exist_ok=True) + output_path = os.path.join(outdir, "question1_emissivity.png") + plt.tight_layout() + plt.savefig(output_path, dpi=300) + plt.close() + return output_path + + +def main(): + thicknesses = [1, 5, 10, 25, 50, 100] + result = compute_emissivity(thicknesses) + outdir = os.path.join(os.path.dirname(__file__), "outputs") + figure_path = plot_emissivity(result, outdir) + + summary_lines = ["thickness_um,avg_emissivity_8_13um"] + for d in thicknesses: + summary_lines.append(f"{d},{result.window_avg[d]:.4f}") + csv_path = os.path.join(outdir, "question1_emissivity_summary.csv") + with open(csv_path, "w", encoding="utf-8") as f: + f.write("\n".join(summary_lines)) + + print(f"Figure saved to: {figure_path}") + print(f"Summary saved to: {csv_path}") + + +if __name__ == "__main__": + main() + diff --git a/org/other/question2_pdms_cooling.py b/org/other/question2_pdms_cooling.py new file mode 100644 index 0000000..072f88b --- /dev/null +++ b/org/other/question2_pdms_cooling.py @@ -0,0 +1,157 @@ +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 + +plt.rcParams["font.sans-serif"] = ["DejaVu Sans", "Arial", "Helvetica"] +plt.rcParams["axes.unicode_minus"] = False + +SIGMA = 5.670374419e-8 + +THICKNESSES = [1, 5, 10, 25, 50, 100] # µm +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 load_emissivity_module(): + here = os.path.dirname(os.path.abspath(__file__)) + target = os.path.join(here, "question1_pdms_emissivity.py") + spec = importlib.util.spec_from_file_location("q1", target) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) # type: ignore[arg-type] + return module + + +def solar_absorptance(thickness_um: float) -> float: + """ + Empirical assumption: In PDMS silver-coated system, visible light absorption + is determined by thin film scattering and internal losses. + At 1 µm, α≈0.05, absorption increases by 0.01 for every 10 µm increase in thickness, capped at 0.15. + """ + 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: + # No sign change, return endpoint with smaller energy consumption + 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 evaluate() -> Tuple[List[CoolingResult], str]: + q1 = load_emissivity_module() + emissivity_data = q1.compute_emissivity(THICKNESSES) + + results: List[CoolingResult] = [] + for d in THICKNESSES: + eps_window = emissivity_data.window_avg[d] + alpha_s = solar_absorptance(d) + 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.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.png") + plot_results(results, fig_path) + + 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", 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") + fig.tight_layout() + plt.savefig(fig_path, dpi=300) + plt.close(fig) + + +def main(): + results, csv_path = evaluate() + print(f"Results written to: {csv_path}") + for r in results: + print( + f"d={r.thickness_um:>3} um, eps_8_13={r.eps_window:.2f}, alpha_sol={r.alpha_solar:.2f}, " + f"q_net(T_amb)={r.net_power_at_amb:.1f} W/m2, T_eq={r.eq_temp_C:.1f} °C " + f"(DeltaT={r.delta_T:.1f} K)" + ) + + +if __name__ == "__main__": + main() diff --git a/org/other/question3_multilayer_optimization.py b/org/other/question3_multilayer_optimization.py new file mode 100644 index 0000000..50d3586 --- /dev/null +++ b/org/other/question3_multilayer_optimization.py @@ -0,0 +1,236 @@ +import itertools +import math +import os +import random +from dataclasses import dataclass +from typing import Dict, List, Sequence, Tuple + +import numpy as np +from matplotlib import pyplot as plt + +plt.rcParams["font.sans-serif"] = ["DejaVu Sans", "Arial", "Helvetica"] +plt.rcParams["axes.unicode_minus"] = False + + +def load_q1_module(): + here = os.path.dirname(os.path.abspath(__file__)) + target = os.path.join(here, "question1_pdms_emissivity.py") + import importlib.util + + spec = importlib.util.spec_from_file_location("q1", target) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) # type: ignore[arg-type] + return module + + +def planck_weight(wavelength_um: np.ndarray, temperature: float = 300.0) -> np.ndarray: + wl_m = wavelength_um * 1e-6 + c1 = 3.7418e-16 + c2 = 1.4388e-2 + spectral = c1 / (wl_m**5 * (np.exp(c2 / (wl_m * temperature)) - 1)) + return spectral + + +def solar_weight(wavelength_um: np.ndarray) -> np.ndarray: + center1, width1 = 0.6, 0.35 + center2, width2 = 1.6, 0.45 + return np.exp(-((wavelength_um - center1) / width1) ** 2) + 0.35 * np.exp( + -((wavelength_um - center2) / width2) ** 2 + ) + + +@dataclass +class Material: + name: str + n_const: float + k_const: float + + def nk(self, wavelength_um: np.ndarray) -> np.ndarray: + n = np.full_like(wavelength_um, self.n_const, dtype=np.complex128) + k = np.full_like(wavelength_um, self.k_const, dtype=np.complex128) + return n - 1j * k + + +def pdms_index(wavelength_um: np.ndarray) -> np.ndarray: + q1 = load_q1_module() + n = q1.cauchy_index(wavelength_um) + k = q1.extinction_coeff(wavelength_um) + return n - 1j * k + + +def ag_index(wavelength_um: np.ndarray) -> np.ndarray: + n = 0.15 + 0.6 * np.exp(-((wavelength_um - 0.5) / 0.4) ** 2) + k = 4.5 + 3.5 * np.exp(-((wavelength_um - 10) / 6) ** 2) + return n - 1j * k + + +def transfer_matrix_stack( + wavelength_um: np.ndarray, + layer_nk: Sequence[np.ndarray], + thickness_um: Sequence[float], + substrate_nk: np.ndarray, + n0: float = 1.0, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + beta = 2 * np.pi / (wavelength_um * 1e-6) + q0 = n0 + qs = substrate_nk + + R = np.zeros_like(wavelength_um) + T = np.zeros_like(wavelength_um) + + for idx, wl in enumerate(wavelength_um): + M = np.identity(2, dtype=complex) + for nk, d in zip(layer_nk, thickness_um): + n_layer = nk[idx] + delta = beta[idx] * n_layer * d * 1e-6 + cos = np.cos(delta) + sin = 1j * np.sin(delta) + q = n_layer + Mj = np.array([[cos, sin / q], [q * sin, cos]], dtype=complex) + M = M @ Mj + + numerator = ( + q0 * M[0, 0] + + q0 * qs[idx] * M[0, 1] + - M[1, 0] + - qs[idx] * M[1, 1] + ) + denominator = ( + q0 * M[0, 0] + + q0 * qs[idx] * M[0, 1] + + M[1, 0] + + qs[idx] * M[1, 1] + ) + r = numerator / denominator + t = 2 * q0 / denominator + R[idx] = np.abs(r) ** 2 + T[idx] = np.real(qs[idx] / q0) * np.abs(t) ** 2 + + A = np.clip(1 - R - T, 0, 1) + return R, T, A + + +def evaluate_stack(design: Dict) -> Dict: + solar_wl = np.linspace(0.35, 2.5, 120) + ir_wl = np.linspace(8, 13, 200) + solar_w = solar_weight(solar_wl) + ir_w = planck_weight(ir_wl) + + substrate = ag_index + + layer_funcs = [] + thickness = [] + for layer in design["layers"]: + material = layer["material"] + thickness.append(layer["thickness"]) + if material == "PDMS": + layer_funcs.append(pdms_index) + else: + mat = MATERIAL_LIBRARY[material] + layer_funcs.append(lambda wl, m=mat: m.nk(wl)) + + solar_nk = [func(solar_wl) for func in layer_funcs] + ir_nk = [func(ir_wl) for func in layer_funcs] + + solar_R, _, solar_A = transfer_matrix_stack( + solar_wl, solar_nk, thickness, substrate(solar_wl) + ) + ir_R, _, ir_A = transfer_matrix_stack(ir_wl, ir_nk, thickness, substrate(ir_wl)) + + alpha = float(np.trapz(solar_A * solar_w, solar_wl) / np.trapz(solar_w, solar_wl)) + epsilon = float(np.trapz(ir_A * ir_w, ir_wl) / np.trapz(ir_w, ir_wl)) + + score = epsilon - 0.3 * alpha + return {"alpha": alpha, "epsilon": epsilon, "score": score} + + +MATERIAL_LIBRARY: Dict[str, Material] = { + "SiO2": Material("SiO2", 1.45, 1e-4), + "Al2O3": Material("Al2O3", 1.76, 1.5e-3), + "TiO2": Material("TiO2", 2.40, 5e-3), + "Si3N4": Material("Si3N4", 2.05, 2e-3), + "HfO2": Material("HfO2", 1.9, 2e-3), +} + + +def random_design() -> Dict: + num_layers = random.choice([2, 3]) + middle_materials = random.sample(list(MATERIAL_LIBRARY.keys()), num_layers) + layers = [{"material": "PDMS", "thickness": random.uniform(10, 50)}] + for mat in middle_materials: + layers.append( + { + "material": mat, + "thickness": random.uniform(0.05, 2.0), + } + ) + return {"layers": layers} + + +def optimize(iterations: int = 800) -> List[Dict]: + best_designs: List[Dict] = [] + for _ in range(iterations): + design = random_design() + metrics = evaluate_stack(design) + design.update(metrics) + best_designs.append(design) + + best_designs.sort(key=lambda x: x["score"], reverse=True) + return best_designs[:15] + + +def write_summary(designs: List[Dict], path: str) -> None: + with open(path, "w", encoding="utf-8") as f: + f.write("rank,score,epsilon,alpha,layers\n") + for idx, design in enumerate(designs, start=1): + layer_desc = ";".join( + f"{layer['material']}@{layer['thickness']:.3f}um" + for layer in design["layers"] + ) + f.write( + f"{idx},{design['score']:.4f},{design['epsilon']:.4f}," + f"{design['alpha']:.4f},{layer_desc}\n" + ) + + +def plot_pareto(designs: List[Dict], path: str) -> None: + eps = [d["epsilon"] for d in designs] + alpha = [d["alpha"] for d in designs] + scores = [d["score"] for d in designs] + fig, ax = plt.subplots(figsize=(6, 5)) + scatter = ax.scatter(alpha, eps, c=scores, cmap="viridis", s=80) + ax.set_xlabel("Solar-weighted Absorption α") + ax.set_ylabel("8-13 µm Emissivity ε") + ax.set_title("Multilayer Design Performance Distribution") + plt.colorbar(scatter, label="Composite Score ε - 0.3α") + for idx, design in enumerate(designs[:5]): + ax.annotate(str(idx + 1), (design["alpha"], design["epsilon"])) + fig.tight_layout() + plt.savefig(path, dpi=300) + plt.close(fig) + + +def main(): + designs = optimize() + outdir = os.path.join(os.path.dirname(__file__), "outputs") + os.makedirs(outdir, exist_ok=True) + summary_path = os.path.join(outdir, "question3_multilayer_summary.csv") + write_summary(designs, summary_path) + plot_path = os.path.join(outdir, "question3_pareto.png") + plot_pareto(designs, plot_path) + + print(f"Optimal designs written to: {summary_path}") + print(f"Performance scatter plot: {plot_path}") + top = designs[0] + layer_desc = "; ".join( + f"{layer['material']}@{layer['thickness']:.2f}um" for layer in top["layers"] + ) + print( + "Best design: score={:.3f}, ε={:.3f}, α={:.3f}, layers={}".format( + top["score"], top["epsilon"], top["alpha"], layer_desc + ) + ) + + +if __name__ == "__main__": + main() diff --git a/org/other/read_pdf_pypdf.py b/org/other/read_pdf_pypdf.py new file mode 100644 index 0000000..39931a8 --- /dev/null +++ b/org/other/read_pdf_pypdf.py @@ -0,0 +1,41 @@ +import os +from PyPDF2 import PdfReader + + +def find_pdf(start_dir: str, filename: str) -> str: + """Return the absolute path to the target PDF, searching downward if needed.""" + candidate = os.path.join(start_dir, filename) + if os.path.exists(candidate): + return candidate + + for root, _, files in os.walk(start_dir): + if filename in files: + return os.path.join(root, filename) + + raise FileNotFoundError(f"Unable to locate {filename} under {start_dir}") + + +def main() -> None: + script_dir = os.path.dirname(os.path.abspath(__file__)) + pdf_name = "2025 APMCM Problem B.pdf" + pdf_path = find_pdf(script_dir, pdf_name) + + reader = PdfReader(pdf_path) + pages = [] + for idx, page in enumerate(reader.pages, start=1): + text = page.extract_text() or "" + pages.append(f"\n=== Page {idx} ===\n{text.strip()}") + + output_text = "".join(pages) + print(output_text) + + output_path = os.path.join(script_dir, "problem_text_pypdf.txt") + with open(output_path, "w", encoding="utf-8") as f: + f.write(output_text) + + print(f"\nText saved to: {output_path}") + + +if __name__ == "__main__": + main() + diff --git a/org/other/read_pdf_simple.py b/org/other/read_pdf_simple.py new file mode 100644 index 0000000..ed900bf --- /dev/null +++ b/org/other/read_pdf_simple.py @@ -0,0 +1,27 @@ +import pdfplumber +import os + +# 使用当前脚本所在目录 +script_dir = os.path.dirname(os.path.abspath(__file__)) +pdf_path = os.path.join(script_dir, "2025 APMCM Problem B.pdf") + +print(f"Looking for PDF at: {pdf_path}") +print(f"File exists: {os.path.exists(pdf_path)}") + +if os.path.exists(pdf_path): + with pdfplumber.open(pdf_path) as pdf: + full_text = "" + for i, page in enumerate(pdf.pages): + text = page.extract_text() + if text: + full_text += f"\n=== Page {i+1} ===\n{text}\n" + print(full_text) + # 保存到文本文件 + output_path = os.path.join(script_dir, "problem_text.txt") + with open(output_path, "w", encoding="utf-8") as f: + f.write(full_text) + print(f"\n文本已保存到: {output_path}") +else: + print(f"文件不存在: {pdf_path}") + print(f"当前工作目录: {os.getcwd()}") + print(f"目录内容: {os.listdir('.')}") diff --git a/org/use/q1_2.py b/org/use/q1_2.py index 8774206..dc626d5 100644 --- a/org/use/q1_2.py +++ b/org/use/q1_2.py @@ -83,7 +83,7 @@ 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] +thicknesses = [0.5, 1.0, 1.5, 2.0,30,40,50] # ----------------------------- diff --git a/org/use/q2.py b/org/use/q2.py new file mode 100644 index 0000000..8553eb8 --- /dev/null +++ b/org/use/q2.py @@ -0,0 +1,385 @@ +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() \ No newline at end of file diff --git a/org/use/q3.py b/org/use/q3.py new file mode 100644 index 0000000..96bb2fd --- /dev/null +++ b/org/use/q3.py @@ -0,0 +1,349 @@ +import itertools +import math +import os +import random +from dataclasses import dataclass +from typing import Dict, List, Sequence, 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 + + +# ----------------------------- +# 你的PDMS光学常数模型 +# ----------------------------- +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): + """ + 解析分块数据格式： + 第一部分：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 + + +# 全局变量存储PDMS插值函数 +_PDMS_CS_N = None +_PDMS_CS_K = None +_PDMS_WL_MIN = None +_PDMS_WL_MAX = None + + +def initialize_pdms_model(data_file_path): + """ + 初始化PDMS光学常数模型 + """ + global _PDMS_CS_N, _PDMS_CS_K, _PDMS_WL_MIN, _PDMS_WL_MAX + + # 读取数据 + wl_all, n_all, k_all = read_split_data(data_file_path) + wl_all, n_all, k_all = make_strictly_increasing(wl_all, n_all, k_all) + + # 创建插值函数 + _PDMS_CS_N = CubicSpline(wl_all, n_all) + _PDMS_CS_K = CubicSpline(wl_all, k_all) + _PDMS_WL_MIN = wl_all.min() + _PDMS_WL_MAX = wl_all.max() + + print(f"PDMS模型初始化完成: 波长范围 [{_PDMS_WL_MIN:.2f}, {_PDMS_WL_MAX:.2f}] μm") + + +def pdms_nk(wavelength_um: np.ndarray) -> np.ndarray: + """ + 使用你的PDMS模型计算复折射率 + """ + global _PDMS_CS_N, _PDMS_CS_K, _PDMS_WL_MIN, _PDMS_WL_MAX + + if _PDMS_CS_N is None: + raise ValueError("PDMS模型未初始化，请先调用initialize_pdms_model()") + + # 确保波长在有效范围内 + wl_clipped = np.clip(wavelength_um, _PDMS_WL_MIN, _PDMS_WL_MAX) + + # 计算n和k + n = _PDMS_CS_N(wl_clipped) + k = _PDMS_CS_K(wl_clipped) + + # 返回复折射率 + return n - 1j * k + + +# 初始化PDMS模型（请修改为你的实际数据文件路径） +initialize_pdms_model('/Users/spasolreisa/IdeaProjects/asiaMath/data.txt') + + +def planck_weight(wavelength_um: np.ndarray, temperature: float = 300.0) -> np.ndarray: + wl_m = wavelength_um * 1e-6 + c1 = 3.7418e-16 + c2 = 1.4388e-2 + spectral = c1 / (wl_m ** 5 * (np.exp(c2 / (wl_m * temperature)) - 1)) + return spectral + + +def solar_weight(wavelength_um: np.ndarray) -> np.ndarray: + center1, width1 = 0.6, 0.35 + center2, width2 = 1.6, 0.45 + return np.exp(-((wavelength_um - center1) / width1) ** 2) + 0.35 * np.exp( + -((wavelength_um - center2) / width2) ** 2 + ) + + +@dataclass +class Material: + name: str + n_const: float + k_const: float + + def nk(self, wavelength_um: np.ndarray) -> np.ndarray: + n = np.full_like(wavelength_um, self.n_const, dtype=np.complex128) + k = np.full_like(wavelength_um, self.k_const, dtype=np.complex128) + return n - 1j * k + + +def pdms_index(wavelength_um: np.ndarray) -> np.ndarray: + """ + 修改为使用你的PDMS模型 + """ + return pdms_nk(wavelength_um) + + +def ag_index(wavelength_um: np.ndarray) -> np.ndarray: + n = 0.15 + 0.6 * np.exp(-((wavelength_um - 0.5) / 0.4) ** 2) + k = 4.5 + 3.5 * np.exp(-((wavelength_um - 10) / 6) ** 2) + return n - 1j * k + + +def transfer_matrix_stack( + wavelength_um: np.ndarray, + layer_nk: Sequence[np.ndarray], + thickness_um: Sequence[float], + substrate_nk: np.ndarray, + n0: float = 1.0, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + beta = 2 * np.pi / (wavelength_um * 1e-6) + q0 = n0 + qs = substrate_nk + + R = np.zeros_like(wavelength_um) + T = np.zeros_like(wavelength_um) + + for idx, wl in enumerate(wavelength_um): + M = np.identity(2, dtype=complex) + for nk, d in zip(layer_nk, thickness_um): + n_layer = nk[idx] + delta = beta[idx] * n_layer * d * 1e-6 + cos = np.cos(delta) + sin = 1j * np.sin(delta) + q = n_layer + Mj = np.array([[cos, sin / q], [q * sin, cos]], dtype=complex) + M = M @ Mj + + numerator = ( + q0 * M[0, 0] + + q0 * qs[idx] * M[0, 1] + - M[1, 0] + - qs[idx] * M[1, 1] + ) + denominator = ( + q0 * M[0, 0] + + q0 * qs[idx] * M[0, 1] + + M[1, 0] + + qs[idx] * M[1, 1] + ) + r = numerator / denominator + t = 2 * q0 / denominator + R[idx] = np.abs(r) ** 2 + T[idx] = np.real(qs[idx] / q0) * np.abs(t) ** 2 + + A = np.clip(1 - R - T, 0, 1) + return R, T, A + + +def evaluate_stack(design: Dict) -> Dict: + solar_wl = np.linspace(0.35, 2.5, 120) + ir_wl = np.linspace(8, 13, 200) + solar_w = solar_weight(solar_wl) + ir_w = planck_weight(ir_wl) + + substrate = ag_index + + layer_funcs = [] + thickness = [] + for layer in design["layers"]: + material = layer["material"] + thickness.append(layer["thickness"]) + if material == "PDMS": + layer_funcs.append(pdms_index) + else: + mat = MATERIAL_LIBRARY[material] + layer_funcs.append(lambda wl, m=mat: m.nk(wl)) + + solar_nk = [func(solar_wl) for func in layer_funcs] + ir_nk = [func(ir_wl) for func in layer_funcs] + + solar_R, _, solar_A = transfer_matrix_stack( + solar_wl, solar_nk, thickness, substrate(solar_wl) + ) + ir_R, _, ir_A = transfer_matrix_stack(ir_wl, ir_nk, thickness, substrate(ir_wl)) + + alpha = float(np.trapz(solar_A * solar_w, solar_wl) / np.trapz(solar_w, solar_wl)) + epsilon = float(np.trapz(ir_A * ir_w, ir_wl) / np.trapz(ir_w, ir_wl)) + + score = epsilon - 0.3 * alpha + return {"alpha": alpha, "epsilon": epsilon, "score": score} + + +MATERIAL_LIBRARY: Dict[str, Material] = { + "SiO2": Material("SiO2", 1.45, 1e-4), + "Al2O3": Material("Al2O3", 1.76, 1.5e-3), + "TiO2": Material("TiO2", 2.40, 5e-3), + "Si3N4": Material("Si3N4", 2.05, 2e-3), + "HfO2": Material("HfO2", 1.9, 2e-3), +} + + +def random_design() -> Dict: + num_layers = random.choice([2, 3]) + middle_materials = random.sample(list(MATERIAL_LIBRARY.keys()), num_layers) + layers = [{"material": "PDMS", "thickness": random.uniform(10, 50)}] + for mat in middle_materials: + layers.append( + { + "material": mat, + "thickness": random.uniform(0.05, 2.0), + } + ) + return {"layers": layers} + + +def optimize(iterations: int = 800) -> List[Dict]: + best_designs: List[Dict] = [] + for _ in range(iterations): + design = random_design() + metrics = evaluate_stack(design) + design.update(metrics) + best_designs.append(design) + + best_designs.sort(key=lambda x: x["score"], reverse=True) + return best_designs[:15] + + +def write_summary(designs: List[Dict], path: str) -> None: + with open(path, "w", encoding="utf-8") as f: + f.write("rank,score,epsilon,alpha,layers\n") + for idx, design in enumerate(designs, start=1): + layer_desc = ";".join( + f"{layer['material']}@{layer['thickness']:.3f}um" + for layer in design["layers"] + ) + f.write( + f"{idx},{design['score']:.4f},{design['epsilon']:.4f}," + f"{design['alpha']:.4f},{layer_desc}\n" + ) + + +def plot_pareto(designs: List[Dict], path: str) -> None: + eps = [d["epsilon"] for d in designs] + alpha = [d["alpha"] for d in designs] + scores = [d["score"] for d in designs] + fig, ax = plt.subplots(figsize=(6, 5)) + scatter = ax.scatter(alpha, eps, c=scores, cmap="viridis", s=80) + ax.set_xlabel("Solar-weighted Absorption α") + ax.set_ylabel("8-13 µm Emissivity ε") + ax.set_title("Multilayer Design Performance Distribution") + plt.colorbar(scatter, label="Composite Score ε - 0.3α") + for idx, design in enumerate(designs[:5]): + ax.annotate(str(idx + 1), (design["alpha"], design["epsilon"])) + fig.tight_layout() + plt.savefig(path, dpi=300) + plt.close(fig) + + +def main(): + designs = optimize() + outdir = os.path.join(os.path.dirname(__file__), "outputs") + os.makedirs(outdir, exist_ok=True) + summary_path = os.path.join(outdir, "question3_multilayer_summary.csv") + write_summary(designs, summary_path) + plot_path = os.path.join(outdir, "question3_pareto.png") + plot_pareto(designs, plot_path) + + print(f"Optimal designs written to: {summary_path}") + print(f"Performance scatter plot: {plot_path}") + top = designs[0] + layer_desc = "; ".join( + f"{layer['material']}@{layer['thickness']:.2f}um" for layer in top["layers"] + ) + print( + "Best design: score={:.3f}, ε={:.3f}, α={:.3f}, layers={}".format( + top["score"], top["epsilon"], top["alpha"], layer_desc + ) + ) + + +if __name__ == "__main__": + main() \ No newline at end of file