np.piecewise函数用法

文章介绍了numpy.piecewise函数的使用方法,该函数用于根据给定的条件列表对数据进行分段处理。在示例中,第一个例子展示了如何将数值分为三段,分别设置为1、0和2。第二个例子则利用该函数生成了不同频率的波形,通过分段显示100Hz、20Hz和50Hz的正弦波。

np.piecewise函数的语法为:

numpy.piecewise(x, condlist, funclist)

主要参数的含义:

x: 表示要进行操作的对象
condlist: 表示要满足的条件列表,可以是多个条件构成的列表
funclist: 执行的操作列表,参数二与参数三是对应的,当参数二为true的时候,则执行相对应的操作函数

np.piecewise函数的举例应用

举例一:

import numpy as np

x = np.arange(0, 20)
print(x)
# 将小于8的值设为1,将大于15的值设为2;条件之外的值默认为0
print(np.piecewise(x, [x < 8, x >= 15], [1, 2])) 

输出的结果为:

[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19]

[1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 2 2 2 2 2]

举例二:
用np.piecewise()函数来分段显示不同频率的波形

import matplotlib.pyplot as plt
import numpy as np

sampling_rate = 1024
t = np.arange(0, 1.0, 1.0 / sampling_rate)
f1 = 100
f2 = 20
f3 = 50
data = np.piecewise(t, [t < 1, t < 0.8, t < 0.3],
                    [lambda t: np.sin(2 * np.pi * f1 * t),
                     lambda t: np.sin(2 * np.pi * f2 * t),
                     lambda t: np.sin(2 * np.pi * f3 * t)])

plt.figure(figsize=(8, 4))
plt.plot(t, data)
plt.xlabel(u"t/s")
plt.title(u"100Hz'20Hz'50Hz 's shipinpu",  fontsize=20)
plt.show()

结果显示为:
在这里插入图片描述

import numpy as np from scipy.signal import savgol_filter, find_peaks from scipy.optimize import least_squares # 物理常数 PHYS_CONST = { "e": 1.60217662e-19, # 电子电荷 (C) "eps0": 8.854187817e-12, # 真空介电常数 (F/m) "m0": 9.10938356e-31, # 电子静止质量 (kg) "c": 2.99792458e8, # 光速 (m/s) } class SiCThicknessAnalyzer: def __init__(self, n_carrier=5e18, R2=0.20338): self.n_carrier = n_carrier self.m_eff = 0.4 * PHYS_CONST["m0"] self.R2 = R2 def sellmeier(self, lambda_um): """使用Sellmeier方程计算4H-SiC折射率""" A = 1.0 B1 = 0.20075 B2 = 5.54861 B3 = 35.650066 C1 = 12.07224 C2 = 35.65066 C3 = 1268.24708 lambda_sq = lambda_um ** 2 n_sq = A + B1 * lambda_sq / (lambda_sq + C1) + \ B2 * lambda_sq / (lambda_sq - C2) + \ B3 * lambda_sq / (lambda_sq - C3) n_sq = np.clip(n_sq, 2.3 ** 2, 2.8 ** 2) return np.sqrt(n_sq) def drude(self, lambda_um, n_int, n_carrier): """Drude模型修正折射率""" lambda_um = np.asarray(lambda_um) n_int = np.asarray(n_int) n_carrier = np.asarray(n_carrier) n_real = n_int.copy() delta_n = np.zeros_like(n_int) relative_change = np.zeros_like(n_int) omega_p = np.zeros_like(n_int) gamma = np.zeros_like(n_int) mask = n_carrier >= 1e15 if np.any(mask): n_carrier_m = n_carrier[mask] * 1e6 mobility = np.zeros_like(n_carrier_m) for i, nc in enumerate(n_carrier_m): mobility[i] = self._get_mobility(nc / 1e6) omega_p_mask = np.sqrt(n_carrier_m * PHYS_CONST["e"] ** 2 / (self.m_eff * PHYS_CONST["eps0"])) mobility_si = mobility * 1e-4 gamma_mask = PHYS_CONST["e"] / (self.m_eff * mobility_si) lambda_m = lambda_um[mask] * 1e-6 omega = 2 * np.pi * PHYS_CONST["c"] / lambda_m epsilon_complex = n_int[mask] ** 2 - (omega_p_mask ** 2) / (omega ** 2 + 1j * omega * gamma_mask) n_complex = np.sqrt(epsilon_complex) n_real_mask = np.real(n_complex) n_real_mask = np.maximum(n_int[mask] * 0.9, n_real_mask) n_real[mask] = n_real_mask delta_n[mask] = n_real_mask - n_int[mask] relative_change[mask] = (delta_n[mask] / n_int[mask]) * 100 omega_p[mask] = omega_p_mask gamma[mask] = gamma_mask return n_real, delta_n, relative_change, omega_p, gamma def _get_mobility(self, n_carrier): """根据载流子浓度获取迁移率""" if n_carrier < 1e16: return 950 elif n_carrier < 5e16: return 800 elif n_carrier < 1e17: return 650 elif n_carrier < 5e17: return 450 elif n_carrier < 1e18: return 300 elif n_carrier < 5e18: return 150 else: return 80 def preprocess_data(self, wave_number, reflectivity): """数据预处理""" wavelength = 1e4 / wave_number valid_mask = (wavelength >= 3.0) & (wavelength <= 8.0) wavelength = wavelength[valid_mask] reflectivity = reflectivity[valid_mask] if len(wavelength) == 0: raise ValueError("无有效波长数据") n_points = len(reflectivity) window_length = max(11, int(n_points * 0.05)) if window_length % 2 == 0: window_length += 1 reflectivity_smooth = savgol_filter(reflectivity, window_length=window_length, polyorder=3) peaks, _ = find_peaks( reflectivity_smooth, prominence=0.2, distance=10, width=3 ) valleys, _ = find_peaks( -reflectivity_smooth, prominence=0.5, distance=20, width=3 ) return wavelength, reflectivity, reflectivity_smooth, peaks, valleys def calculate_thickness(self, wavelength, reflectivity, extrema_indices, theta_deg, n_carrier=5e18): """使用干涉模型计算厚度""" theta_rad = np.radians(theta_deg) initial_thicknesses = [] refractive_indices = [] phase_correction = True for i in range(len(extrema_indices) - 1): idx1, idx2 = extrema_indices[i], extrema_indices[i + 1] lambda1, lambda2 = wavelength[idx1], wavelength[idx2] if abs(lambda2 - lambda1) < 0.1: continue lambda_avg = (lambda1 + lambda2) / 2 n_int = self.sellmeier(lambda_avg) n, _, _, _, _ = self.drude(lambda_avg, n_int, n_carrier) refractive_indices.append(n) sin_theta_i = np.sin(theta_rad) / n if abs(sin_theta_i) >= 1.0: continue theta_i = np.arcsin(sin_theta_i) if phase_correction: d = (lambda1 * lambda2) / (2 * n * np.cos(theta_i) * abs(lambda2 - lambda1)) else: d = 1 / (2 * n * np.cos(theta_i) * abs(1 / lambda1 - 1 / lambda2)) if 2 < d < 20: initial_thicknesses.append(d) if not initial_thicknesses: raise ValueError("未找到有效的厚度估计") avg_refractive_index = np.mean(refractive_indices) if refractive_indices else 2.65 optimized_params = self._optimize_params( wavelength, reflectivity, np.median(initial_thicknesses), theta_rad, n_carrier ) d_opt, A_opt, B_opt, n_scale_opt = optimized_params fit_quality = self.calculate_fit_quality( wavelength, reflectivity, d_opt, A_opt, B_opt, theta_deg, n_carrier, n_scale_opt ) fit_quality['refractive_index'] = avg_refractive_index * n_scale_opt return d_opt, A_opt, B_opt, fit_quality def _optimize_params(self, wavelength, reflectivity, initial_d, theta_rad, n_carrier): """优化厚度和反射率常数""" def objective(params): d_opt, A, B, phi, n_scale = params residuals = [] for i, wl in enumerate(wavelength): n_int = self.sellmeier(wl) n_real, *_ = self.drude(wl, n_int, n_carrier) n_opt = n_real * n_scale sin_theta_i = np.sin(theta_rad) / n_opt if abs(sin_theta_i) >= 1.0: return np.full(len(wavelength), 1000) theta_i = np.arcsin(sin_theta_i) R1 = ((1 - n_opt) / (1 + n_opt)) ** 2 phase = 4 * np.pi * n_opt * d_opt * np.cos(theta_i) / wl + phi + np.pi numerator = R1 + self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase) denominator = 1 + R1 * self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase) R_theoretical = numerator / denominator R_theoretical *= 100 R_theoretical = np.clip(R_theoretical, 0, 100) residuals.append((reflectivity[i] - R_theoretical) ** 2) return np.array(residuals) A_guess = np.clip(np.mean(reflectivity), 15, 35) B_guess = np.clip((np.max(reflectivity) - np.min(reflectivity)) / 2, 5, 20) initial_guess = [initial_d, A_guess, B_guess, 0, 1.0] bounds = ( [2, 10, 2, -2 * np.pi, 0.95], [20, 40, 25, 2 * np.pi, 1.05] ) result = least_squares( objective, initial_guess, bounds=bounds, max_nfev=3000, ftol=1e-8, xtol=1e-8 ) return result.x[0], result.x[1], result.x[2], result.x[4] def calculate_fit_quality(self, wavelength, reflectivity, d, A, B, theta_deg, n_carrier, n_scale=1.0): """计算拟合质量指标""" theta_rad = np.radians(theta_deg) lambda_avg = np.mean(wavelength) n_int = self.sellmeier(lambda_avg) n_real, *_ = self.drude(lambda_avg, n_int, n_carrier) n_opt = n_real * n_scale sin_theta_i = np.sin(theta_rad) / n_opt if abs(sin_theta_i) >= 1.0: return {"error": "Invalid refraction angle"} theta_i = np.arcsin(sin_theta_i) theoretical = [] for wl in wavelength: R1 = ((1 - n_opt) / (1 + n_opt)) ** 2 phase = 4 * np.pi * n_opt * d * np.cos(theta_i) / wl numerator = R1 + self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase) denominator = 1 + R1 * self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase) R_theoretical = numerator / denominator R_theoretical *= 100 theoretical.append(R_theoretical) theoretical = np.array(theoretical) residuals = reflectivity - theoretical ss_res = np.sum(residuals ** 2) ss_tot = np.sum((reflectivity - np.mean(reflectivity)) ** 2) return { 'r_squared': 1 - (ss_res / ss_tot) if ss_tot > 0 else 0, 'rmse': np.sqrt(np.mean(residuals ** 2)), 'mae': np.mean(np.abs(residuals)), 'refractive_index': n_opt, 'residuals_std': np.std(residuals) } import os import pandas as pd import numpy as np import matplotlib.pyplot as plt from matplotlib import rcParams import core # 设置中文字体 rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans'] rcParams['axes.unicode_minus'] = False def plot_results(self, file_path, wavelength, raw_refl, smooth_refl,peaks, valleys, thickness, A, B, theta_deg, n_carrier, quality): """绘制分析结果(包含质量指标)""" plt.figure(figsize=(15, 12)) # 增加高度以容纳质量指标 # 原始和平滑数据 plt.subplot(3, 1, 1) # 改为3行1列布局 plt.plot(wavelength, raw_refl, 'b-', alpha=0.6, label='原始数据') plt.plot(wavelength, smooth_refl, 'r-', label='平滑数据') # 标记极值点 if len(peaks) > 0: plt.plot(wavelength[peaks], smooth_refl[peaks], 'go', markersize=6, label='峰值') if len(valleys) > 0: plt.plot(wavelength[valleys], smooth_refl[valleys], 'ro', markersize=6, label='谷值') plt.xlabel('波长 (μm)') plt.ylabel('反射率') plt.title(f'{os.path.basename(file_path)} - 预处理结果') plt.legend() plt.grid(True, alpha=0.3) # 理论拟合 plt.subplot(3, 1, 2) plt.plot(wavelength, raw_refl, 'bo', alpha=0.6, markersize=4, label='实验数据') # 计算理论曲线(使用实际波长点计算) theta_rad = np.radians(theta_deg) n_int = self.sellmeier(wavelength) # 改为向量化计算 n_opt, *_ = self.drude(wavelength, n_int, n_carrier) # 改为向量化计算 # 折射角计算(向量化) sin_theta_i = np.sin(theta_rad) / n_opt valid_idx = np.abs(sin_theta_i) < 1.0 theta_i = np.arcsin(np.where(valid_idx, sin_theta_i, 0)) # 理论曲线计算(使用新的反射率公式) R1 = ((1 - n_opt) / (1 + n_opt)) ** 2 phase = 4 * np.pi * n_opt * thickness * np.cos(theta_i) / wavelength numerator = R1 + self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase) denominator = 1 + R1 * self.R2 - 2 * np.sqrt(R1 * self.R2) * np.cos(phase) theory_curve = numerator / denominator theory_curve *= 100 # 转换为百分比 # 只绘制有效的理论曲线部分 plt.plot(wavelength[valid_idx], theory_curve[valid_idx], 'r-', label='理论拟合') plt.xlabel('波长 (μm)') plt.ylabel('反射率') plt.title(f'厚度: {thickness:.3f}μm, 折射率: {np.mean(n_opt):.3f}, R2: {self.R2:.3f}') plt.legend() plt.grid(True, alpha=0.3) # 新增质量指标展示区域 plt.subplot(3, 1, 3) # 隐藏坐标轴 plt.axis('off') # 展示质量指标表格 quality_text = ( f"拟合质量评估:\n\n" f"• 相关系数 R² = {quality['r_squared']:.4f}\n" f"• 均方根误差 RMSE = {quality['rmse']:.4f}%\n" f"• 折射率 = {quality['refractive_index']:.4f}\n" f"• 入射角 = {theta_deg}°\n" f"• 载流子浓度 = {n_carrier:.1e} cm⁻³\n" f"• R2 (常数近似) = {self.R2:.4f}" ) plt.text(0.5, 0.5, quality_text, fontsize=12, ha='center', va='center', bbox=dict(facecolor='white', alpha=0.8)) plt.tight_layout() # 创建规范化的保存路径 save_dir = os.path.dirname(file_path) save_name = os.path.splitext(os.path.basename(file_path))[0] + '_analysis.png' save_path = os.path.join(save_dir, "analysis_results", save_name) os.makedirs(os.path.dirname(save_path), exist_ok=True) plt.savefig(save_path, dpi=300, bbox_inches='tight') plt.close() # 关闭图形以避免内存泄漏 def analyze_file(self, file_path, theta_deg, n_carrier, plot=True): """分析单个光谱文件""" try: # 读取数据 data = pd.read_excel(file_path) wave_number = data.iloc[:, 0].values reflectivity = data.iloc[:, 1].values # 预处理 wavelength, raw_refl, smooth_refl, peaks, valleys = self.preprocess_data( wave_number, reflectivity ) print(f"找到峰值: {len(peaks)}, 谷值: {len(valleys)}") # 合并极值点 extrema = np.sort(np.concatenate((peaks, valleys))) if len(extrema) < 2: raise ValueError("极值点不足,无法计算厚度") # 计算厚度 thickness, A, B, quality = self.calculate_thickness( wavelength, smooth_refl, extrema, theta_deg, n_carrier ) print( f"计算结果: 厚度={thickness:.4f}μm, 折射率={quality['refractive_index']:.4f}, R²={quality['r_squared']:.4f}") # 绘制结果 if plot: self.plot_results(file_path, wavelength, raw_refl, smooth_refl, peaks, valleys, thickness, A, B, theta_deg, n_carrier, quality) return { 'thickness': thickness, 'A': A, 'B': B, 'quality': quality, 'wavelength': wavelength, 'reflectivity': raw_refl } except Exception as e: print(f"分析文件 {file_path} 出错: {str(e)}") import traceback traceback.print_exc() return None def sensitivity_analysis(self, file_path, theta_deg, n_carrier_range, param='carrier'): """ 灵敏度分析: 分析厚度对某一参数的敏感性 :param file_path: 数据文件路径 :param theta_deg: 入射角度 :param n_carrier_range: 参数变化范围 :param param: 分析参数 ('carrier' 或 'angle') :return: 厚度和折射率随参数变化的结果 """ thickness_results = [] refractive_results = [] r_squared_results = [] # 读取和分析数据 data = pd.read_excel(file_path) wave_number = data.iloc[:, 0].values reflectivity = data.iloc[:, 1].values wavelength, raw_refl, smooth_refl, peaks, valleys = self.preprocess_data(wave_number, reflectivity) extrema = np.sort(np.concatenate((peaks, valleys))) if len(extrema) < 2: raise ValueError("极值点不足,无法进行灵敏度分析") for value in n_carrier_range: try: if param == 'carrier': # 分析载流子浓度变化对厚度的影响 thickness, A, B, quality = self.calculate_thickness( wavelength, smooth_refl, extrema, theta_deg, value ) elif param == 'angle': # 分析入射角度变化对厚度的影响 thickness, A, B, quality = self.calculate_thickness( wavelength, smooth_refl, extrema, value, self.n_carrier ) else: raise ValueError(f"未知参数类型: {param}") thickness_results.append(thickness) refractive_results.append(quality['refractive_index']) r_squared_results.append(quality['r_squared']) print(f"{param}={value:.2e} -> 厚度={thickness:.4f}μm, R²={quality['r_squared']:.4f}") except Exception as e: print(f"在{param}={value:.2e}时出错: {str(e)}") thickness_results.append(np.nan) refractive_results.append(np.nan) r_squared_results.append(np.nan) return { 'param_values': n_carrier_range, 'thickness': np.array(thickness_results), 'refractive_index': np.array(refractive_results), 'r_squared': np.array(r_squared_results), 'param_name': param } # 主程序 def main(): """主分析流程""" files_config = [ {"path": "附件1.xlsx", "angle": 10, "carrier": 1e18, "R2": 0.20388}, {"path": "附件2.xlsx", "angle": 15, "carrier": 1e18, "R2": 0.20388} ] results = {} print("=" * 60) print("碳化硅外延层厚度分析程序(基于精确干涉模型)") print("=" * 60) for config in files_config: file_path = config["path"] theta_deg = config["angle"] n_carrier = config["carrier"] R2 = config["R2"] if not os.path.exists(file_path): print(f"文件不存在: {file_path}") continue # 创建分析器实例,传入R2值 analyzer = core.SiCThicknessAnalyzer(n_carrier=n_carrier, R2=R2) print(f"分析文件: {file_path}") print(f" - 入射角: {theta_deg}°") print(f" - 载流子浓度: {n_carrier:.1e} cm⁻³") print(f" - R2 (常数近似): {R2:.3f}") # 分析文件 result = analyzer.analyze_file(file_path, theta_deg, n_carrier) if result: results[file_path] = result # 保留4位小数,符合题目精度要求 print( f"分析完成: 厚度 = {result['thickness']:.4f} μm, 折射率 = {result['quality']['refractive_index']:.4f}") else: print(f"分析失败: {file_path}") # 结果比较和分析 if len(results) >= 2: print("\n" + "=" * 60) print("结果一致性分析") print("=" * 60) # 提取厚度结果 thicknesses = [] refractive_indices = [] angles = [] files = [] for file_path, result in results.items(): thicknesses.append(result['thickness']) refractive_indices.append(result['quality']['refractive_index']) angles.append(next(cfg["angle"] for cfg in files_config if cfg["path"] == file_path)) files.append(os.path.basename(file_path)) # 计算统计量 mean_thickness = np.mean(thicknesses) std_thickness = np.std(thicknesses) max_diff = max(thicknesses) - min(thicknesses) relative_diff = (max_diff / mean_thickness) * 100 if mean_thickness > 0 else 0 mean_refractive = np.mean(refractive_indices) std_refractive = np.std(refractive_indices) # 输出结果 print(" 厚度计算结果:") for i, (file, angle, thickness, refractive) in enumerate(zip(files, angles, thicknesses, refractive_indices)): print(f" {file} ({angle}°): {thickness:.4f} μm, 折射率: {refractive:.4f}") print(f"\n 统计指标:") print(f" 平均厚度: {mean_thickness:.4f} μm") print(f" 厚度标准差: {std_thickness:.4f} μm") print(f" 最大差异: {max_diff:.4f} μm") print(f" 相对差异: {relative_diff:.2f}%") print(f" 平均折射率: {mean_refractive:.4f}") print(f" 折射率标准差: {std_refractive:.4f}") # 可靠性评估 print(f"\n 可靠性评估:") if relative_diff < 3: print(" 优秀 - 结果高度一致") elif relative_diff < 5: print(" 良好 - 结果基本一致") elif relative_diff < 10: print(" 一般 - 存在一定差异") else: print(" 较差 - 差异较大,需要检查") # 拟合质量评估 print(f"\n 拟合质量:") for file_path, result in results.items(): quality = result['quality'] print(f" {os.path.basename(file_path)}: R² = {quality['r_squared']:.4f}, " f"RMSE = {quality['rmse']:.3f}%") elif len(results) == 1: print("\n 只有一个文件分析成功,无法进行一致性比较") file_path, result = list(results.items())[0] quality = result['quality'] print(f"厚度: {result['thickness']:.4f} μm") print(f"折射率: {quality['refractive_index']:.4f}") print(f"拟合质量: R² = {quality['r_squared']:.4f}, RMSE = {quality['rmse']:.3f}%") else: print("\n 没有文件分析成功") return # 灵敏度分析 print("\n" + "=" * 60) print("灵敏度分析") print("=" * 60) sensitivity_results = {} # 分析载流子浓度对厚度的影响 if "附件1.xlsx" in results: print("\n分析附件1的厚度对载流子浓度的敏感性...") n_carrier_range = np.logspace(17, 19, 15) # 从1e17到1e19,取15个点(对数尺度) analyzer = core.SiCThicknessAnalyzer(n_carrier=1e18, R2=0.20388) carrier_sensitivity = analyzer.sensitivity_analysis( "附件1.xlsx", theta_deg=10, n_carrier_range=n_carrier_range, param='carrier' ) sensitivity_results['carrier'] = carrier_sensitivity # 绘制灵敏度曲线 plt.figure(figsize=(12, 8)) # 厚度变化曲线 plt.subplot(2, 1, 1) plt.semilogx(carrier_sensitivity['param_values'], carrier_sensitivity['thickness'], 'bo-') plt.axvline(1e18, color='r', linestyle='--', label=f'标称值 (n=1e18 cm⁻³)') # 修改标签,避免使用上标 plt.xlabel('载流子浓度 n_carrier (cm⁻³)') plt.ylabel('计算出的厚度 d (μm)') plt.title('厚度对载流子浓度的灵敏度分析') plt.grid(True, which="both", ls="--", alpha=0.5) plt.legend() # 折射率变化曲线 plt.subplot(2, 1, 2) plt.semilogx(carrier_sensitivity['param_values'], carrier_sensitivity['refractive_index'], 'go-') plt.axvline(1e18, color='r', linestyle='--', label=f'标称值 (n=1e18 cm⁻³)') # 修改标签,避免使用上标 plt.xlabel('载流子浓度 n_carrier (cm⁻³)') plt.ylabel('折射率 n') plt.title('折射率对载流子浓度的灵敏度分析') plt.grid(True, which="both", ls="--", alpha=0.5) plt.legend() plt.tight_layout() plt.savefig("analysis_results/载流子浓度灵敏度分析.png", dpi=300) plt.show() # 分析入射角度对厚度的影响 if "附件2.xlsx" in results: print("\n分析附件2的厚度对入射角度的敏感性...") angle_range = np.linspace(5, 25, 11) # 从5°到25°,取11个点 analyzer = core.SiCThicknessAnalyzer(n_carrier=1e18, R2=0.20388) angle_sensitivity = analyzer.sensitivity_analysis( "附件2.xlsx", theta_deg=15, n_carrier_range=angle_range, param='angle' ) sensitivity_results['angle'] = angle_sensitivity # 绘制灵敏度曲线 plt.figure(figsize=(12, 8)) # 厚度变化曲线 plt.subplot(2, 1, 1) plt.plot(angle_sensitivity['param_values'], angle_sensitivity['thickness'], 'bo-') plt.axvline(15, color='r', linestyle='--', label=f'标称值 (θ=15°)') plt.xlabel('入射角度 θ (°)') plt.ylabel('计算出的厚度 d (μm)') plt.title('厚度对入射角度的灵敏度分析') plt.grid(True, ls="--", alpha=0.5) plt.legend() # 拟合优度变化曲线 plt.subplot(2, 1, 2) plt.plot(angle_sensitivity['param_values'], angle_sensitivity['r_squared'], 'mo-') plt.axvline(15, color='r', linestyle='--', label=f'标称值 (θ=15°)') plt.xlabel('入射角度 θ (°)') plt.ylabel('拟合优度 R²') plt.title('拟合质量对入射角度的灵敏度分析') plt.grid(True, ls="--", alpha=0.5) plt.legend() plt.tight_layout() plt.savefig("analysis_results/入射角度灵敏度分析.png", dpi=300) plt.show() # 保存灵敏度分析结果 print("\n保存灵敏度分析结果...") for param, result in sensitivity_results.items(): df = pd.DataFrame({ result['param_name']: result['param_values'], '厚度(μm)': result['thickness'], '折射率': result['refractive_index'], 'R²': result['r_squared'] }) df.to_csv(f"analysis_results/{param}_灵敏度分析.csv", index=False) # 保存详细结果 save_detailed_results(results, files_config) print("\n" + "=" * 60) print("分析完成!") print("=" * 60) def save_detailed_results(results, files_config): """保存详细分析结果到文件""" try: # 创建结果目录 os.makedirs("analysis_results", exist_ok=True) # 保存汇总结果 summary_data = [] for config in files_config: file_path = config["path"] if file_path in results: result = results[file_path] # 添加字段存在性检查 summary_data.append({ '文件': os.path.basename(file_path), '入射角(°)': config.get("angle", "N/A"), '载流子浓度(cm⁻³)': config.get("carrier", "N/A"), 'R2': config.get("R2", "N/A"), '厚度(μm)': result.get('thickness', "N/A"), '折射率': result.get('quality', {}).get('refractive_index', "N/A"), 'R²': result.get('quality', {}).get('r_squared', "N/A"), 'RMSE(%)': result.get('quality', {}).get('rmse', "N/A"), }) if summary_data: summary_df = pd.DataFrame(summary_data) summary_file = "analysis_results/厚度分析汇总.csv" summary_df.to_csv(summary_file, index=False, encoding='utf-8-sig') print(f"汇总结果已保存至: {summary_file}") # 将详细数据保存移入循环内部 for file_path, result in results.items(): filename = os.path.splitext(os.path.basename(file_path))[0] # 安全移除扩展名 detail_file = f"analysis_results/{filename}_详细数据.csv" detail_data = { '波长(μm)': result.get('wavelength', []), '反射率(%)': result.get('reflectivity', []) } detail_df = pd.DataFrame(detail_data) detail_df.to_csv(detail_file, index=False, encoding='utf-8-sig') print(f"文件 {filename} 详细数据已保存") # 生成分析报告 generate_report(results, files_config) except Exception as e: print(f"保存过程中发生错误: {str(e)}") # 可添加日志记录或错误上报 def generate_report(results, files_config): """生成分析报告""" report_file = "analysis_results/分析报告.txt" with open(report_file, 'w', encoding='utf-8') as f: f.write("碳化硅外延层厚度分析报告\n") f.write("=" * 50 + "\n\n") f.write("一、分析概述\n") f.write("-" * 30 + "\n") f.write(f"分析时间: {pd.Timestamp.now().strftime('%Y-%m-%d %H:%M:%S')}\n") f.write(f"分析文件数: {len(results)}/{len(files_config)}\n\n") f.write("二、详细结果\n") f.write("-" * 30 + "\n") for config in files_config: file_path = config["path"] if file_path in results: result = results[file_path] f.write(f"文件: {os.path.basename(file_path)}\n") f.write(f" 入射角: {config['angle']}°\n") f.write(f" 载流子浓度: {config['carrier']:.1e} cm⁻³\n") f.write(f" R2 (常数近似): {config['R2']:.4f}\n") f.write(f" 计算厚度: {result['thickness']:.4f} μm\n") f.write(f" 折射率: {result['quality']['refractive_index']:.4f}\n") f.write(f" 拟合优度 R²: {result['quality']['r_squared']:.6f}\n") f.write(f" 均方根误差: {result['quality']['rmse']:.4f}%\n") f.write("\n") # 一致性分析 if len(results) >= 2: f.write("三、一致性分析\n") f.write("-" * 30 + "\n") thicknesses = [result['thickness'] for result in results.values()] refractive_indices = [result['quality']['refractive_index'] for result in results.values()] mean_thickness = np.mean(thicknesses) std_thickness = np.std(thicknesses) max_diff = max(thicknesses) - min(thicknesses) relative_diff = (max_diff / mean_thickness) * 100 if mean_thickness > 0 else 0 mean_refractive = np.mean(refractive_indices) std_refractive = np.std(refractive_indices) f.write(f"平均厚度: {mean_thickness:.4f} μm\n") f.write(f"厚度标准差: {std_thickness:.4f} μm\n") f.write(f"最大差异: {max_diff:.4f} μm\n") f.write(f"相对差异: {relative_diff:.2f}%\n") f.write(f"平均折射率: {mean_refractive:.4f}\n") f.write(f"折射率标准差: {std_refractive:.4f}\n\n") f.write("四、可靠性评估\n") f.write("-" * 30 + "\n") if relative_diff < 3: f.write("优秀 - 不同入射角下的结果高度一致,可靠性很高\n") elif relative_diff < 5: f.write("良好 - 结果基本一致,可靠性较好\n") elif relative_diff < 10: f.write("一般 - 存在一定差异,建议进一步分析\n") else: f.write("较差 - 差异较大,需要检查数据质量和模型参数\n") f.write("\n五、结论\n") f.write("-" * 30 + "\n") if len(results) >= 2: f.write("基于红外干涉法成功测得了碳化硅外延层的厚度和折射率。") f.write("不同入射角下的测量结果表现出良好的一致性,") f.write("表明所建立的数学模型和算法是可靠的。\n") else: f.write("成功完成了碳化硅外延层厚度和折射率的测量。") f.write("建议补充其他入射角的数据以进行交叉验证。\n") print(f" 分析报告已保存至: {report_file}") # 运行示例 if __name__ == "__main__": # 检查文件是否存在 required_files = ["附件1.xlsx", "附件2.xlsx"] missing_files = [f for f in required_files if not os.path.exists(f)] if missing_files: print(" 缺少必要文件:") for f in missing_files: print(f" - {f}") print("\n请确保所有附件文件都在当前目录下") else: main()
最新发布
09-08
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Dream_Bri

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值