clear; clc; close all

%% 参数区
file1= '附件3.xlsx';
file2 = '附件4.xlsx';
n1    = 3.416;              % 仅用于初步近似
angles_degs = [10, 15];
choose_fit_angles = [10, 15]; % 全谱拟合角度

% 预处理
sgolay_order = 3;
sgolay_win_cfg   = [];
prominence   = 0.5;
min_pk_dist_auto_ratio = 0.05;

% IRLS 稳健拟合
irls_loop = 100;      % 多数函数3-5步即可
tukey_c = 4.653;

% 厚度边界（cm）
d_lb = 1e-6; d_ub = 1e-3;          % 0.01–10 μm
A_lb = -15; A_ub =  15;            % 振幅允许正负，降低与相位的耦合

% ------------------ 物理模型开关与参数 ------------------
USE_FRESNEL_PHYSICS = true;   % 打开物理 Fresnel + 相位 + 色散/Drude + s/p 平均
POLARIZATION_AVG    = true;     % 无偏振平均（默认 true）

% 外延层（膜层）与衬底的基体色散（以 Cauchy 为例，可替换为 Sellmeier）
% n_base(λ) = A + B/λ^2 + C/λ^4,  λ 单位 μm；Im(n_base) 默认 0（K≈0）
cauchy_layer.A = 3.416; cauchy_layer.B = 0.18;  cauchy_layer.C = 0.014; % 硅晶圆片
% cauchy_layer.A = 2.554; cauchy_layer.B = -0.035;  cauchy_layer.C = 0.00; % 碳化硅晶圆片
cauchy_sub  .A = 3.416; cauchy_sub  .B = 0.18;  cauchy_sub  .C = 0.014; % 硅晶圆片
% cauchy_sub  .A = 2.554; cauchy_sub  .B = -0.035;  cauchy_sub  .C = 0.00; % 碳化硅晶圆片

% 自由载流子（Drude）参数（若不考虑，置 N=0 或 add_drude=false）
fc_layer.add_drude = true;              % 外延层是否加 Drude
fc_layer.N_cm3     = 3*10^16;                 % 掺杂浓度（cm^-3）
fc_layer.mstar_m0  = 0.42;              % 有效质量（m*/m0）
fc_layer.mu_cm2Vs  = 100;               % 迁移率（cm^2/Vs）

fc_sub.add_drude   = true;              % 衬底 Drude；常规 K≈0
fc_sub.N_cm3       = 10^18;
fc_sub.mstar_m0    = 0.42;
fc_sub.mu_cm2Vs    = 100;

% 相位微调（额外常量相位，吸收入 e^{2i(delta) + i*phi}），用于小的系统相位偏差
phi_init = 0;        % 初值
phi_lb   = -pi;      % 边界
phi_ub   =  pi;

%% ------------------ 读取数据 ------------------
T1 = readtable(file1, 'ReadVariableNames', false);
T2 = readtable(file2, 'ReadVariableNames', false);
T1.Properties.VariableNames = {'bs','fsl'};
T2.Properties.VariableNames = {'bs','fsl'};

% 按波数排序并去重
T1 = sortrows(T1, 'bs'); T1 = unique(T1, 'rows');
T2 = sortrows(T2, 'bs'); T2 = unique(T2, 'rows');

datasets = {T1, T2};

%% ------------------ 主流程：极值差值法评估厚度 + 稳健全谱拟合 ------------------
all_results = struct([]);

for kk = 1:numel(datasets)
    rushe_deg = angles_degs(kk);
    tab = datasets{kk};

    % 归一化以及预处理
    x = tab.bs(:);            % 单位为cm^-1
    y = tab.fsl(:) / 100.0;   % 分布：0–1
    x = fillmissing(x,'linear'); y = fillmissing(y,'linear');

    n = numel(x);
    sgolay_win = sgolay_win_cfg;
    if isempty(sgolay_win)
        sgolay_win = min(51, max(11, 2*floor(n/10)+1)); % 满足奇数
    end
    y_s = sgolayfilt(y, sgolay_order, sgolay_win);      % 实现平滑操作

    is_out = isoutlier(y_s,'movmedian',max(11, floor(sgolay_win/2)));
    y_s(is_out) = NaN; y_s = fillmissing(y_s,'linear','EndValues','nearest');   %剔除并线性补

    % ---------- 以自动蜂距峰谷检测 ----------
    min_pk_dist = max(5, round(n * min_pk_dist_auto_ratio));
    [~, peaks]   = findpeaks(y_s,  'MinPeakDistance',min_pk_dist, 'MinPeakProminence',prominence/100);
    [~, valleys] = findpeaks(-y_s, 'MinPeakDistance',min_pk_dist, 'MinPeakProminence',prominence/100);

    % ---------- 极值法厚度 ----------
    d_list_um = [];
    delta_nu_peaks = [];  % 该参数用于存储相邻峰之间的波数差
    delta_nu_valleys = []; % 该参数用于存储相邻谷之间的波数差
    
    if numel(peaks) >= 2
        [d_peaks_um, d_nu_peaks] = boshu_distance(peaks, x, n1, rushe_deg);
        d_list_um = [d_list_um; d_peaks_um];
        delta_nu_peaks = [delta_nu_peaks; d_nu_peaks];
    end
    if numel(valleys) >= 2
        [d_valleys_um, d_nu_valleys] = boshu_distance(valleys, x, n1, rushe_deg);
        d_list_um = [d_list_um; d_valleys_um];
        delta_nu_valleys = [delta_nu_valleys; d_nu_valleys];
    end
    
    % 合并所有相邻极值点之间的距离，
    % 其意义是：综合所有可用的干涉条纹信息来获得更可靠的厚度初值估计。
    all_delta_nu = [delta_nu_peaks; delta_nu_valleys];
    
    d0_um = median(d_list_um,'omitnan');
    
    % 计算厚度d的不确定度
    if numel(d_list_um) > 1
        d_std_um = std(d_list_um);
        d_mean_um = mean(d_list_um);
        d_uncertainty_um = d_std_um / sqrt(numel(d_list_um)); % 标准误差
    else
        d_std_um = NaN;
        d_mean_um = mean(d_list_um, 'omitnan');
        d_uncertainty_um = NaN;
    end

    % 考虑厚度无效情况：若极值不足或 d0 异常，使用 FFT 估计周期给出 d0，返回。
    if isempty(d0_um) || ~isfinite(d0_um) || d0_um<=0
        d0_um = fft_um(x, y_s, n1, rushe_deg);
    end
    if ~isfinite(d0_um) || d0_um<=0
        % 给一个更加moderate初值，不妨假设为1 μm
        d0_um = 1.0;
    end
    d0_cm = d0_um * 1e-4;

    % 指定角度做全谱拟合
    will_fit = ismember(rushe_deg, choose_fit_angles);
    fit_car = struct();         %本质上是数据包

    if will_fit
        % 取中位数作为厚度
        d_median_um_loc = median(d_list_um,'omitnan');
        has_valid_median = isfinite(d_median_um_loc) && d_median_um_loc>0;
        if has_valid_median
            d_med_cm = d_median_um_loc * 1e-4;
        end

        % 使用一次背景模型
        model_order = 1;
        [fit_lin, ok_lin] = robust_full_fit(x, y_s, n1, rushe_deg, d0_cm, ...
                                            model_order, tukey_c, irls_loop, ...
                                            d_lb, d_ub, A_lb, A_ub, ...
                                            USE_FRESNEL_PHYSICS, POLARIZATION_AVG, ...
                                            cauchy_layer, fc_layer, cauchy_sub, fc_sub, ...
                                            phi_init, phi_lb, phi_ub);

        fit_car = fit_lin; fit_car.model = 'linear bg + osc/Poly4';
        fit_car.selection = struct();
        fit_car.selection.ok_linear   = ok_lin;
        fit_car.selection.d_median_um = d_median_um_loc;
        if has_valid_median
            fit_car.selection.delta_lin_um = abs(fit_lin.d_cm - d_med_cm)*1e4;
        else
            fit_car.selection.delta_lin_um = NaN;
        end

        % 计算全谱拟合的不确定度
        [fit_uncertainty, fit_error_metrics] = cal_unct(x, y_s, fit_car);
        fit_car.uncertainty = fit_uncertainty;
        fit_car.error_metrics = fit_error_metrics;


        if has_valid_median
            fprintf('角度 %d°：模型 %s；|d_fit-d_median| (μm)=%.3f\n', ...
                rushe_deg, fit_car.model, fit_car.selection.delta_lin_um);
            fprintf('  d敏感度: rel=%.2e, abs=%.2e, dd=%.3gcm; IRLS: %s 于第%d轮\n', ...
                fit_car.sens_d_rel, fit_car.sens_d_abs, fit_car.dd_test, ...
                ternary(fit_car.irls_converged,'收敛','未收敛'), fit_car.irls_iters_run);
            fprintf('  停止原因: %s\n', fit_car.irls_stop_reason);
            fprintf('  拟合不确定度: %.3f μm (%.2f%%)\n', ...
                fit_uncertainty.d_um_uncertainty, fit_uncertainty.d_relative_uncertainty*100);
        else
            fprintf('角度 %d°：模型 %s\n', rushe_deg, fit_car.model);
        end
    end

    % ---------- 汇总 ----------
    res.rushe_deg = rushe_deg;
    res.d_list_um = d_list_um;
    res.d_mean_um = d_mean_um;
    res.d_std_um = d_std_um;
    res.d_uncertainty_um = d_uncertainty_um;
    res.all_delta_nu = all_delta_nu;  % res.all_delta_nu为所有相邻极值点之间的波数差，这和上述波数差作用不完全相同
    res.d0_um     = d0_um;

    if will_fit
        res.fit_model     = fit_car.model;
        res.fit_params    = fit_car.params;
        res.fit_d_um      = fit_car.d_cm * 1e4;
        res.rmse          = fit_car.rmse;
        res.rmse_tail     = fit_car.rmse_tail;
        res.mse           = fit_car.rmse^2;
        res.mse_tail      = fit_car.rmse_tail^2;
        res.fit_uncertainty = fit_car.uncertainty;
        res.error_metrics = fit_car.error_metrics;

        figure('Position',[100,100,1200,1000]);
        
        % 步骤1，先完成：原始/平滑 + 极值
        subplot(5,1,1);
        plot(x, y*100, 'b-', 'LineWidth',1.2, 'DisplayName','原始(%)'); hold on;
        plot(x, y_s*100, 'k-', 'LineWidth',1.2, 'DisplayName','平滑(%)');
        if ~isempty(peaks),   scatter(x(peaks),   y_s(peaks)*100,   25,'ro','DisplayName','极大值'); end
        if ~isempty(valleys), scatter(x(valleys), y_s(valleys)*100, 25,'gx','DisplayName','极小值'); end
        ylabel('反射率(%)'); legend('Location','best'); grid on;
        title(sprintf('角度 %d°：平滑与极值',rushe_deg));

        % 步骤2：拟合（归一）
        subplot(5,1,2);
        plot(x, y_s, 'k-', 'LineWidth',1.2,'DisplayName','数据(归一)');
        hold on; plot(x, fit_car.full_fit, 'r--','LineWidth',1.5,'DisplayName','拟合');
        if ~isempty(fit_car.bg_fit)
            plot(x, fit_car.bg_fit, 'g:','LineWidth',1.2,'DisplayName','背景');
        end
        xlabel('波数(cm^{-1})'); ylabel('反射率(归一)'); legend('Location','best'); grid on;
        title(sprintf('全谱拟合（%s）', res.fit_model));

        % 步骤3： 完善残差分析
        subplot(5,1,3);
        plot(x, fit_car.residual, 'b-', 'LineWidth',1.2);
        hold on; 
        yline(0,'r--');
        yline(mean(fit_car.residual), 'g--', 'LineWidth',1.5, 'DisplayName','残差均值');
        yline(mean(fit_car.residual) + std(fit_car.residual), 'm:', 'LineWidth',1.5, 'DisplayName','±1σ');
        yline(mean(fit_car.residual) - std(fit_car.residual), 'm:', 'LineWidth',1.5);
        xlabel('波数(cm^{-1})'); ylabel('残差'); 
        legend('残差', '零线', '残差均值', 'Location','best'); 
        grid on;
        title(sprintf('RMSE=%.4f, 尾段RMSE=%.4f, 残差标准差=%.4f', res.rmse, res.rmse_tail, std(fit_car.residual)));

        % 步骤4：更新稳健权重
        subplot(5,1,4);
        plot(x, fit_car.final_weights, 'm-','LineWidth',1.5);
        xlabel('波数(cm^{-1})'); ylabel('权重'); grid on;
        title(sprintf('最终权重（稳健）- 平均权重=%.3f', mean(fit_car.final_weights)));

        % 步骤5：画出残差分布直方图
        subplot(5,1,5);
        histogram(fit_car.residual, 30, 'Normalization','pdf', 'FaceColor','blue', 'FaceAlpha',0.7);
        hold on;
        x_values = linspace(min(fit_car.residual), max(fit_car.residual), 100);
        pdf_normal = normpdf(x_values, mean(fit_car.residual), std(fit_car.residual));
        plot(x_values, pdf_normal, 'r-', 'LineWidth',2, 'DisplayName','正态分布');
        xlabel('残差'); ylabel('概率密度'); 
        title(sprintf('残差分布 (偏度=%.3f, 峰度=%.3f)', ...
            res.error_metrics.residual_skewness, res.error_metrics.residual_kurtosis));
        legend('Location','best');
        grid on;

        saveas(gcf, sprintf('角度%d_厚度拟合结果.png', rushe_deg));
    end

    all_results = [all_results; res]; %#ok<AGROW>
end

%% ------------------ 文本部分输出 ------------------
fid = fopen('厚度计算_鲁棒版结果.txt','w','n','UTF-8');
for ii=1:numel(all_results)
    R = all_results(ii);
    fprintf(fid, '角度 %d°：\n', R.rushe_deg);
    fprintf(fid, '  极值法平均厚度：%.3f μm（初值(中位数) %.3f μm）\n', R.d_mean_um, R.d0_um);
    fprintf(fid, '  厚度标准差：%.3f μm\n', R.d_std_um);
    fprintf(fid, '  厚度不确定度（标准误差）：%.3f μm\n', R.d_uncertainty_um);
    fprintf(fid, '  相邻极值点波数差（cm⁻¹）:\n');
    for jj = 1:numel(R.all_delta_nu)
        fprintf(fid, '    Δν%d = %.3f\n', jj, R.all_delta_nu(jj));
    end
    
    if isfield(R,'fit_d_um')
        fprintf(fid, '  全谱拟合模型：%s\n', R.fit_model);
        fprintf(fid, '  全谱拟合厚度：%.3f μm\n', R.fit_d_um);
        fprintf(fid, '  RMSE=%.4f, 尾段RMSE=%.4f\n', R.rmse, R.rmse_tail);
        
        % 以下为拟合误差指标
        fprintf(fid, '  拟合误差指标：\n');
        fprintf(fid, '    残差均值：%.4f\n', R.error_metrics.residual_mean);
        fprintf(fid, '    残差标准差：%.4f\n', R.error_metrics.residual_std);
        fprintf(fid, '    残差偏度：%.3f\n', R.error_metrics.residual_skewness);
        fprintf(fid, '    残差峰度：%.3f\n', R.error_metrics.residual_kurtosis);
        fprintf(fid, '    权重均值：%.3f\n', R.error_metrics.weight_mean);
        fprintf(fid, '    权重标准差：%.3f\n', R.error_metrics.weight_std);
        fprintf(fid, '    低权重比例：%.1f%%\n', R.error_metrics.low_weight_ratio*100);
        
        % 以下为拟合不确定度
        fprintf(fid, '  拟合不确定度：\n');
        fprintf(fid, '    厚度不确定度：%.3f μm\n', R.fit_uncertainty.d_um_uncertainty);
        fprintf(fid, '    相对不确定度：%.2f%%\n', R.fit_uncertainty.d_relative_uncertainty*100);
        fprintf(fid, '    误差传播贡献：%.3f μm\n', R.fit_uncertainty.error_propagation_uncertainty);
        fprintf(fid, '    模型误差贡献：%.3f μm\n', R.fit_uncertainty.model_error_uncertainty);
        fprintf(fid, '    权重分布贡献：%.3f μm\n', R.fit_uncertainty.weight_distribution_uncertainty);
        fprintf(fid, '    置信区间(95%%): [%.3f, %.3f] μm\n', ...
            R.fit_uncertainty.confidence_interval(1), R.fit_uncertainty.confidence_interval(2));
    end
    fprintf(fid, '\n');
end
fclose(fid);
fprintf('\n结果已保存：厚度计算_鲁棒版结果.txt 以及各角度对应的 PNG 图。\n');

%% ================== 函数区 ==================

function [d_um, delta_nu] = boshu_distance(ext, nu, n1_loc, rushe_deg)
% 相邻极值波数间距，我们规定以微米为单位，用于初值估计
% 波束以cm^-1为单位，后面不再提示！
    rushe = rushe_deg*pi/180;
    cos_rushe_t = sqrt(1 - (sin(rushe)/n1_loc)^2);
    delta_nu = diff(nu(ext));         
    d_um = 1 ./ (2*n1_loc*cos_rushe_t*delta_nu) * 1e4;
    % 这里移除无效厚度
    valid_idx = isfinite(d_um) & (d_um > 0);
    d_um = d_um(valid_idx);
    delta_nu = delta_nu(valid_idx);
end

function d0_um = fft_um(x, y, n1_loc, rushe_deg)
% 用 FFT 估计振荡主频，换算膜厚
    x = x(:); y = y(:);
    dx = median(diff(x));
    y = y - sgolayfilt(y, 3, min(101, max(21, 2*floor(numel(y)/8)+1))); % 去趋势
    y = y - mean(y);
    N = numel(y);
    Y = abs(fft(y));
    Y(1:3) = 0;                 % 去近直流，避免直流分量主导FFT结果；突出振荡信号
    [~, idx] = max(Y(1:floor(N/2)));
    f = (idx-1) / (N*dx);       % 周期/单位波数
    if f<=0
        d0_um = NaN; return;
    end
    % 对于 cos(4π n1 d cosθ * x)，x 的周期 Δν = 1/(2 n1 d cosθ)
    theta = rushe_deg*pi/180;
    cth = sqrt(1 - (sin(theta)/n1_loc)^2);
    d0_um = (f) / (2*n1_loc*cth) * 1e4;   %即数学公式： d ≈ f/(2 n1 cosθ_t)
    if ~isfinite(d0_um) || d0_um<=0
        d0_um = NaN;
    end
end

function [uncertainty, error_metrics] = cal_unct(x, y, fit_car)
% 计算全谱拟合的不确定度
    residuals = fit_car.residual;
    weights = fit_car.final_weights;
    d_fit_um = fit_car.d_cm * 1e4;
    
    % 计算误差指标
    error_metrics.residual_mean = mean(residuals);
    error_metrics.residual_std = std(residuals);
    error_metrics.residual_skewness = skewness(residuals);
    error_metrics.residual_kurtosis = kurtosis(residuals);
    error_metrics.weight_mean = mean(weights);
    error_metrics.weight_std = std(weights);
    error_metrics.low_weight_ratio = mean(weights < 0.5); % 低权重比例
    
    %------以下都是误差分析中不确定度-----
    % 1. RMSE的误差传播不确定度
    n_params = 10; % 模型参数共10个
    n_points = length(x);
    dof = n_points - n_params; % 自由度
    
    if dof > 0
        rmse_uncertainty = fit_car.rmse / sqrt(dof);
    else
        rmse_uncertainty = fit_car.rmse;
    end
    
    % 2. 残差分布的不确定度
    residual_uncertainty = error_metrics.residual_std;
    
    % 3. 权重分布的不确定度
    weight_uncertainty = (1 - error_metrics.weight_mean) * fit_car.rmse;
    
    % 4. 模型误差不确定度，这里需要考虑残差的非正态性
    non_normality_factor = 1 + abs(error_metrics.residual_skewness) + max(0, error_metrics.residual_kurtosis - 3)/10;
    model_error_uncertainty = fit_car.rmse * non_normality_factor;
    
    % 综合不确定度（各项的均方根）
    combined_uncertainty = sqrt(rmse_uncertainty^2 + residual_uncertainty^2 + ...
                               weight_uncertainty^2 + model_error_uncertainty^2);
    
    % 转换为厚度不确定度
    % 这里使用经验系数，实际应用中可能需要根据具体问题调整，具体问题具体分析
    d_uncertainty_um = combined_uncertainty * 0.1 * d_fit_um;
    
    % 计算95%置信区间
    confidence_level = 0.95;
    t_value = tinv(1 - (1-confidence_level)/2, dof);
    if ~isfinite(t_value) || dof <= 0
        t_value = 1.96; % 使用正态分布近似
    end
    ci_lower = d_fit_um - t_value * d_uncertainty_um;
    ci_upper = d_fit_um + t_value * d_uncertainty_um;
    
    uncertainty.d_um_uncertainty = d_uncertainty_um;
    uncertainty.d_relative_uncertainty = d_uncertainty_um / d_fit_um;
    uncertainty.error_propagation_uncertainty = rmse_uncertainty * 0.1 * d_fit_um;
    uncertainty.model_error_uncertainty = model_error_uncertainty * 0.1 * d_fit_um;
    uncertainty.weight_distribution_uncertainty = weight_uncertainty * 0.1 * d_fit_um;
    uncertainty.confidence_interval = [ci_lower, ci_upper];
end

function [out, ok] = robust_full_fit(x, y, n1_const, rushe_deg, d0_cm, ...
                                     model_order, tukey_c, irls_iters, ...
                                     d_lb, d_ub, A_lb, A_ub, ...
                                     USE_FRESNEL, POLARIZATION_AVG, ...
                                     cauchy_layer, fc_layer, cauchy_sub, fc_sub, ...
                                     phi_init, phi_lb, phi_ub)
% 稳健全谱拟合 + 回退判据 + d敏感度快检 + 更严格的IRLS收敛检查（含 Poly4 分母）
    x = x(:); y = y(:); n = numel(x);
    ok = false;

    % 背景基函数
    x_mu = mean(x); x_std = std(x); if x_std==0, x_std = 1; end
    x_poly = (x - x_mu) / x_std;

    % ===== 参数布局（新增 Q4：q0..q4） =====
    % 使用十个参数来拟合函数 params = [a1, a0, q0..q4, A, d, phi]
    switch model_order
        case 1
            p0 = [0, mean(y), 1, 0, 0, 0, 0, 1.0, d0_cm, phi_init];
            lb = [-1, 0,  0.1,-2,-2,-2,-2, A_lb, d_lb,  phi_lb];
            ub = [ 1, 1, 10.0, 2, 2, 2, 2, A_ub, d_ub,  phi_ub];
            idx = struct('a1',1,'a0',2,'q0',3,'q1',4,'q2',5,'q3',6,'q4',7,'A',8,'d',9,'phi',10);
    end

    % 目标模型
    model_fun = @(p) moxing(x, x_poly, p, n1_const, rushe_deg, model_order, ...
                                USE_FRESNEL, POLARIZATION_AVG, ...
                                cauchy_layer, fc_layer, cauchy_sub, fc_sub, idx);

    % ---- d 敏感度快检（初值附近）----
    sens_d_rel = NaN; sens_d_abs = NaN; dd_test = NaN; d_insensitive = false;
    if USE_FRESNEL
        theta0 = rushe_deg*pi/180;
        dd_test = max(1e-7, 0.02*abs(d0_cm)); if dd_test==0, dd_test = 1e-6; end
        R0    = reflectance_single_layer(x, theta0, d0_cm,        phi_init, POLARIZATION_AVG, cauchy_layer, fc_layer, cauchy_sub, fc_sub);
        Rplus = reflectance_single_layer(x, theta0, d0_cm+dd_test, phi_init, POLARIZATION_AVG, cauchy_layer, fc_layer, cauchy_sub, fc_sub);
        Rminus= reflectance_single_layer(x, theta0, d0_cm-dd_test, phi_init, POLARIZATION_AVG, cauchy_layer, fc_layer, cauchy_sub, fc_sub);
        rms_ = @(v) sqrt(mean(abs(v).^2));
        sens_d_abs = rms_(Rplus - Rminus);
        sens_d_rel = sens_d_abs / max(1e-12, rms_(R0));
        d_insensitive = (sens_d_rel < 1e-4) || (sens_d_abs < 1e-6);
    end

    % 初始权重（全 1）
    robust_w = ones(n,1);
    weights  = robust_w;

    % 内层 lsqnonlin 数值选项（更严格）
    typicalX = abs(p0); typicalX(typicalX==0) = 1;
    opts = optimoptions('lsqnonlin','Display','off', ...
                        'MaxIterations', 600, ...
                        'MaxFunctionEvaluations', 6000, ...
                        'FunctionTolerance',    1e-13, ...
                        'StepTolerance',        1e-14, ...
                        'OptimalityTolerance',  1e-12, ...
                        'TypicalX', typicalX, ...
                        'FiniteDifferenceType','forward');

    % ---- IRLS 更严格的收敛判据 ----
    param_tol_rel = 1e-8;      % 参数相对变化
    param_tol_abs = 1e-12;     % 参数绝对变化（防 p≈0 失效）
    weight_tol    = 5e-9;      % 权重最大变化
    obj_tol_rel   = 5e-13;     % 目标相对改进

    % 厚度 d 专属阈值（更严）
    d_abs_tol = 5e-8;          % cm（≈0.5 nm）
    d_rel_tol = 1e-7;

    stable_needed = 2;         % 需要连续满足的轮数
    stable_count  = 0;

    p = p0;
    J_prev = inf;
    p_prev = p;
    w_prev = weights;

    irls_converged = false;
    stop_reason = 'max_iters';

    for tt = 1:irls_iters
        % 内层：固定权重的加权最小二乘
        sqrt_w = sqrt(max(weights, eps));
        res_fun = @(pp) sqrt_w .* (model_fun(pp) - y);
        try
            p = lsqnonlin(res_fun, p, lb, ub, opts);
        catch ME
            stop_reason = ['lsqnonlin failure: ' ME.message];
            break;
        end

        % 残差与新权重（仅稳健）
        r = model_fun(p) - y;
        robust_w = tukey_w(r, tukey_c);
        weights  = robust_w;

        % 目标
        J = sum(weights .* (model_fun(p) - y).^2);

        % 变化量
        dp = p - p_prev;
        rel_p_change   = norm(dp) / max(1e-12, norm(p_prev));
        abs_p_change   = norm(dp, Inf);
        max_w_change   = max(abs(weights - w_prev));
        rel_obj_improv = (J_prev - J) / max(1e-12, J_prev);

        % d 专属变化
        d_change_abs = abs(p(idx.d) - p_prev(idx.d));
        d_change_rel = d_change_abs / max(1e-16, abs(p_prev(idx.d)));

        % 是否"稳定一轮"
        param_stable = (rel_p_change < param_tol_rel) || (abs_p_change < param_tol_abs);
        weight_stable= (max_w_change < weight_tol);
        obj_stable   = (rel_obj_improv < obj_tol_rel);
        d_stable     = (d_change_abs < d_abs_tol) || (d_change_rel < d_rel_tol);

        if (param_stable && weight_stable && d_stable) || obj_stable
            stable_count = stable_count + 1;
        else
            stable_count = 0;
        end

        if stable_count >= stable_needed
            irls_converged = true;
            if obj_stable && ~(param_stable && weight_stable && d_stable) == 1
                stop_reason = 'obj_improv_small';
            else
                stop_reason = 'param_weight_d_stable';
            end
            break;
        end

        % 次轮准备
        p_prev = p; w_prev = weights; J_prev = J;
    end

    % 输出
    yhat = model_fun(p);
    resid = y - yhat;
    rmse = sqrt(mean(resid.^2));
    n_last = max(1, floor(n/3));
    rmse_tail = sqrt(mean(resid(end-n_last+1:end).^2));

    [bg_fit, osc_fit] = split_osc(x, x_poly, p, n1_const, rushe_deg, model_order, ...
                                     USE_FRESNEL, POLARIZATION_AVG, ...
                                     cauchy_layer, fc_layer, cauchy_sub, fc_sub, idx);

    out.params = p;
    out.full_fit = yhat;
    out.bg_fit = bg_fit;
    out.osc_fit = osc_fit;
    out.residual = resid;
    out.rmse = rmse;
    out.rmse_tail = rmse_tail;
    out.final_weights = weights;
    out.d_cm = p(idx.d);

    % 敏感度 & 收敛 & 目标
    out.sens_d_rel   = sens_d_rel;
    out.sens_d_abs   = sens_d_abs;
    out.dd_test      = dd_test;
    out.d_insensitive= d_insensitive;

    out.irls_converged   = irls_converged;
    out.irls_iters_run   = tt;
    out.irls_stop_reason = stop_reason;
    out.J_final          = sum((weights .* (model_fun(p)-y)).^2);
    
    % RMSE参数小于预设值说明拟合完美
    ok = (rmse <= 0.01) && (rmse_tail <= 0.012);
end

function yhat = moxing(x, x_poly, p, n1_const, rushe_deg, model_order, ...
                           USE_FRESNEL, POLARIZATION_AVG, ...
                           cauchy_layer, fc_layer, cauchy_sub, fc_sub, idx)
% 一次背景 + 余弦(或 Fresnel) / 四次多项式
    theta0 = rushe_deg*pi/180;

    switch model_order
        case 1
            a1 = p(idx.a1); a0 = p(idx.a0);
            q  = [p(idx.q0), p(idx.q1), p(idx.q2), p(idx.q3), p(idx.q4)];
            A  = p(idx.A);  d  = p(idx.d);  phi = p(idx.phi);
            bg = a1*x_poly + a0;
        otherwise
            error('本版仅支持一次背景模型（model_order=1）');
    end

    % 分母四次多项式（在标准化 x 上）
    denom = eval_poly4(x_poly, q);
    denom = max(denom, 1e-4);  % 数值下限，避免除数为零

    if ~USE_FRESNEL
        % 简化余弦模型（常数 n1）
        cth = sqrt(1 - (sin(theta0)/n1_const)^2);
        osc = cos(4*pi*n1_const*d*cth*x + phi);
    else
        % 物理 Fresnel + 相位 + 色散/Drude + s/p 平均都已考虑模型
        osc = reflectance_single_layer(x, theta0, d, phi, POLARIZATION_AVG, ...
                                       cauchy_layer, fc_layer, cauchy_sub, fc_sub);
    end

    yhat = bg + A .* (osc ./ denom);
end

function [bg, osc_fit] = split_osc(x, x_poly, p, n1_const, rushe_deg, model_order, ...
                                      USE_FRESNEL, POLARIZATION_AVG, ...
                                      cauchy_layer, fc_layer, cauchy_sub, fc_sub, idx)
% 用于输出背景/振荡分量供作图
    theta0 = rushe_deg*pi/180;

    switch model_order
        case 1
            a1 = p(idx.a1); a0 = p(idx.a0);
            q  = [p(idx.q0), p(idx.q1), p(idx.q2), p(idx.q3), p(idx.q4)];
            A  = p(idx.A);  d  = p(idx.d);  phi = p(idx.phi);
            bg = a1*x_poly + a0;
        otherwise
            error('本版仅支持一次背景模型（model_order=1）');
    end

    denom = eval_poly4(x_poly, q); denom = max(denom, 1e-4);

    if ~USE_FRESNEL
        cth = sqrt(1 - (sin(theta0)/n1_const)^2);
        osc = cos(4*pi*n1_const*d*cth*x + phi);
    else
        osc = reflectance_single_layer(x, theta0, d, phi, POLARIZATION_AVG, ...
                                       cauchy_layer, fc_layer, cauchy_sub, fc_sub);
    end

    osc_fit = A .* (osc ./ denom);
end

function w = tukey_w(r, c)
% Tukey biweight（基于MAD的尺度），返回 w ，[0,1]
    r = r(:);
    s = 1.4826 * median(abs(r - median(r))); % MAD尺度
    if s<=eps, w = ones(size(r)); return; end
    u = r / (c*s);
    w = (abs(u)<1) .* (1 - u.^2).^2;
end

%% ------------------ 物理子函数 ------------------

function R = reflectance_single_layer(nu_cm1, theta0, d_cm, phi, POLARIZATION_AVG, ...
                                      cauchy_layer, fc_layer, cauchy_sub, fc_sub)
% 计算空气/单层/衬底结构在斜入射下的反射率（s/p 平均）
% nu_cm1: 波数(cm^-1) 列向量；d_cm: 厚度(cm)；phi: 额外相位（标量）
    nu_cm1 = nu_cm1(:);
    n0 = 1; % 空气

    % 波长换算
    lambda_um = 1e4 ./ nu_cm1;
    lambda_m  = 1e-2 ./ nu_cm1;

    % 复折射率：外延层、衬底
    n1 = n_complex_cauchy_drude(lambda_um, lambda_m, cauchy_layer, fc_layer);
    n2 = n_complex_cauchy_drude(lambda_um, lambda_m, cauchy_sub,   fc_sub);

    % 各层倾斜角的 cos（复数）
    cos0 = cos(theta0) * ones(size(nu_cm1));
    sin0 = sin(theta0);
    cos1 = sqrt(1 - (n0./n1 .* sin0).^2);
    cos2 = sqrt(1 - (n0./n2 .* sin0).^2);
    % 分支稳定化（确保与入射半空间相连的主值）
    cos1 = sign(real(cos1)).*cos1;
    cos2 = sign(real(cos2)).*cos2;

    % 界面 Fresnel 振幅系数
    r01_s = (n0.*cos0 - n1.*cos1) ./ (n0.*cos0 + n1.*cos1);
    r12_s = (n1.*cos1 - n2.*cos2) ./ (n1.*cos1 + n2.*cos2);

    r01_p = (n1.*cos0 - n0.*cos1) ./ (n1.*cos0 + n0.*cos1);
    r12_p = (n2.*cos1 - n1.*cos2) ./ (n2.*cos1 + n1.*cos2);

    % 往返相位（使用波数更稳健）：δ = 2π * n1 * d * cosθ1 * ν
    delta = 2*pi .* n1 .* d_cm .* cos1 .* nu_cm1;  % 复数允许

    % 额外相位：e^{2iδ + iφ}
    expo = exp(2*1i*delta + 1i*phi);

    % 单层 Airy 公式
    rs = (r01_s + r12_s .* expo) ./ (1 + r01_s .* r12_s .* expo);
    rp = (r01_p + r12_p .* expo) ./ (1 + r01_p .* r12_p .* expo);

    if POLARIZATION_AVG
        R = 0.5*(abs(rs).^2 + abs(rp).^2);
    else
        % 若有极化信息可改为返回 rs 或 rp
        R = abs(rs).^2;
    end
end

function n_c = n_complex_cauchy_drude(lambda_um, lambda_m, cauchy, fc)
% 基体色散（Cauchy） + 自由载流子 Drude → 复折射率
% lambda_um: μm；lambda_m: m
    % 基体：Cauchy（实数部分），可按需扩展出 k_base(λ)
    n_base = cauchy.A + cauchy.B./(lambda_um.^2) + cauchy.C./(lambda_um.^4);
    eps_lat = n_base.^2;                % 介电常数（实数）

    if isfield(fc,'add_drude') && fc.add_drude && fc.N_cm3>0
        % 常数
        eps0 = 8.854187817e-12;         % F/m
        e    = 1.602176634e-19;         % C
        m0   = 9.1093837015e-31;        % kg
        c    = 299792458;               % m/s

        N_m3 = fc.N_cm3 * 1e6;                          % cm^-3 → m^-3
        mstar= fc.mstar_m0 * m0;
        mu   = fc.mu_cm2Vs * 1e-4;                      % cm^2/Vs → m^2/Vs

        omega = 2*pi*c ./ lambda_m;                     % rad/s
        omega_p2 = N_m3 * e^2 / (eps0 * mstar);         % rad^2/s^2
        gamma    = e / (mstar * mu);                    % s^-1

        deps = - omega_p2 ./ (omega.^2 + 1i*gamma.*omega);
        eps_total = eps_lat + deps;
    else
        eps_total = eps_lat;
    end

    n_c = sqrt(eps_total);  % 复数主值
    % 保证正的虚部（吸收非负）
    n_c = (real(n_c)>=0).*n_c + (real(n_c)<0).*(-n_c);
end

%% ------------------ 小工具 ------------------
function val = eval_poly4(xn, q)
% q = [q0 q1 q2 q3 q4]，在标准化 x 上计算 q0+q1*x+q2*x^2+q3*x^3+q4*x^4
    val = q(1) + q(2).*xn + q(3).*xn.^2 + q(4).*xn.^3 + q(5).*xn.^4;
end

function s = ternary(cond, a, b)
    if cond, s = a; else, s = b; end
end