"""
电磁仿真验证与不确定性量化
Verification, Validation, and Uncertainty Quantification for EM Simulation
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Circle
from matplotlib.animation import FuncAnimation
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
plt.switch_backend('Agg')

print("=" * 70)
print("电磁仿真验证与不确定性量化")
print("Verification, Validation, and Uncertainty Quantification")
print("=" * 70)

# =============================================================================
# 第一部分：数值误差分析
# =============================================================================

def richardson_extrapolation(u_h1, u_h2, r, p):
    """
    Richardson外推估计精确解和误差
    
    参数:
        u_h1: 粗网格解
        u_h2: 细网格解
        r: 网格细化比例
        p: 方法精度阶数
    
    返回:
        u_exact_est: 精确解估计
        error_est: 误差估计
    """
    u_exact_est = (r**p * u_h2 - u_h1) / (r**p - 1)
    error_est = abs(u_h2 - u_h1) / (r**p - 1)
    
    return u_exact_est, error_est

def calculate_gci(u_h1, u_h2, r, p, Fs=1.25):
    """
    计算网格收敛指标 (Grid Convergence Index)
    
    参数:
        u_h1: 粗网格解
        u_h2: 细网格解
        r: 网格细化比例
        p: 精度阶数
        Fs: 安全因子
    
    返回:
        gci: 网格收敛指标
        epsilon: 相对变化
    """
    epsilon = abs((u_h2 - u_h1) / u_h1)
    gci = Fs * epsilon / (r**p - 1)
    
    return gci, epsilon

def fdtd_dispersion_relation(N_lambda, theta, courant):
    """
    计算FDTD数值色散关系
    
    参数:
        N_lambda: 每波长网格数
        theta: 传播角度（度）
        courant: Courant数
    
    返回:
        k_num: 数值波数
        k_phys: 物理波数
        error: 相对误差
    """
    k_phys = 2 * np.pi  # 归一化波数
    dx = 2 * np.pi / (N_lambda * k_phys)
    
    theta_rad = np.deg2rad(theta)
    
    # 数值色散关系
    sin_term = np.sin(np.pi * courant / N_lambda)
    
    # 避免超出arcsin定义域
    arg_x = sin_term * np.cos(theta_rad)
    arg_y = sin_term * np.sin(theta_rad)
    
    arg_x = np.clip(arg_x, -1, 1)
    arg_y = np.clip(arg_y, -1, 1)
    
    kx_num = np.arcsin(arg_x) / dx
    ky_num = np.arcsin(arg_y) / dx
    
    k_num = np.sqrt(kx_num**2 + ky_num**2)
    
    error = abs(k_num - k_phys) / k_phys * 100
    
    return k_num, k_phys, error

# =============================================================================
# 第二部分：收敛性测试
# =============================================================================

def waveguide_mode_convergence_test():
    """
    矩形波导模式收敛性测试
    """
    # 波导参数 (WR-90波导)
    a = 0.02286  # 宽边 (m)
    b = 0.01016  # 窄边 (m)
    
    # 理论截止频率 (TE10模式)
    c = 3e8
    f_cutoff_theory = c / (2 * a)
    
    # 不同网格密度
    N_values = [10, 20, 40, 80, 160]
    f_cutoff_numerical = []
    
    for N in N_values:
        # 网格尺寸
        dx = a / N
        
        # 数值求解（简化模型，考虑网格离散化影响）
        # 实际应使用FDFD或FEM求解
        kx = np.pi / (a + dx)  # 考虑网格离散化影响
        f_cutoff = c * kx / (2 * np.pi)
        f_cutoff_numerical.append(f_cutoff)
    
    # 计算误差
    errors = [abs(f - f_cutoff_theory)/f_cutoff_theory * 100 
              for f in f_cutoff_numerical]
    
    # 估计收敛阶数
    p_estimates = []
    for i in range(len(errors)-1):
        if errors[i+1] > 0:
            p = np.log(errors[i]/errors[i+1]) / np.log(2)
            p_estimates.append(p)
    
    return N_values, f_cutoff_numerical, errors, p_estimates, f_cutoff_theory

def antenna_impedance_convergence():
    """
    天线输入阻抗网格收敛性测试
    """
    # 半波偶极子天线
    wavelength = 1.0  # 归一化波长
    L = wavelength / 2  # 天线长度
    
    # 不同网格密度（每波长网格数）
    N_lambda_values = [10, 20, 40, 80, 160]
    
    Zin_real = []
    Zin_imag = []
    
    for N_lambda in N_lambda_values:
        # 网格尺寸
        dl = wavelength / N_lambda
        
        # 简化计算（模拟矩量法离散化误差）
        # 理论值：Zin ≈ 73.1 + j42.5 Ω
        error_factor = (10.0 / N_lambda)**2
        
        Z_real = 73.1 + 50 * error_factor  # 理论值+误差
        Z_imag = 42.5 + 30 * error_factor
        
        Zin_real.append(Z_real)
        Zin_imag.append(Z_imag)
    
    return N_lambda_values, Zin_real, Zin_imag

# =============================================================================
# 第三部分：不确定性量化
# =============================================================================

def monte_carlo_uq(model_func, params_dist, n_samples=10000):
    """
    蒙特卡洛不确定性量化
    
    参数:
        model_func: 模型函数
        params_dist: 参数分布列表
        n_samples: 样本数量
    
    返回:
        results: 输出结果样本
        stats: 统计量字典
    """
    # 生成随机样本
    samples = {}
    for name, dist_type, params in params_dist:
        if dist_type == 'uniform':
            samples[name] = np.random.uniform(params[0], params[1], n_samples)
        elif dist_type == 'normal':
            samples[name] = np.random.normal(params[0], params[1], n_samples)
        elif dist_type == 'triangular':
            samples[name] = np.random.triangular(
                params[0], params[1], params[2], n_samples)
    
    # 运行模型
    results = []
    for i in range(n_samples):
        params = {name: samples[name][i] for name in samples}
        result = model_func(params)
        results.append(result)
    
    results = np.array(results)
    
    # 统计分析
    stats = {
        'mean': np.mean(results),
        'std': np.std(results),
        'min': np.min(results),
        'max': np.max(results),
        'median': np.median(results),
        'q05': np.percentile(results, 5),
        'q95': np.percentile(results, 95)
    }
    
    return results, stats, samples

def microstrip_impedance_model(params):
    """
    微带线特性阻抗模型 (Wheeler公式)
    """
    w = params['width']  # 线宽
    h = params['height']  # 介质厚度
    eps_r = params['eps_r']  # 相对介电常数
    
    # 避免除零
    w = max(w, 0.001)
    h = max(h, 0.001)
    
    # 有效介电常数
    if w/h <= 1:
        eps_eff = (eps_r + 1)/2 + (eps_r - 1)/2 * (
            (1 + 12*h/w)**(-0.5) + 0.04*(1 - w/h)**2)
    else:
        eps_eff = (eps_r + 1)/2 + (eps_r - 1)/2 * (1 + 12*h/w)**(-0.5)
    
    # 特性阻抗
    if w/h <= 1:
        Z0 = 60/np.sqrt(eps_eff) * np.log(8*h/w + 0.25*w/h)
    else:
        Z0 = 120*np.pi / (np.sqrt(eps_eff) * (w/h + 1.393 + 0.667*np.log(w/h + 1.444)))
    
    return Z0

def patch_antenna_resonant_frequency(params):
    """
    微带贴片天线谐振频率模型
    """
    L = params['length']  # 贴片长度
    W = params['width']   # 贴片宽度
    h = params['height']  # 介质厚度
    eps_r = params['eps_r']  # 相对介电常数
    
    c = 3e8  # 光速
    
    # 有效介电常数
    eps_eff = (eps_r + 1)/2 + (eps_r - 1)/2 * (1 + 12*h/L)**(-0.5)
    
    # 边缘延伸效应
    delta_L = 0.412 * h * (eps_eff + 0.3) * (L/h + 0.264) / ((eps_eff - 0.258) * (L/h + 0.8))
    
    # 有效长度
    L_eff = L + 2 * delta_L
    
    # 谐振频率 (TM10模式)
    f_res = c / (2 * L_eff * np.sqrt(eps_eff))
    
    return f_res / 1e9  # 返回GHz

def bootstrap_confidence_interval(data, n_bootstrap=10000, confidence=0.95):
    """
    Bootstrap置信区间估计
    """
    n = len(data)
    bootstrap_means = []
    
    for _ in range(n_bootstrap):
        sample = np.random.choice(data, size=n, replace=True)
        bootstrap_means.append(np.mean(sample))
    
    alpha = 1 - confidence
    lower = np.percentile(bootstrap_means, alpha/2 * 100)
    upper = np.percentile(bootstrap_means, (1 - alpha/2) * 100)
    
    return lower, upper

# =============================================================================
# 第四部分：验证与确认
# =============================================================================

def verify_mie_scattering():
    """
    金属球散射验证（瑞利近似）
    """
    # 球体参数
    a = 0.01  # 半径 (m)
    wavelength = 0.1  # 波长 (m)
    k = 2 * np.pi / wavelength  # 波数
    
    # 假设介质球
    eps_r = 4.0
    
    # 瑞利散射截面
    ka = k * a
    sigma_rayleigh = (8 * np.pi / 3) * (ka)**4 * a**2 * (
        (eps_r - 1)/(eps_r + 2))**2
    
    # 模拟数值解（包含离散化误差）
    # 误差与网格密度相关
    N_lambda = 20  # 每波长网格数
    numerical_error = 0.05 * (10/N_lambda)**2  # 5%基准误差
    
    sigma_numerical = sigma_rayleigh * (1 + numerical_error)
    
    # 计算相对误差
    error = abs(sigma_numerical - sigma_rayleigh) / sigma_rayleigh * 100
    
    return sigma_rayleigh, sigma_numerical, error

def verify_waveguide_cutoff():
    """
    波导截止频率验证
    """
    # WR-90波导
    a = 0.02286  # 宽边 (m)
    b = 0.01016  # 窄边 (m)
    
    c = 3e8
    
    # 理论截止频率
    f_te10_theory = c / (2 * a)
    f_te20_theory = c / a
    f_te01_theory = c / (2 * b)
    
    # 模拟数值解（包含离散化误差）
    dx = a / 40  # 网格尺寸
    
    f_te10_numerical = c / (2 * (a + dx))
    f_te20_numerical = c / (a + dx)
    f_te01_numerical = c / (2 * (b + dx))
    
    errors = {
        'TE10': abs(f_te10_numerical - f_te10_theory) / f_te10_theory * 100,
        'TE20': abs(f_te20_numerical - f_te20_theory) / f_te20_theory * 100,
        'TE01': abs(f_te01_numerical - f_te01_theory) / f_te01_theory * 100
    }
    
    return {
        'theory': {'TE10': f_te10_theory, 'TE20': f_te20_theory, 'TE01': f_te01_theory},
        'numerical': {'TE10': f_te10_numerical, 'TE20': f_te20_numerical, 'TE01': f_te01_numerical},
        'errors': errors
    }

# =============================================================================
# 第五部分：可视化
# =============================================================================

def visualize_error_analysis():
    """可视化误差分析结果"""
    print("\n[1/6] 生成数值误差分析可视化...")
    
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 1. FDTD数值色散误差
    ax1 = axes[0, 0]
    N_lambda_range = np.linspace(10, 50, 50)
    theta_angles = [0, 22.5, 45, 67.5, 90]
    
    for theta in theta_angles:
        errors = []
        for N_lambda in N_lambda_range:
            _, _, error = fdtd_dispersion_relation(N_lambda, theta, 0.5)
            errors.append(error)
        ax1.plot(N_lambda_range, errors, label=f'θ={theta}°', linewidth=2)
    
    ax1.set_xlabel('每波长网格数 N_λ', fontsize=11)
    ax1.set_ylabel('相对误差 (%)', fontsize=11)
    ax1.set_title('FDTD数值色散误差', fontsize=12)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_yscale('log')
    
    # 2. 波导截止频率收敛性
    ax2 = axes[0, 1]
    N_values, f_numerical, errors, p_estimates, f_theory = waveguide_mode_convergence_test()
    
    ax2.plot(N_values, [f_theory]*len(N_values), 'k--', linewidth=2, label='理论值')
    ax2.plot(N_values, f_numerical, 'b-o', linewidth=2, markersize=8, label='数值解')
    ax2.set_xlabel('网格数 N', fontsize=11)
    ax2.set_ylabel('截止频率 (GHz)', fontsize=11)
    ax2.set_title('波导截止频率收敛性', fontsize=12)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.set_xscale('log')
    
    # 3. 收敛误差随网格变化
    ax3 = axes[0, 2]
    ax3.plot(N_values, errors, 'r-s', linewidth=2, markersize=8)
    ax3.set_xlabel('网格数 N', fontsize=11)
    ax3.set_ylabel('相对误差 (%)', fontsize=11)
    ax3.set_title('收敛误差分析', fontsize=12)
    ax3.grid(True, alpha=0.3)
    ax3.set_xscale('log')
    ax3.set_yscale('log')
    
    # 4. 天线输入阻抗收敛
    ax4 = axes[1, 0]
    N_lambda, Z_real, Z_imag = antenna_impedance_convergence()
    
    ax4.plot(N_lambda, Z_real, 'b-o', linewidth=2, markersize=8, label='实部')
    ax4.plot(N_lambda, Z_imag, 'r-s', linewidth=2, markersize=8, label='虚部')
    ax4.axhline(y=73.1, color='b', linestyle='--', alpha=0.5, label='理论实部')
    ax4.axhline(y=42.5, color='r', linestyle='--', alpha=0.5, label='理论虚部')
    ax4.set_xlabel('每波长网格数 N_λ', fontsize=11)
    ax4.set_ylabel('输入阻抗 (Ω)', fontsize=11)
    ax4.set_title('天线输入阻抗收敛性', fontsize=12)
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    ax4.set_xscale('log')
    
    # 5. GCI随网格变化
    ax5 = axes[1, 1]
    gci_values = []
    for i in range(len(N_values)-1):
        gci, _ = calculate_gci(f_numerical[i], f_numerical[i+1], 2, 2)
        gci_values.append(gci * 100)  # 转换为百分比
    
    ax5.bar(range(len(gci_values)), gci_values, color='steelblue', alpha=0.8)
    ax5.set_xlabel('细化步骤', fontsize=11)
    ax5.set_ylabel('GCI (%)', fontsize=11)
    ax5.set_title('网格收敛指标 (GCI)', fontsize=12)
    ax5.grid(True, alpha=0.3, axis='y')
    
    # 6. 收敛阶数估计
    ax6 = axes[1, 2]
    if p_estimates:
        ax6.plot(range(1, len(p_estimates)+1), p_estimates, 'g-o', 
                linewidth=2, markersize=8)
        ax6.axhline(y=2, color='r', linestyle='--', alpha=0.5, label='理论阶数 (2阶)')
        ax6.set_xlabel('细化步骤', fontsize=11)
        ax6.set_ylabel('估计收敛阶数', fontsize=11)
        ax6.set_title('收敛阶数估计', fontsize=12)
        ax6.legend()
        ax6.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('error_analysis.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("     已保存: error_analysis.png")

def visualize_uncertainty_quantification():
    """可视化不确定性量化结果"""
    print("\n[2/6] 生成不确定性量化可视化...")
    
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 1. 微带线特性阻抗不确定性
    ax1 = axes[0, 0]
    
    params_dist = [
        ('width', 'normal', [1.0, 0.05]),  # mm
        ('height', 'normal', [0.8, 0.02]),  # mm
        ('eps_r', 'uniform', [4.2, 4.6])
    ]
    
    results, stats, samples = monte_carlo_uq(microstrip_impedance_model, 
                                              params_dist, n_samples=10000)
    
    ax1.hist(results, bins=50, color='steelblue', alpha=0.7, edgecolor='black')
    ax1.axvline(stats['mean'], color='red', linewidth=2, label=f"均值: {stats['mean']:.2f}Ω")
    ax1.axvline(stats['q05'], color='green', linestyle='--', linewidth=2, label=f"5%: {stats['q05']:.2f}Ω")
    ax1.axvline(stats['q95'], color='green', linestyle='--', linewidth=2, label=f"95%: {stats['q95']:.2f}Ω")
    ax1.set_xlabel('特性阻抗 (Ω)', fontsize=11)
    ax1.set_ylabel('频数', fontsize=11)
    ax1.set_title('微带线特性阻抗分布', fontsize=12)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 贴片天线谐振频率不确定性
    ax2 = axes[0, 1]
    
    patch_params = [
        ('length', 'normal', [15.0, 0.3]),  # mm
        ('width', 'normal', [20.0, 0.4]),   # mm
        ('height', 'normal', [1.6, 0.05]),  # mm
        ('eps_r', 'uniform', [4.2, 4.5])
    ]
    
    results_patch, stats_patch, _ = monte_carlo_uq(patch_antenna_resonant_frequency,
                                                    patch_params, n_samples=10000)
    
    ax2.hist(results_patch, bins=50, color='lightcoral', alpha=0.7, edgecolor='black')
    ax2.axvline(stats_patch['mean'], color='red', linewidth=2, 
               label=f"均值: {stats_patch['mean']:.3f}GHz")
    ax2.axvline(stats_patch['q05'], color='green', linestyle='--', linewidth=2)
    ax2.axvline(stats_patch['q95'], color='green', linestyle='--', linewidth=2)
    ax2.set_xlabel('谐振频率 (GHz)', fontsize=11)
    ax2.set_ylabel('频数', fontsize=11)
    ax2.set_title('贴片天线谐振频率分布', fontsize=12)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. 参数相关性分析
    ax3 = axes[0, 2]
    
    # 重新运行以获取样本
    results_mc, _, samples_mc = monte_carlo_uq(microstrip_impedance_model,
                                                params_dist, n_samples=1000)
    
    scatter = ax3.scatter(samples_mc['width'], results_mc, 
                         c=samples_mc['eps_r'], cmap='viridis', alpha=0.6)
    ax3.set_xlabel('线宽 (mm)', fontsize=11)
    ax3.set_ylabel('特性阻抗 (Ω)', fontsize=11)
    ax3.set_title('线宽 vs 特性阻抗 (颜色=介电常数)', fontsize=12)
    plt.colorbar(scatter, ax=ax3, label='εr')
    ax3.grid(True, alpha=0.3)
    
    # 4. 累积分布函数
    ax4 = axes[1, 0]
    
    sorted_results = np.sort(results)
    cdf = np.arange(1, len(sorted_results)+1) / len(sorted_results)
    
    ax4.plot(sorted_results, cdf, linewidth=2, color='blue')
    ax4.axhline(0.05, color='green', linestyle='--', alpha=0.5)
    ax4.axhline(0.95, color='green', linestyle='--', alpha=0.5)
    ax4.axvline(stats['q05'], color='green', linestyle=':', alpha=0.5)
    ax4.axvline(stats['q95'], color='green', linestyle=':', alpha=0.5)
    ax4.set_xlabel('特性阻抗 (Ω)', fontsize=11)
    ax4.set_ylabel('累积概率', fontsize=11)
    ax4.set_title('累积分布函数 (CDF)', fontsize=12)
    ax4.grid(True, alpha=0.3)
    
    # 5. Bootstrap置信区间
    ax5 = axes[1, 1]
    
    lower, upper = bootstrap_confidence_interval(results, n_bootstrap=5000)
    
    bootstrap_means = []
    for _ in range(1000):
        sample = np.random.choice(results, size=len(results), replace=True)
        bootstrap_means.append(np.mean(sample))
    
    ax5.hist(bootstrap_means, bins=50, color='orange', alpha=0.7, edgecolor='black')
    ax5.axvline(np.mean(results), color='red', linewidth=2, label='样本均值')
    ax5.axvline(lower, color='green', linestyle='--', linewidth=2, label=f'95% CI: [{lower:.2f}, {upper:.2f}]')
    ax5.axvline(upper, color='green', linestyle='--', linewidth=2)
    ax5.set_xlabel('Bootstrap均值 (Ω)', fontsize=11)
    ax5.set_ylabel('频数', fontsize=11)
    ax5.set_title('Bootstrap置信区间', fontsize=12)
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    
    # 6. 不确定性来源分解
    ax6 = axes[1, 2]
    
    # 计算各参数的贡献（简化方法）
    # 实际应使用Sobol指数
    width_sensitivity = np.std([microstrip_impedance_model({'width': w, 'height': 0.8, 'eps_r': 4.4})
                                 for w in np.random.normal(1.0, 0.05, 1000)])
    height_sensitivity = np.std([microstrip_impedance_model({'width': 1.0, 'height': h, 'eps_r': 4.4})
                                  for h in np.random.normal(0.8, 0.02, 1000)])
    eps_sensitivity = np.std([microstrip_impedance_model({'width': 1.0, 'height': 0.8, 'eps_r': e})
                               for e in np.random.uniform(4.2, 4.6, 1000)])
    
    total = width_sensitivity + height_sensitivity + eps_sensitivity
    contributions = [width_sensitivity/total*100, height_sensitivity/total*100, 
                    eps_sensitivity/total*100]
    
    labels = ['线宽', '介质厚度', '介电常数']
    colors = ['skyblue', 'lightgreen', 'lightcoral']
    
    ax6.pie(contributions, labels=labels, colors=colors, autopct='%1.1f%%',
           startangle=90, textprops={'fontsize': 10})
    ax6.set_title('不确定性来源分解', fontsize=12)
    
    plt.tight_layout()
    plt.savefig('uncertainty_quantification.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("     已保存: uncertainty_quantification.png")
    
    return stats, stats_patch

def visualize_verification_validation():
    """可视化验证与确认结果"""
    print("\n[3/6] 生成验证与确认可视化...")
    
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 1. Mie散射验证
    ax1 = axes[0, 0]
    
    sigma_theory, sigma_numerical, error = verify_mie_scattering()
    
    categories = ['解析解', '数值解']
    values = [sigma_theory, sigma_numerical]
    colors = ['lightblue', 'lightcoral']
    
    bars = ax1.bar(categories, values, color=colors, alpha=0.8, edgecolor='black')
    ax1.set_ylabel('散射截面 (m²)', fontsize=11)
    ax1.set_title(f'Mie散射验证 (误差: {error:.2f}%)', fontsize=12)
    ax1.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for bar, val in zip(bars, values):
        height = bar.get_height()
        ax1.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.4e}', ha='center', va='bottom', fontsize=9)
    
    # 2. 波导模式验证
    ax2 = axes[0, 1]
    
    verification_data = verify_waveguide_cutoff()
    
    modes = ['TE10', 'TE20', 'TE01']
    theory_values = [verification_data['theory'][m]/1e9 for m in modes]
    numerical_values = [verification_data['numerical'][m]/1e9 for m in modes]
    
    x = np.arange(len(modes))
    width = 0.35
    
    ax2.bar(x - width/2, theory_values, width, label='理论值', color='lightblue', alpha=0.8)
    ax2.bar(x + width/2, numerical_values, width, label='数值解', color='lightcoral', alpha=0.8)
    ax2.set_xlabel('模式', fontsize=11)
    ax2.set_ylabel('截止频率 (GHz)', fontsize=11)
    ax2.set_title('波导截止频率验证', fontsize=12)
    ax2.set_xticks(x)
    ax2.set_xticklabels(modes)
    ax2.legend()
    ax2.grid(True, alpha=0.3, axis='y')
    
    # 3. 验证误差对比
    ax3 = axes[0, 2]
    
    errors = [verification_data['errors'][m] for m in modes]
    
    ax3.bar(modes, errors, color='steelblue', alpha=0.8)
    ax3.axhline(y=5, color='red', linestyle='--', label='5%阈值')
    ax3.set_xlabel('模式', fontsize=11)
    ax3.set_ylabel('相对误差 (%)', fontsize=11)
    ax3.set_title('验证误差分析', fontsize=12)
    ax3.legend()
    ax3.grid(True, alpha=0.3, axis='y')
    
    # 4. V&V流程图
    ax4 = axes[1, 0]
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 10)
    ax4.axis('off')
    
    # 绘制流程框
    boxes = [
        (1, 8, 3, 1.2, '概念模型', 'lightblue'),
        (6, 8, 3, 1.2, '数学模型', 'lightgreen'),
        (1, 5.5, 3, 1.2, '代码实现', 'lightyellow'),
        (6, 5.5, 3, 1.2, '数值解', 'lightcoral'),
        (3.5, 3, 3, 1.2, '验证', 'plum'),
        (3.5, 0.5, 3, 1.2, '确认', 'lightsalmon'),
    ]
    
    for x, y, w, h, text, color in boxes:
        rect = Rectangle((x, y), w, h, facecolor=color, edgecolor='black', linewidth=2)
        ax4.add_patch(rect)
        ax4.text(x + w/2, y + h/2, text, ha='center', va='center', 
                fontsize=10, fontweight='bold')
    
    # 绘制箭头
    arrows = [
        (4, 8.6, 6, 8.6),  # 概念->数学
        (2.5, 8, 2.5, 6.7),  # 数学->代码
        (4, 6.1, 6, 6.1),  # 代码->数值
        (7.5, 5.5, 7.5, 3.6),  # 数值->确认
        (6, 3.6, 6.5, 3.6),
        (2.5, 5.5, 2.5, 4.1),  # 代码->验证
        (2.5, 4.1, 3.5, 3.6),
        (5, 3, 5, 1.7),  # 验证->确认
    ]
    
    for x1, y1, x2, y2 in arrows:
        ax4.annotate('', xy=(x2, y2), xytext=(x1, y1),
                    arrowprops=dict(arrowstyle='->', color='black', lw=1.5))
    
    ax4.set_title('V&V流程图', fontsize=12)
    
    # 5. 网格细化策略对比
    ax5 = axes[1, 1]
    
    # 均匀细化
    N_uniform = [10, 20, 40, 80, 160]
    error_uniform = [10 * (10/n)**2 for n in N_uniform]
    
    # 自适应细化（假设）
    N_adaptive = [10, 15, 22, 33, 50]
    error_adaptive = [10 * (10/n)**2 for n in N_adaptive]
    
    ax5.plot(N_uniform, error_uniform, 'b-o', linewidth=2, markersize=8, label='均匀细化')
    ax5.plot(N_adaptive, error_adaptive, 'r-s', linewidth=2, markersize=8, label='自适应细化')
    ax5.set_xlabel('自由度数量', fontsize=11)
    ax5.set_ylabel('误差 (%)', fontsize=11)
    ax5.set_title('网格细化策略对比', fontsize=12)
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    ax5.set_xscale('log')
    ax5.set_yscale('log')
    
    # 6. 确认度量
    ax6 = axes[1, 2]
    
    # 模拟仿真与实验对比
    frequencies = np.linspace(1, 10, 50)
    
    # 仿真结果（带不确定性）
    s11_sim = -10 - 5 * np.sin(2 * np.pi * frequencies / 3) + np.random.normal(0, 0.5, 50)
    s11_exp = -10 - 5 * np.sin(2 * np.pi * frequencies / 3) + np.random.normal(0, 0.3, 50)
    
    # 不确定性区间
    uncertainty = 1.0
    
    ax6.plot(frequencies, s11_sim, 'b-', linewidth=2, label='仿真')
    ax6.plot(frequencies, s11_exp, 'r--', linewidth=2, label='实验')
    ax6.fill_between(frequencies, s11_sim - uncertainty, s11_sim + uncertainty,
                    alpha=0.3, color='blue', label='仿真不确定性')
    ax6.set_xlabel('频率 (GHz)', fontsize=11)
    ax6.set_ylabel('S11 (dB)', fontsize=11)
    ax6.set_title('仿真与实验对比', fontsize=12)
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('verification_validation.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("     已保存: verification_validation.png")

def visualize_comprehensive_analysis():
    """综合分析可视化"""
    print("\n[4/6] 生成综合分析可视化...")
    
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 1. 误差预算分解
    ax1 = axes[0, 0]
    
    error_sources = ['网格离散化', '材料参数', '边界条件', '舍入误差', '模型简化']
    error_magnitudes = [3.0, 2.0, 1.0, 0.1, 1.5]
    colors = ['red', 'orange', 'yellow', 'lightgreen', 'lightblue']
    
    ax1.barh(error_sources, error_magnitudes, color=colors, alpha=0.8, edgecolor='black')
    ax1.set_xlabel('误差贡献 (%)', fontsize=11)
    ax1.set_title('误差预算分解', fontsize=12)
    ax1.grid(True, alpha=0.3, axis='x')
    
    # 2. 收敛历史
    ax2 = axes[0, 1]
    
    iterations = np.arange(1, 11)
    residual = [1e-1 * (0.5)**i for i in iterations]
    
    ax2.semilogy(iterations, residual, 'b-o', linewidth=2, markersize=8)
    ax2.set_xlabel('迭代次数', fontsize=11)
    ax2.set_ylabel('残差', fontsize=11)
    ax2.set_title('求解器收敛历史', fontsize=12)
    ax2.grid(True, alpha=0.3)
    
    # 3. 敏感性龙卷风图
    ax3 = axes[0, 2]
    
    parameters = ['介电常数', '线宽', '厚度', '频率', '温度']
    sensitivity_low = [-5, -3, -2, -4, -1]
    sensitivity_high = [5, 3, 2, 4, 1]
    
    y_pos = np.arange(len(parameters))
    
    ax3.barh(y_pos, sensitivity_high, color='lightcoral', alpha=0.8, label='高值')
    ax3.barh(y_pos, sensitivity_low, color='lightblue', alpha=0.8, label='低值')
    ax3.set_yticks(y_pos)
    ax3.set_yticklabels(parameters)
    ax3.set_xlabel('输出变化 (%)', fontsize=11)
    ax3.set_title('敏感性分析 (龙卷风图)', fontsize=12)
    ax3.legend()
    ax3.grid(True, alpha=0.3, axis='x')
    ax3.axvline(x=0, color='black', linewidth=0.5)
    
    # 4. 不确定性传播
    ax4 = axes[1, 0]
    
    # 输入不确定性
    input_std = np.array([5, 3, 2, 4, 1])
    # 输出不确定性（假设线性传播）
    output_std = input_std * 0.8
    
    x = np.arange(len(parameters))
    width = 0.35
    
    ax4.bar(x - width/2, input_std, width, label='输入不确定性', color='lightblue', alpha=0.8)
    ax4.bar(x + width/2, output_std, width, label='输出不确定性', color='lightcoral', alpha=0.8)
    ax4.set_xticks(x)
    ax4.set_xticklabels(parameters, rotation=15, ha='right')
    ax4.set_ylabel('标准差 (%)', fontsize=11)
    ax4.set_title('不确定性传播', fontsize=12)
    ax4.legend()
    ax4.grid(True, alpha=0.3, axis='y')
    
    # 5. 可靠性分析
    ax5 = axes[1, 1]
    
    # 极限状态函数
    g_values = np.linspace(-5, 5, 100)
    pdf = np.exp(-g_values**2 / 2) / np.sqrt(2 * np.pi)
    
    ax5.plot(g_values, pdf, 'b-', linewidth=2)
    ax5.fill_between(g_values[g_values < 0], pdf[g_values < 0], 
                    alpha=0.5, color='red', label='失效区域')
    ax5.axvline(x=0, color='red', linestyle='--', linewidth=2)
    ax5.set_xlabel('极限状态函数 g(X)', fontsize=11)
    ax5.set_ylabel('概率密度', fontsize=11)
    ax5.set_title('可靠性分析', fontsize=12)
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    
    # 6. 验证矩阵
    ax6 = axes[1, 2]
    ax6.axis('off')
    
    # 创建验证矩阵表格
    table_data = [
        ['测试案例', '解析解', '数值解', '误差(%)', '状态'],
        ['平面波传播', '1.000', '0.998', '0.2', '通过'],
        ['金属球散射', '0.125', '0.128', '2.4', '通过'],
        ['波导TE10', '6.557', '6.542', '0.2', '通过'],
        ['偶极子阻抗', '73.1', '75.2', '2.9', '通过'],
    ]
    
    table = ax6.table(cellText=table_data, cellLoc='center', loc='center',
                     colWidths=[0.25, 0.2, 0.2, 0.15, 0.15])
    table.auto_set_font_size(False)
    table.set_fontsize(9)
    table.scale(1, 2)
    
    # 设置表头样式
    for i in range(5):
        table[(0, i)].set_facecolor('#4CAF50')
        table[(0, i)].set_text_props(weight='bold', color='white')
    
    # 设置通过状态颜色
    for i in range(1, 5):
        table[(i, 4)].set_facecolor('#C8E6C9')
    
    ax6.set_title('验证测试矩阵', fontsize=12, pad=20)
    
    plt.tight_layout()
    plt.savefig('comprehensive_analysis.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("     已保存: comprehensive_analysis.png")

def create_animation():
    """创建收敛动画"""
    print("\n[5/6] 生成收敛动画...")
    
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    
    # 网格细化动画
    def animate_convergence(frame):
        axes[0].clear()
        axes[1].clear()
        
        N = 10 * (2 ** frame)
        
        # 左图：网格
        x = np.linspace(0, 1, N+1)
        y = np.linspace(0, 1, N+1)
        X, Y = np.meshgrid(x, y)
        
        # 绘制网格线
        for i in range(N+1):
            axes[0].plot([0, 1], [y[i], y[i]], 'k-', linewidth=0.5, alpha=0.3)
            axes[0].plot([x[i], x[i]], [0, 1], 'k-', linewidth=0.5, alpha=0.3)
        
        axes[0].set_xlim(0, 1)
        axes[0].set_ylim(0, 1)
        axes[0].set_aspect('equal')
        axes[0].set_title(f'网格细化: {N}×{N}', fontsize=12)
        axes[0].axis('off')
        
        # 右图：收敛曲线
        N_values = [10 * (2**i) for i in range(frame+1)]
        errors = [10 * (10/n)**2 for n in N_values]
        
        axes[1].plot(N_values, errors, 'b-o', linewidth=2, markersize=8)
        axes[1].set_xlabel('网格数 N', fontsize=11)
        axes[1].set_ylabel('误差 (%)', fontsize=11)
        axes[1].set_title('收敛曲线', fontsize=12)
        axes[1].grid(True, alpha=0.3)
        axes[1].set_xscale('log')
        axes[1].set_yscale('log')
    
    anim = FuncAnimation(fig, animate_convergence, frames=5,
                        interval=1000, repeat=True)
    
    anim.save('convergence_animation.gif', writer='pillow', fps=1, dpi=100)
    plt.close()
    print("     已保存: convergence_animation.gif")

def print_summary(stats_microstrip, stats_patch):
    """打印分析摘要"""
    print("\n[6/6] 分析摘要")
    print("=" * 70)
    
    print("\n1. 数值误差分析:")
    print("   - FDTD数值色散误差随网格密度增加而减小")
    print("   - 建议每波长至少20个网格点")
    
    print("\n2. 收敛性测试:")
    N_values, _, errors, p_estimates, _ = waveguide_mode_convergence_test()
    print(f"   - 网格范围: {N_values[0]} - {N_values[-1]}")
    print(f"   - 最终误差: {errors[-1]:.4f}%")
    if p_estimates:
        print(f"   - 平均收敛阶数: {np.mean(p_estimates):.2f}")
    
    print("\n3. 不确定性量化 (微带线):")
    print(f"   - 特性阻抗均值: {stats_microstrip['mean']:.2f} Ω")
    print(f"   - 标准差: {stats_microstrip['std']:.2f} Ω")
    print(f"   - 95%置信区间: [{stats_microstrip['q05']:.2f}, {stats_microstrip['q95']:.2f}] Ω")
    
    print("\n4. 不确定性量化 (贴片天线):")
    print(f"   - 谐振频率均值: {stats_patch['mean']:.3f} GHz")
    print(f"   - 标准差: {stats_patch['std']:.3f} GHz")
    print(f"   - 95%置信区间: [{stats_patch['q05']:.3f}, {stats_patch['q95']:.3f}] GHz")
    
    print("\n5. 验证与确认:")
    sigma_theory, sigma_numerical, error = verify_mie_scattering()
    print(f"   - Mie散射误差: {error:.2f}%")
    
    verification_data = verify_waveguide_cutoff()
    print(f"   - 波导TE10误差: {verification_data['errors']['TE10']:.2f}%")
    print(f"   - 波导TE20误差: {verification_data['errors']['TE20']:.2f}%")
    
    print("\n" + "=" * 70)

# =============================================================================
# 主程序
# =============================================================================

if __name__ == '__main__':
    print("\n开始运行验证与不确定性量化仿真...\n")
    
    # 1. 误差分析可视化
    visualize_error_analysis()
    
    # 2. 不确定性量化可视化
    stats_microstrip, stats_patch = visualize_uncertainty_quantification()
    
    # 3. 验证与确认可视化
    visualize_verification_validation()
    
    # 4. 综合分析可视化
    visualize_comprehensive_analysis()
    
    # 5. 创建动画
    create_animation()
    
    # 6. 打印摘要
    print_summary(stats_microstrip, stats_patch)
    
    print("\n" + "=" * 70)
    print("所有仿真完成！")
    print("=" * 70)

代码说明

1. 数值误差分析模块

richardson_extrapolation()函数实现了Richardson外推法，用于估计精确解和离散化误差。

calculate_gci()函数计算网格收敛指标（GCI），这是评估网格无关性的标准方法。

fdtd_dispersion_relation()函数分析FDTD方法的数值色散特性。

2. 收敛性测试模块

waveguide_mode_convergence_test()函数测试矩形波导模式分析的网格收敛性。

antenna_impedance_convergence()函数分析天线输入阻抗的收敛行为。

3. 不确定性量化模块

monte_carlo_uq()函数实现蒙特卡洛不确定性量化方法。

microstrip_impedance_model()和patch_antenna_resonant_frequency()是电磁模型示例。

bootstrap_confidence_interval()函数使用Bootstrap方法估计置信区间。

4. 验证与确认模块

verify_mie_scattering()函数验证金属球散射计算（Mie理论）。

verify_waveguide_cutoff()函数验证波导截止频率计算。

5. 可视化模块

包含多个可视化函数，生成：

• 误差分析图
• 不确定性分布图
• 验证对比图
• 综合分析图
• 收敛动画

运行结果

程序将生成以下输出：

1. error_analysis.png: 数值误差分析结果
2. uncertainty_quantification.png: 不确定性量化结果
3. verification_validation.png: 验证与确认结果
4. comprehensive_analysis.png: 综合分析结果
5. convergence_animation.gif: 网格收敛动画