循环调用修正sic86

  

create or replace procedure rebuild_sic86_wyl(pi_aac001 in number,
                                              po_fhz    out varchar2,
                                              po_msg    out varchar2) is
  --1.根据账户类型来判断是本地的还是转入的,
  --2.如果是本地的,
  --    写一个游标存放sic86.截止上年末缴费月数jzsnm,本年度缴费月数bn,本年累计缴费月数bnnj,
  --    从最小的年份开始,清空最小年费的下一年jzsnm,bnlj
  --    update sic86 set jzsnm = null ,bnlj = null where aae001 >min(aae001)
  --    update 最小年费下一年的 jzsnm = 最小年份bn+bnlj
  v_count      number(4);
  v_zhlx       sic86.zhlx%type;
  v_aae001_min sic86.aae001%type;
  v_aae001_max sic86.aae001%type;
  v_nf         sic86.aae001%type;
  v_bn         sic86.bn%type;
  v_bnlj       sic86.bnlj%type;
  v_jzsnm      sic86.jzsnm%type;
  cursor c_zhlx_list is
    select zhlx from sic86 where aac001 = pi_aac001;
begin
  po_fhz := '1';
  po_msg := '修正成功';
  for rec_zhlx in c_zhlx_list loop
    --选取最小年份和最小年份
    select min(aae001)
      into v_aae001_min
      from sic86
     where aac001 = pi_aac001;
    select min(aae001)
      into v_aae001_max
      from sic86
     where aac001 = pi_aac001;
    if rec_zhlx.zhlx = '0' then
      for nf in v_aae001_min .. v_aae001_max loop
      --去最小年份的本年缴费月数
        select bn
          into v_bn
          from sic86
         where aac001 = pi_aac001
           and aae001 = nf; 
      --取最小年份的本年累计月数
        select bnlj
          into v_bnlj
          from sic86
         where aac001 = pi_aac001
           and aae001 = nf; 
      --最小年份下一年 v_nf
        v_nf := nf + 1; 
      
        -- 修正1 最小年份下一年的的  jzsnm = 最小年度bn+最小年度bnlj
        update sic86
           set jzsnm =
               (v_bn + v_bnlj)
         where aac001 = pi_aac001
           and aae001 = v_nf;
      --取最小年份下一年的的截至上年末月数
        select jzsnm
          into v_jzsnm 
          from sic86
         where aac001 = pi_aac001
           and aae001 = v_nf;
        -- po_msg :=v_nf||'年份的jzsnm='||v_jzsnm;
        --return;         --调试用,正式用的时候注释掉
        --去最小年份下一年的本年缴费月数   
        select bn
          into v_bn
          from sic86
         where aac001 = pi_aac001
           and aae001 = v_nf; 
      
        --清空最小年份下一年的 截至上年末月数 和 本年累计(其实清空不情况无所谓,不影响)
        /*
          update sic86
            set jzsnm = null, bnlj = null
          where aac001 = pi_aac001
            and aae001 = v_nf;
        */
      
        -- 修正2 最小年份下一年的 本年累计月数 bnlj = 本年的jzsnm+本年的bn
        update sic86
           set bnlj =
               (v_jzsnm + v_bn)
         where aac001 = pi_aac001
           and aae001 = v_nf;
      end loop;
      po_fhz := '1';
      po_msg := '循环修正成功';
    elsif rec_zhlx.zhlx = '1' then
      po_fhz := '-1';
      po_msg := '转入的暂不处理';
      return;
    end if;
  end loop;
end;

  以上只能修正最小年度的下一年。

以下是建表语句:

-- Create table
create table SIC86
(
  jzsnm  NUMBER(4),
  bn     NUMBER(4),
  bnlj   NUMBER(4),
  zhlx   VARCHAR2(2),
  aac001 NUMBER(10),
  aae001 NUMBER(4)
)
tablespace USERS
  pctfree 10
  initrans 1
  maxtrans 255
  storage
  (
    initial 64K
    next 1M
    minextents 1
    maxextents unlimited
  );
-- Add comments to the table 
comment on table SIC86
  is '模拟养老年度账户';
-- Add comments to the columns 
comment on column SIC86.jzsnm
  is '截至上年末缴费月数';
comment on column SIC86.bn
  is '本年缴费月数';
comment on column SIC86.bnlj
  is '截止本年累计月数';
comment on column SIC86.zhlx
  is '账户类型,0:本地; 1:转入生成';
comment on column SIC86.aac001
  is '个人编号';

  最终的修正程序:

可以正常使用,要分两段,第一:最小年份的要单独更新 本年累计bnlj = 截至上年末 jzsnm+本年bn ;

            第二:最大年份的时候同上。

CREATE OR REPLACE PROCEDURE REBUILD_SIC86_WYL(PI_AAC001 IN NUMBER,
                                              PO_FHZ    OUT VARCHAR2,
                                              PO_MSG    OUT VARCHAR2) IS
  --1.根据账户类型来判断是本地的还是转入的,
  --2.如果是本地的,
  --    写一个游标存放sic86.截止上年末缴费月数jzsnm,本年度缴费月数bn,本年累计缴费月数bnnj,
  --    从最小的年份开始,清空最小年费的下一年jzsnm,bnlj
  --    update sic86 set jzsnm = null ,bnlj = null where aae001 >min(aae001)
  --    update 最小年费下一年的 jzsnm = 最小年份bn+bnlj
  V_ZHLX       SIC86.ZHLX%TYPE;
  V_AAE001_MIN SIC86.AAE001%TYPE;
  V_AAE001_MAX SIC86.AAE001%TYPE;
  V_NF         SIC86.AAE001%TYPE;
  V_BN         SIC86.BN%TYPE;
  V_BNLJ       SIC86.BNLJ%TYPE;
  V_JZSNM      SIC86.JZSNM%TYPE;
  CURSOR C_ZHLX_LIST IS
    SELECT ZHLX FROM SIC86 WHERE AAC001 = PI_AAC001;
BEGIN
  PO_FHZ := '1';
  PO_MSG := '修正成功';

  FOR REC_ZHLX IN C_ZHLX_LIST LOOP
    --选取最小年份和最小年份
    SELECT MIN(AAE001)
      INTO V_AAE001_MIN
      FROM SIC86
     WHERE AAC001 = PI_AAC001;
    SELECT MAX(AAE001)
      INTO V_AAE001_MAX
      FROM SIC86
     WHERE AAC001 = PI_AAC001;
  
    FOR NF IN V_AAE001_MIN .. V_AAE001_MAX LOOP
        --修正01  最小年份的本年累计,只修正最小年份,只修正一次
        IF (NF = V_AAE001_MIN) THEN
          --取最小年份的截止上年末月数,这个只用一次 
          UPDATE SIC86
             SET BNLJ =
                 (NVL(JZSNM, 0) + NVL(BN, 0))
           WHERE AAC001 = PI_AAC001
             AND AAE001 = NF;
        ELSE
          NULL;
        END IF;
      
        --去最小年份的本年缴费月数
        SELECT BN
          INTO V_BN
          FROM SIC86
         WHERE AAC001 = PI_AAC001
           AND AAE001 = NF;
        --取最小年份的本年累计月数
        SELECT BNLJ
          INTO V_BNLJ
          FROM SIC86
         WHERE AAC001 = PI_AAC001
           AND AAE001 = NF;
      
        --最小年份下一年 v_nf
        V_NF := NF + 1;
      
        -- 修正1 最小年份下一年的的  jzsnm = 最小年度bnlj
        UPDATE SIC86
           SET JZSNM = V_BNLJ
         WHERE AAC001 = PI_AAC001
           AND AAE001 = V_NF;
        --加入判断截至上年末下一年v_nf是否为最后一年,如果直接退出 
        IF V_NF = V_AAE001_MAX THEN
          --取最小年份的截止上年末月数,这个只用一次 
          UPDATE SIC86
             SET BNLJ =
                 (NVL(JZSNM, 0) + NVL(BN, 0))
           WHERE AAC001 = PI_AAC001
             AND AAE001 = V_NF;
          PO_FHZ := '1';
          PO_MSG := '修复成功';
          RETURN;
        END IF;
        --取最小年份下一年的的截至上年末月数
        SELECT JZSNM
          INTO V_JZSNM
          FROM SIC86
         WHERE AAC001 = PI_AAC001
           AND AAE001 = V_NF;
        -- po_msg :=v_nf||'年份的jzsnm='||v_jzsnm;
        --return;         --调试用,正式用的时候注释掉
        --去最小年份下一年的本年缴费月数
        SELECT BN
          INTO V_BN
          FROM SIC86
         WHERE AAC001 = PI_AAC001
           AND AAE001 = V_NF;
      
        -- 修正2 最小年份下一年的 本年累计月数 bnlj = 本年的jzsnm+本年的bn
        UPDATE SIC86
           SET BNLJ =
               (V_JZSNM + V_BN)
         WHERE AAC001 = PI_AAC001
           AND AAE001 = V_NF;
      END LOOP;
      PO_FHZ := '1';
      PO_MSG := '循环修正成功';
  END LOOP;
END;

  

我正在编辑【python】代码,遇到了 【Unresolved attribute reference 'analyze_file' for class 'SiCThicknessAnalyzer'】 ,请帮我检查并改正错误点。我的原始代码如下: 【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 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()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) - 真空中的光速 } # 碳化硅外延层厚度计算模型 # 基于红外干涉法和物理光学模型计算4H-SiC外延层的厚度 class SiCThicknessAnalyzer: # 初始化函数,设定默认参数 # n_carrier: 载流子浓度,单位cm⁻³,默认5×10¹⁸ # R2: SiC-SiC界面反射率,默认0.20338(根据文献资料) def __init__(self, n_carrier=5e18, R2=0.20338): self.n_carrier = n_carrier # 载流子浓度 (cm^-3)影响Drude模型修正 self.m_eff = 0.4 * PHYS_CONST["m0"] # 4H-SiC电子有效质量 (kg) self.R2 = R2 # 常数近似法中的R2值(SiC-SiC界面反射率)- 衬底反射率 # Sellmeier方程计算4H-SiC折射率 # lambda_um: 波长,单位微米(μm) # 返回: 对应波长的折射率 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 # Sellmeier方程计算折射率平方 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) # 物理限制:4H-SiC折射率范围(2.3~2.8) # 确保计算结果在物理合理范围内,避免数值误差导致的异常值 n_sq = np.clip(n_sq, 2.3 ** 2, 2.8 ** 2) return np.sqrt(n_sq) # 返回折射率 # Drude模型修正折射率 # 考虑载流子浓度对折射率的影响(自由载流子吸收效应) # lambda_um: 波长数组,单位微米 # n_int: 本征折射率数组(Sellmeier方程计算结果) # n_carrier: 载流子浓度数组 # 返回: 修正后的折射率及相关参数 def drude(self, lambda_um, n_int, n_carrier): # 将输入转换为数组,确保处理标量和数组的一致性 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) # 散射率 # 找出需要Drude修正的点 # 只有载流子浓度足够高时才需要修正(n_carrier >= 1e15) mask = n_carrier >= 1e15 if np.any(mask): # 单位转换 cm⁻³ → m⁻³ # 1 cm⁻³ = 10⁶ m⁻³(因为1 m³ = 10⁶ cm³) 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) # 转换回cm⁻³单位输入 # 等离子体频率计算 omega_p_mask = np.sqrt(n_carrier_m * PHYS_CONST["e"] ** 2 / (self.m_eff * PHYS_CONST["eps0"])) # 散射率计算 mobility_si = mobility * 1e-4 # 单位转换,cm²/V·s → m²/V·s gamma_mask = PHYS_CONST["e"] / (self.m_eff * mobility_si) # 光波角频率计算 lambda_m = lambda_um[mask] * 1e-6 # μm → m omega = 2 * np.pi * PHYS_CONST["c"] / lambda_m # Drude模型修正介电函数 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) # 取实部(折射率) # 物理合理性检查:如果修正后的折射率大于本征折射率,则进行限制 # 选择将修正后的折射率限制在本征折射率的90%以上 # Drude修正通常使折射率减小,但数值计算可能产生异常 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 # 根据载流子浓度获取迁移率 # 迁移率随载流子浓度增加而减小(散射增强) # n_carrier: 载流子浓度,单位cm⁻³ # 返回: 迁移率,单位cm²/V·s 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 # 高浓度时迁移率低 # 数据预处理(筛选有效干涉条纹) # wave_number: 波数数据,单位cm⁻¹ # reflectivity: 反射率数据,单位% # 返回: 处理后的波长、反射率数据及极值点位置 def preprocess_data(self, wave_number, reflectivity): # 波数转波长 wavelength = 1e4 / wave_number # 限制波长为3~8μm(碳化硅透明窗口) # 此范围外数据可能受吸收影响,不适合干涉分析 valid_mask = (wavelength >= 3.0) & (wavelength <= 8.0) wavelength = wavelength[valid_mask] reflectivity = reflectivity[valid_mask] if len(wavelength) == 0: raise ValueError("无有效波长数据(需3~8μm范围)") # 平滑处理 - 使用Savitzky-Golay滤波器去除噪声 # 保持信号特征的同时减少随机噪声 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, # 最小 prominence 确保是显著峰值 distance=10, # 最小 peak 间距 width=3 # 最小 peak 宽度 ) valleys, _ = find_peaks( -reflectivity_smooth, # 找反射率的谷值(即-reflectivity的峰值) prominence=0.5, # 谷值需要更大的 prominence distance=20, # 谷值间距通常较大 width=3 ) return wavelength, reflectivity, reflectivity_smooth, peaks, valleys # 使用干涉模型计算厚度 # wavelength: 波长数组 # reflectivity: 反射率数组 # extrema_indices: 极值点索引数组 # theta_deg: 入射角度,单位度 # n_carrier: 载流子浓度 # 返回: 优化后的厚度和拟合质量信息 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 # 计算平均波长处的折射率(Sellmeier+Drude) 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)) # 碳化硅外延层工艺厚度范围(2~20μm) # 过滤不合理的结果 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, n_scale_opt = optimized_params # 计算拟合质量 fit_quality = self.calculate_fit_quality( wavelength, reflectivity, d_opt, theta_deg, n_carrier, n_scale_opt ) # 添加折射率到质量评估中 fit_quality['refractive_index'] = avg_refractive_index * n_scale_opt return d_opt, fit_quality # 优化厚度和反射率常数 # 使用最小二乘法优化参数,使理论反射率最接近实验值 def _optimize_params(self, wavelength, reflectivity, initial_d, theta_rad, n_carrier): def objective(params): # 参数: [厚度d, 相位φ, 折射率缩放因子n_scale] d_opt, 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) # Drude修正 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(空气-外延层界面反射率) # 菲涅尔反射公式 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 # 反射率物理限制:0~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) # 反射率基线15~35%(符合空气-SiC反射特性) B_guess = np.clip((np.max(reflectivity) - np.min(reflectivity)) / 2, 5, 20) # 振幅5~20% initial_guess = [initial_d, A_guess, B_guess, 0, 1.0] # [厚度, A, B, 相位, 折射率缩放] # 参数边界:严格贴合碳化硅工程范围 bounds = ( [2, 10, 2, -2 * np.pi, 0.95], # 下限:厚度≥2μm,A≥10%,B≥2%,n_scale≥0.95 [20, 40, 25, 2 * np.pi, 1.05] # 上限:厚度≤20μm,A≤40%,B≤25%,n_scale≤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(空气-外延层界面反射率) 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, # 决定系数,越接近1越好 'rmse': np.sqrt(np.mean(residuals ** 2)), # 均方根误差,越小越好 'mae': np.mean(np.abs(residuals)), # 平均绝对误差,越小越好 'refractive_index': n_opt, # 优化后的折射率 'residuals_std': np.std(residuals) # 残差标准差 } 】
09-08
import os import pandas as pd import numpy as np from scipy.signal import savgol_filter, find_peaks from scipy.optimize import least_squares import matplotlib.pyplot as plt from matplotlib import rcParams # 设置中文字体 rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans'] rcParams['axes.unicode_minus'] = False # 物理常数 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): self.n_carrier = n_carrier # 载流子浓度 (cm^-3) self.m_eff = 0.4 * PHYS_CONST["m0"] # 4H-SiC电子有效质量 (kg) def sellmeier(self, lambda_um): """使用Sellmeier方程计算4H-SiC折射率""" # Sellmeier系数 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) # 物理限制:4H-SiC折射率范围 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模型修正折射率""" if n_carrier < 1e15: return n_int, 0, 0, 0, 0 try: # 单位转换 n_carrier_m = n_carrier * 1e6 # cm⁻³ → m⁻³ # 迁移率与载流子浓度关系 mobility = self._get_mobility(n_carrier) # 等离子体频率 omega_p = np.sqrt(n_carrier_m * PHYS_CONST["e"] ** 2 / (self.m_eff * PHYS_CONST["eps0"])) # 散射率 mobility_si = mobility * 1e-4 # cm²/V·s → m²/V·s gamma = PHYS_CONST["e"] / (self.m_eff * mobility_si) # 光波角频率 lambda_m = lambda_um * 1e-6 omega = 2 * np.pi * PHYS_CONST["c"] / lambda_m # Drude模型修正 epsilon_complex = n_int ** 2 - (omega_p ** 2) / (omega ** 2 + 1j * omega * gamma) n_complex = np.sqrt(epsilon_complex) n_real = np.real(n_complex) # 物理合理性检查 if n_real > n_int: n_real = max(n_int * 0.9, n_real) delta_n = n_real - n_int relative_change = (delta_n / n_int) * 100 if n_int > 1e-10 else 0 return n_real, delta_n, relative_change, omega_p, gamma except Exception as e: print(f"Drude修正错误: {e}") return n_int, 0, 0, 0, 0 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): """数据预处理(增加B题专用的干涉条纹周期性验证,筛选有效条纹)""" # 波数转波长(B题附件1/2第1列波数→波长,单位适配) wavelength = 1e4 / wave_number # B题专用:限制波长为3~8μm(碳化硅透明窗口,排除非干涉区域) valid_mask = (wavelength >= 3.0) & (wavelength <= 8.0) wavelength = wavelength[valid_mask] reflectivity = reflectivity[valid_mask] if len(wavelength) == 0: raise ValueError("B题:无有效波长数据(需3~8μm范围)") # 平滑处理:适配B题实测数据噪声,窗口长度调整为数据量的5%(动态适配) n_points = len(reflectivity) window_length = max(11, int(n_points * 0.05)) # 原固定11,改为动态窗口 if window_length % 2 == 0: window_length += 1 reflectivity_smooth = savgol_filter(reflectivity, window_length=window_length, polyorder=3) # 极值点提取:B题专用参数(提升条纹识别准确性) peaks, _ = find_peaks( reflectivity_smooth, prominence=0.5, # 原0.2,提高阈值过滤小噪声 distance=20, # 原10,增加距离筛选非连续条纹 width=3 # 新增宽度约束,排除尖锐噪声 ) valleys, _ = find_peaks( -reflectivity_smooth, prominence=0.5, distance=20, width=3 ) # B题核心:验证干涉条纹周期性(双光束干涉的关键特性,问题1隐含要求) def validate_periodicity(extrema_lambdas): if len(extrema_lambdas) < 3: return True # 至少3个点才验证周期性 # 双光束干涉周期性:1/λ与干涉级次k线性相关(R²≥0.95为有效) k = np.arange(len(extrema_lambdas)) wave_number = 1 / extrema_lambdas from scipy.stats import linregress lin_res = linregress(k, wave_number) return lin_res.rvalue ** 2 >= 0.95 # 周期性R²阈值,符合B题可靠性要求 # 合并并验证极值点周期性 all_extrema = np.sort(np.concatenate((peaks, valleys))) if not validate_periodicity(wavelength[all_extrema]): print("B题:警告!部分极值点不满足周期性,已过滤非干涉条纹") # 保留前N个满足周期性的极值点 for n in range(len(all_extrema), 2, -1): if validate_periodicity(wavelength[all_extrema[:n]]): all_extrema = all_extrema[:n] break # 拆分峰值和谷值(基于筛选后的极值点) peaks = [p for p in peaks if p in all_extrema] valleys = [v for v in valleys if v in all_extrema] return wavelength, reflectivity, reflectivity_smooth, peaks, valleys def calculate_thickness(self, wavelength, reflectivity, extrema_indices, theta_deg, n_carrier=5e18): """使用干涉模型计算厚度(显式修正相位差,贴合B题问题1双光束模型)""" theta_rad = np.radians(theta_deg) initial_thicknesses = [] refractive_indices = [] # 存储每个极值点对计算的折射率 # B题问题1:双光束干涉的相位差修正(空气→外延层界面反射有π相位突变,等效λ/2光程差) phase_correction = True # 启用相位差修正,完全匹配B题模型 for i in range(len(extrema_indices) - 1): idx1, idx2 = extrema_indices[i], extrema_indices[i + 1] lambda1, lambda2 = wavelength[idx1], wavelength[idx2] # B题实测数据特性:过滤过近的极值点(避免噪声虚假条纹) if abs(lambda2 - lambda1) < 0.1: # 从0.05μm收紧到0.1μm,贴合B题有效干涉条纹间隔 continue # 计算平均波长处的折射率(Sellmeier+Drude,符合B题"折射率与波长、载流子相关"设定) lambda_avg = (lambda1 + lambda2) / 2 n_int = self.sellmeier(lambda_avg) n, _, _, _, _ = self.drude(lambda_avg, n_int, n_carrier) refractive_indices.append(n) # 记录折射率 # 斯涅尔定律计算折射角(适配B题非垂直入射场景) sin_theta_i = np.sin(theta_rad) / n if abs(sin_theta_i) >= 1.0: continue theta_i = np.arcsin(sin_theta_i) # B题问题1核心公式:显式包含相位差修正 if phase_correction: # 修正后公式:2ndcosθ_i + λ/2 = kλ → d = λ1λ2/(2ncosθ_i|λ2-λ1|)(与原公式形式一致,但显式标注相位差) 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)) d -= 10 # B题物理约束:碳化硅外延层工艺厚度范围(2~20μm,问题2隐含要求) 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): """优化厚度和反射率常数(贴合B题物理约束,避免非物理解)""" def objective(params): d_opt, A, B, phi, n_scale = params residuals = [] # B题优化:每个波长单独计算折射率(避免平均波长导致的色散误差) 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) # 显式包含相位差(π),匹配B题问题1模型 phase = 4 * np.pi * n_opt * d_opt * np.cos(theta_i) / wl + phi + np.pi R_theoretical = A + B * np.cos(phase) # B题物理约束:反射率0~100% R_theoretical = np.clip(R_theoretical, 0, 100) residuals.append((reflectivity[i] - R_theoretical) ** 2) return np.array(residuals) # B题专用初始值:基于碳化硅特性设定 A_guess = np.clip(np.mean(reflectivity), 15, 35) # 反射率基线15~35%(符合空气-SiC反射特性) B_guess = np.clip((np.max(reflectivity) - np.min(reflectivity)) / 2, 5, 20) # 振幅5~20% initial_guess = [initial_d, A_guess, B_guess, 0, 1.0] # 相位初始值0,折射率缩放1.0 # B题专用参数边界:严格贴合碳化硅工程范围 bounds = ( [2, 10, 2, -2 * np.pi, 0.95], # 下限:厚度≥2μm,A≥10%,B≥2%,n_scale≥0.95 [20, 40, 25, 2 * np.pi, 1.05] # 上限:厚度≤20μm,A≤40%,B≤25%,n_scale≤1.05 ) # 优化算法:增加迭代次数,确保收敛(适配B题复杂实测数据) result = least_squares( objective, initial_guess, bounds=bounds, max_nfev=3000, # 原1000,增加迭代次数 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: phase = 4 * np.pi * n_opt * d * np.cos(theta_i) / wl R_theoretical = A + B * np.cos(phase) 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) } def analyze_file(self, file_path, theta_deg, n_carrier, plot=True): """分析单个光谱文件,增加plot参数控制绘图""" 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(): """主分析流程(完全适配B题附件1/2处理)""" # B题附件配置:严格对应题目"附件1=10°,附件2=15°"设定 files_config = [ {"path": "附件1.xlsx", "angle": 10, "carrier": 1e18}, # 可根据B题实际载流子浓度调整 {"path": "附件2.xlsx", "angle": 15, "carrier": 1e18} ] results = {} # 创建分析器实例 - 使用默认载流子浓度 analyzer = SiCThicknessAnalyzer(n_carrier=1e18) print("=" * 60) print("B题 碳化硅外延层厚度分析程序(基于双光束干涉模型)") print("=" * 60) for config in files_config: file_path = config["path"] theta_deg = config["angle"] n_carrier = config["carrier"] if not os.path.exists(file_path): print(f"❌ B题:文件不存在: {file_path}(需与附件1/2文件名一致)") continue print(f"\n🔍 B题:正在分析 {file_path}") print(f" - 入射角: {theta_deg}°(匹配题目附件设定)") print(f" - 载流子浓度: {n_carrier:.1e} cm⁻³(贴合B题外延层特性)") # 分析文件:调用修改后的B题专用逻辑 result = analyzer.analyze_file(file_path, theta_deg, n_carrier) if result: results[file_path] = result # B题结果输出:保留4位小数,符合题目精度要求 print( f"✅ B题:分析完成: 厚度 = {result['thickness']:.4f} μm, 折射率 = {result['quality']['refractive_index']:.4f}") else: print(f"❌ B题:分析失败: {file_path}") # 结果比较和分析 if len(results) >= 2: # 灵敏度分析(以附件1为例) 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个点(对数尺度) 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:.1e} 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:.1e} 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个点 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) 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["angle"], '载流子浓度(cm⁻³)': config["carrier"], '厚度(μm)': result['thickness'], '折射率': result['quality']['refractive_index'], 'R²': result['quality']['r_squared'], 'RMSE(%)': result['quality']['rmse'], }) 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.basename(file_path).replace('.xlsx', '') detail_file = f"analysis_results/{filename}_详细数据.csv" detail_data = { '波长(μm)': result['wavelength'], '反射率(%)': result['reflectivity'] } detail_df = pd.DataFrame(detail_data) detail_df.to_csv(detail_file, index=False, encoding='utf-8-sig') # 生成分析报告 generate_report(results, files_config) except Exception as e: print(f"保存结果时出错: {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" 计算厚度: {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
import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.signal import find_peaks, savgol_filter, welch from scipy.fft import fft, fftfreq, fftshift from scipy.optimize import minimize_scalar, curve_fit from scipy.interpolate import interp1d import warnings warnings.filterwarnings('ignore') plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False class EnhancedInterferenceAnalyzer: """增强的干涉测厚分析器 - 支持动态折射率和优化的傅里叶分析""" def __init__(self): print("增强干涉测厚分析器") print("=" * 60) print("新增功能:") print("1. 动态折射率计算:n = f(λ, 载流子浓度)") print("2. 优化傅里叶分析:窗函数、零填充、功率谱密度") print("3. 多波段分析:分波段计算提高精度") print("4. 不确定性评估:计算测量不确定度") print("=" * 60) # SiC基本参数 self.n0_sic = 2.6 # 基础折射率 self.wavelength_ref = 10.0 # 参考波长 (μm) # 色散参数 (Cauchy方程: n = A + B/λ² + C/λ⁴) self.cauchy_A = 2.6 self.cauchy_B = 0.05 # μm² self.cauchy_C = 0.001 # μm⁴ def calculate_dynamic_refractive_index(self, wavelength_um, carrier_concentration=1e16): """ 计算动态折射率 参数: wavelength_um: 波长 (μm) carrier_concentration: 载流子浓度 (cm^-3) 返回: n: 折射率 """ # 基础色散 (Cauchy方程) n_base = (self.cauchy_A + self.cauchy_B / (wavelength_um ** 2) + self.cauchy_C / (wavelength_um ** 4)) # 载流子对折射率的影响 (Drude模型) plasma_frequency = np.sqrt(carrier_concentration * 1.6e-19 ** 2 / (8.85e-12 * 0.5 * 9.11e-31)) wavelength_m = wavelength_um * 1e-6 frequency = 3e8 / wavelength_m # 简化的载流子修正 carrier_correction = -plasma_frequency ** 2 / (2 * frequency ** 2) * 0.1 n_corrected = n_base + carrier_correction return np.clip(n_corrected, 1.5, 4.0) # 限制在合理范围内 def enhanced_fourier_analysis(self, wavenumbers, reflectances): """ 增强的傅里叶分析 改进: 1. 使用汉宁窗减少频谱泄露 2. 零填充提高频率分辨率 3. 功率谱密度分析 4. 多峰检测 """ print("\n增强傅里叶分析:") # 确保数据等间距 k_min, k_max = wavenumbers.min(), wavenumbers.max() n_points = len(wavenumbers) k_uniform = np.linspace(k_min, k_max, n_points) # 插值到等间距网格 interp_func = interp1d(wavenumbers, reflectances, kind='cubic', bounds_error=False, fill_value='extrapolate') r_uniform = interp_func(k_uniform) # 预处理:去趋势和去均值 r_detrended = self.detrend_data(k_uniform, r_uniform) r_centered = r_detrended - np.mean(r_detrended) # 应用汉宁窗 window = np.hanning(len(r_centered)) r_windowed = r_centered * window # 零填充 (提高频率分辨率) padding_factor = 4 n_padded = n_points * padding_factor r_padded = np.zeros(n_padded) r_padded[:n_points] = r_windowed # FFT fft_result = fft(r_padded) freqs = fftfreq(n_padded, d=(k_max - k_min) / n_points) # 功率谱密度 power_spectrum = np.abs(fft_result) ** 2 # 归一化 power_spectrum = power_spectrum / np.max(power_spectrum) # 只考虑正频率 positive_mask = freqs > 0 positive_freqs = freqs[positive_mask] positive_power = power_spectrum[positive_mask] # 平滑功率谱以减少噪声 if len(positive_power) > 20: positive_power_smooth = savgol_filter(positive_power, 11, 3) else: positive_power_smooth = positive_power # 多峰检测 peak_indices, peak_props = find_peaks( positive_power_smooth, prominence=0.1, # 相对于最大值的10% height=0.05, # 至少5%的高度 distance=len(positive_power_smooth) // 20 # 最小间距 ) dominant_frequencies = [] dominant_powers = [] if len(peak_indices) > 0: # 按功率排序 sorted_indices = peak_indices[np.argsort(positive_power_smooth[peak_indices])[::-1]] for i, peak_idx in enumerate(sorted_indices[:3]): # 取前3个主要峰 freq = positive_freqs[peak_idx] power = positive_power_smooth[peak_idx] period = 1 / freq if freq > 0 else 0 dominant_frequencies.append(freq) dominant_powers.append(power) print(f"主要频率 {i + 1}: {freq:.6f} cm, 功率: {power:.3f}") print(f"对应周期: {period:.2f} cm^-1") return dominant_frequencies, dominant_powers, positive_freqs, positive_power_smooth def detrend_data(self, x, y): """去除数据趋势""" # 使用二次多项式拟合趋势 coeffs = np.polyfit(x, y, 2) trend = np.polyval(coeffs, x) return y - trend def multi_band_analysis(self, wavenumbers, reflectances, theta_deg): """ 多波段分析 将光谱分为多个波段分别分析,提高精度 """ print(f"\n多波段分析:") # 转换为波长 wavelengths = 10000 / wavenumbers # cm^-1 to μm # 定义波段 bands = [ (2.5, 5.0, "短波红外"), (5.0, 10.0, "中波红外"), (10.0, 25.0, "长波红外") ] band_results = [] for wl_min, wl_max, band_name in bands: # 选择波段数据 band_mask = (wavelengths >= wl_min) & (wavelengths <= wl_max) if np.sum(band_mask) < 20: # 数据点太少 continue band_wavenumbers = wavenumbers[band_mask] band_reflectances = reflectances[band_mask] band_wavelengths = wavelengths[band_mask] print(f"\n{band_name} ({wl_min}-{wl_max} μm):") print(f"数据点数: {len(band_wavenumbers)}") # 计算该波段的平均折射率 avg_wavelength = np.mean(band_wavelengths) avg_refractive_index = self.calculate_dynamic_refractive_index(avg_wavelength) print(f"平均波长: {avg_wavelength:.2f} μm") print(f"平均折射率: {avg_refractive_index:.3f}") # 傅里叶分析 dominant_freqs, powers, freqs, power_spectrum = self.enhanced_fourier_analysis( band_wavenumbers, band_reflectances ) # 计算厚度 if dominant_freqs: for i, freq in enumerate(dominant_freqs[:2]): # 取前两个主要频率 period = 1 / freq if freq > 0 else 0 if period > 0: thickness = self.calculate_thickness_from_period_dynamic( period, theta_deg, avg_refractive_index ) band_results.append({ 'band': band_name, 'wavelength_range': (wl_min, wl_max), 'avg_wavelength': avg_wavelength, 'refractive_index': avg_refractive_index, 'period': period, 'thickness': thickness, 'power': powers[i] if i < len(powers) else 0 }) print(f"厚度 (频率{i + 1}): {thickness:.3f} μm") return band_results def calculate_thickness_from_period_dynamic(self, period_cm_inv, theta_deg, refractive_index): """使用动态折射率计算厚度""" theta_rad = np.radians(theta_deg) # 计算层内传播角度 sin_theta1 = np.sin(theta_rad) / refractive_index if sin_theta1 >= 1: cos_theta1 = 0.1 # 避免全反射情况 else: cos_theta1 = np.sqrt(1 - sin_theta1 ** 2) # 厚度计算 thickness_cm = 1 / (2 * refractive_index * cos_theta1 * period_cm_inv) thickness_um = thickness_cm * 1e4 # cm → μm return thickness_um def uncertainty_analysis(self, results_list): """不确定性分析""" print(f"\n不确定性分析:") if not results_list: return None thicknesses = [r['thickness'] for r in results_list] powers = [r['power'] for r in results_list] # 加权平均 (用功率作为权重) total_power = sum(powers) if total_power > 0: weights = [p / total_power for p in powers] weighted_mean = sum(t * w for t, w in zip(thicknesses, weights)) else: weighted_mean = np.mean(thicknesses) # 标准差 std_dev = np.std(thicknesses) # 相对不确定度 relative_uncertainty = std_dev / weighted_mean * 100 if weighted_mean > 0 else 0 print(f"测量结果数: {len(thicknesses)}") print(f"厚度范围: {min(thicknesses):.3f} - {max(thicknesses):.3f} μm") print(f"加权平均: {weighted_mean:.3f} μm") print(f"标准差: {std_dev:.3f} μm") print(f"相对不确定度: {relative_uncertainty:.2f}%") return { 'weighted_mean': weighted_mean, 'std_dev': std_dev, 'relative_uncertainty': relative_uncertainty, 'count': len(thicknesses) } def load_data(self, filename): """加载数据""" try: df = pd.read_excel(filename) wavenumbers = df.iloc[:, 0].values # cm^-1 reflectances = df.iloc[:, 1].values # % # 基本清理 valid_mask = np.isfinite(wavenumbers) & np.isfinite(reflectances) wavenumbers = wavenumbers[valid_mask] reflectances = reflectances[valid_mask] # 排序 sort_idx = np.argsort(wavenumbers) wavenumbers = wavenumbers[sort_idx] reflectances = reflectances[sort_idx] print(f"\n数据加载: {filename}") print(f"数据点数: {len(wavenumbers)}") print(f"波数范围: {wavenumbers.min():.1f} - {wavenumbers.max():.1f} cm^-1") print(f"波长范围: {10000 / wavenumbers.max():.2f} - {10000 / wavenumbers.min():.2f} μm") return wavenumbers, reflectances except Exception as e: print(f"数据加载失败: {e}") return None, None def comprehensive_analysis(self, filename, theta_deg): """综合分析""" print(f"\n" + "=" * 80) print(f"增强综合分析: {filename} (入射角: {theta_deg}°)") print("=" * 80) # 加载数据 wavenumbers, reflectances = self.load_data(filename) if wavenumbers is None: return None # 多波段分析 band_results = self.multi_band_analysis(wavenumbers, reflectances, theta_deg) # 不确定性分析 uncertainty_result = self.uncertainty_analysis(band_results) # 整体傅里叶分析 print(f"\n整体傅里叶分析:") dominant_freqs, powers, freqs, power_spectrum = self.enhanced_fourier_analysis( wavenumbers, reflectances ) # 使用平均折射率计算整体厚度 avg_wavelength = np.mean(10000 / wavenumbers) avg_n = self.calculate_dynamic_refractive_index(avg_wavelength) overall_thicknesses = [] if dominant_freqs: for freq in dominant_freqs[:2]: period = 1 / freq if freq > 0 else 0 if period > 0: thickness = self.calculate_thickness_from_period_dynamic( period, theta_deg, avg_n ) overall_thicknesses.append(thickness) # 结果汇总 results = { 'filename': filename, 'theta_deg': theta_deg, 'band_results': band_results, 'overall_thicknesses': overall_thicknesses, 'uncertainty': uncertainty_result, 'wavenumbers': wavenumbers, 'reflectances': reflectances, 'avg_refractive_index': avg_n } print(f"\n" + "=" * 60) print("结果汇总:") print("=" * 60) if uncertainty_result: print(f"多波段加权平均: {uncertainty_result['weighted_mean']:.3f} ± {uncertainty_result['std_dev']:.3f} μm") print(f"相对不确定度: {uncertainty_result['relative_uncertainty']:.2f}%") if overall_thicknesses: print(f"整体分析结果: {np.mean(overall_thicknesses):.3f} μm") print(f"平均折射率: {avg_n:.3f}") # 绘图 self.plot_enhanced_results(results) return results def plot_enhanced_results(self, results): """绘制增强结果""" fig = plt.figure(figsize=(16, 12)) wavenumbers = results['wavenumbers'] reflectances = results['reflectances'] wavelengths = 10000 / wavenumbers # 1. 原始光谱 (双坐标轴) ax1 = plt.subplot(3, 3, 1) ax1.plot(wavenumbers, reflectances, 'b-', linewidth=1) ax1.set_xlabel('波数 (cm^-1)') ax1.set_ylabel('反射率 (%)', color='b') ax1.set_title(f'{results["filename"]} - 原始光谱') ax1.grid(True, alpha=0.3) # 波长坐标轴 ax1_twin = ax1.twiny() ax1_twin.plot(wavelengths, reflectances, 'r-', alpha=0.6) ax1_twin.set_xlabel('波长 (μm)', color='r') ax1_twin.tick_params(axis='x', labelcolor='r') # 2. 动态折射率 ax2 = plt.subplot(3, 3, 2) wl_range = np.linspace(wavelengths.min(), wavelengths.max(), 100) n_range = [self.calculate_dynamic_refractive_index(wl) for wl in wl_range] ax2.plot(wl_range, n_range, 'g-', linewidth=2) ax2.set_xlabel('波长 (μm)') ax2.set_ylabel('折射率') ax2.set_title('动态折射率') ax2.grid(True, alpha=0.3) # 3. 多波段厚度结果 ax3 = plt.subplot(3, 3, 3) if results['band_results']: bands = [r['band'] for r in results['band_results']] thicknesses = [r['thickness'] for r in results['band_results']] colors = ['red', 'green', 'blue', 'orange', 'purple'][:len(bands)] bars = ax3.bar(range(len(bands)), thicknesses, color=colors, alpha=0.7) ax3.set_xticks(range(len(bands))) ax3.set_xticklabels(bands, rotation=45) ax3.set_ylabel('厚度 (μm)') ax3.set_title('多波段厚度结果') # 添加数值标签 for bar, thickness in zip(bars, thicknesses): height = bar.get_height() ax3.text(bar.get_x() + bar.get_width() / 2., height + 0.01, f'{thickness:.3f}', ha='center', va='bottom', fontsize=8) # 4. 傅里叶功率谱 ax4 = plt.subplot(3, 3, 4) freqs, power_spectrum = self.enhanced_fourier_analysis(wavenumbers, reflectances)[2:4] ax4.semilogy(1 / freqs[freqs > 0], power_spectrum[freqs > 0]) ax4.set_xlabel('周期 (cm^-1)') ax4.set_ylabel('功率谱密度') ax4.set_title('傅里叶功率谱') ax4.grid(True, alpha=0.3) # 5. 不确定性分析 ax5 = plt.subplot(3, 3, 5) if results['uncertainty']: unc = results['uncertainty'] mean_val = unc['weighted_mean'] std_val = unc['std_dev'] # 误差条图 ax5.errorbar([0], [mean_val], yerr=[std_val], fmt='ro', capsize=10, capthick=2, markersize=8) ax5.set_xlim(-0.5, 0.5) ax5.set_ylabel('厚度 (μm)') ax5.set_title(f'测量不确定度\n{mean_val:.3f} ± {std_val:.3f} μm') ax5.grid(True, alpha=0.3) ax5.set_xticks([]) # 6. 波段功率对比 ax6 = plt.subplot(3, 3, 6) if results['band_results']: wavelength_centers = [r['avg_wavelength'] for r in results['band_results']] powers = [r['power'] for r in results['band_results']] ax6.scatter(wavelength_centers, powers, s=100, alpha=0.7, c=colors[:len(wavelength_centers)]) ax6.set_xlabel('波长 (μm)') ax6.set_ylabel('功率') ax6.set_title('各波段功率分布') ax6.grid(True, alpha=0.3) # 7-9. 分波段光谱 band_plots = [ax for ax in [plt.subplot(3, 3, 7), plt.subplot(3, 3, 8), plt.subplot(3, 3, 9)]] bands = [(2.5, 5.0, "短波"), (5.0, 10.0, "中波"), (10.0, 25.0, "长波")] for i, ((wl_min, wl_max, band_name), ax) in enumerate(zip(bands, band_plots)): band_mask = (wavelengths >= wl_min) & (wavelengths <= wl_max) if np.sum(band_mask) > 0: ax.plot(wavelengths[band_mask], reflectances[band_mask], 'b-', linewidth=1) ax.set_xlabel('波长 (μm)') ax.set_ylabel('反射率 (%)') ax.set_title(f'{band_name}红外段') ax.grid(True, alpha=0.3) plt.tight_layout() plt.show() def analyze_both_files(self): """分析两个文件""" print("增强分析附件1和附件2") print("=" * 60) # 分析附件1 result1 = self.comprehensive_analysis("C:\\Users\\Virgo\\Desktop\\数学建模大赛\\B题资料\\cleaned_碳化硅_10度入射.xlsx", 10) # 分析附件2 result2 = self.comprehensive_analysis("C:\\Users\\Virgo\\Desktop\\数学建模大赛\\B题资料\\cleaned_碳化硅_15度入射.xlsx", 15) # 综合分析 if result1 and result2: print(f"\n" + "=" * 80) print("最终综合结果") print("=" * 80) unc1 = result1.get('uncertainty') unc2 = result2.get('uncertainty') if unc1 and unc2: t1 = unc1['weighted_mean'] std1 = unc1['std_dev'] t2 = unc2['weighted_mean'] std2 = unc2['std_dev'] # 加权平均 (权重与不确定度成反比) w1 = 1 / (std1 ** 2 + 1e-6) w2 = 1 / (std2 ** 2 + 1e-6) final_thickness = (w1 * t1 + w2 * t2) / (w1 + w2) final_uncertainty = 1 / np.sqrt(w1 + w2) print(f"附件1厚度: {t1:.3f} ± {std1:.3f} μm") print(f"附件2厚度: {t2:.3f} ± {std2:.3f} μm") print(f"加权平均厚度: {final_thickness:.3f} ± {final_uncertainty:.3f} μm") thickness_diff = abs(t1 - t2) print(f"厚度差异: {thickness_diff:.3f} μm") print(f"相对差异: {thickness_diff / final_thickness * 100:.2f}%") if thickness_diff / final_thickness < 0.05: print("两个角度结果一致性优秀") elif thickness_diff / final_thickness < 0.1: print("两个角度结果一致性良好") else: print("两个角度结果存在差异,需进一步分析") return result1, result2 def main(): """主程序""" print("问题2:碳化硅外延层厚度确定算法 (增强版)") print("=" * 80) analyzer = EnhancedInterferenceAnalyzer() result1, result2 = analyzer.analyze_both_files() print("\n增强算法完成!") if __name__ == "__main__": main()(对代码进行全面的优化和检查,得出一个全新的代码,)
09-06
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值