python Chapter 1:实测粒子下落速度与经典值之间的差异大于5m/s的掩膜

 

import numpy as np
import pandas as pd


def calculate_velocity_mask_and_save_excel():
    """
    计算速度差异掩膜（一次性计算，所有分钟共用）
    返回32×32的布尔数组，True表示速度差异>5m/s
    并将掩膜矩阵保存到Excel文件
    行=直径通道，列=速度通道
    """

    # ==================== 第一步：计算经典下落速度 ====================
    # 32个直径通道的中心直径 (mm)
    diameters_center = np.array([
        0.062, 0.187, 0.312, 0.437, 0.562, 0.687, 0.812, 0.937,
        1.062, 1.187, 1.375, 1.625, 1.875, 2.125, 2.375, 2.75,
        3.25, 3.75, 4.25, 4.75, 5.5, 6.5, 7.5, 8.5, 9.5,
        11, 13, 15, 17, 19, 21.5, 24.5
    ])

    # 32个速度通道的中心速度 (m/s)
    measured_velocity = np.array([
        0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95,
        1.1, 1.3, 1.5, 1.7, 1.9, 2.2, 2.6, 3, 3.4, 3.8,
        4.4, 5.2, 6, 6.8, 7.6, 8.8, 10.4, 12, 13.6, 15.2, 17.6, 20.8
    ])

    # 计算经典下落速度 (1×32数组)
    def atlas_velocity(diameters_mm):
        """计算经典下落速度 (m/s)，输入直径单位为mm"""
        return 9.65 - 10.3 * np.exp(-0.6 * diameters_mm)

    classical_velocities = atlas_velocity(diameters_center)  # 1×32

    print("第一步：经典下落速度 (1×32)")
    print(f"直径[11] = {diameters_center[10]:.3f} mm")
    print(f"经典速度[11] = {classical_velocities[10]:.3f} m/s")

    # ==================== 第二步：扩展经典速度网格 ====================
    # 将1×32数组沿行方向复制32次 → 32×32数组
    # classical_grid[i,j] = 直径通道i的经典速度 (行=直径, 列=速度)
    classical_grid = np.tile(classical_velocities.reshape(-1, 1), (1, 32))  # 32×32

    print("\n第二步：经典速度网格 (32×32)")
    print(f"形状: {classical_grid.shape}")
    print(f"第11行第1列: {classical_grid[10, 0]:.3f} m/s (应该是直径11的经典速度)")
    print(f"第1行第11列: {classical_grid[0, 10]:.3f} m/s (应该是直径1的经典速度)")

    # ==================== 第三步：扩展实测速度网格 ====================
    # 将1×32数组沿列方向复制32次 → 32×32数组
    # measured_grid[i,j] = 速度通道j的实测速度 (行=直径, 列=速度)
    measured_grid = np.tile(measured_velocity, (32, 1))  # 32×32

    print("\n第三步：实测速度网格 (32×32)")
    print(f"形状: {measured_grid.shape}")
    print(f"第1行第11列: {measured_grid[0, 10]:.3f} m/s (应该是速度通道11的实测速度)")
    print(f"第11行第1列: {measured_grid[10, 0]:.3f} m/s (应该是速度通道1的实测速度)")

    # ==================== 第四步：计算速度差异 ====================
    # (32×32) - (32×32) = 32×32
    velocity_diff = measured_grid - classical_grid

    print("\n第四步：速度差异计算")
    print(f"第11行第1列差异: {measured_grid[10, 0]:.3f} - {classical_grid[10, 0]:.3f} = {velocity_diff[10, 0]:.3f} m/s")

    # ==================== 第五步：标记异常位置 ====================
    velocity_mask = np.abs(velocity_diff) > 5.0  # 32×32布尔数组

    print("\n第五步：标记异常位置")
    print(f"第11行第1列: |{velocity_diff[10, 0]:.3f}| > 5? {velocity_mask[10, 0]}")

    # 统计信息
    total_cells = velocity_mask.size
    flagged_cells = np.sum(velocity_mask)

    print(f"\n统计信息:")
    print(f"总网格数: {total_cells} (32×32)")
    print(f"应剔除数: {flagged_cells}")
    print(f"剔除比例: {flagged_cells / total_cells * 100:.2f}%")

    # ==================== 将掩膜矩阵写入Excel ====================
    # 创建列名（速度通道1-32）
    columns = [f'速度通道{i + 1}' for i in range(32)]

    # 创建索引名（直径通道1-32）
    index = [f'直径通道{i + 1}' for i in range(32)]

    # 创建DataFrame
    df_mask = pd.DataFrame(velocity_mask.astype(int),
                           columns=columns,
                           index=index)

    # 保存到Excel
    excel_path = "D:/yanmo.xlsx"
    df_mask.to_excel(excel_path, sheet_name='速度差异掩膜')

    print(f"\n速度差异掩膜已保存到: {excel_path}")

    # 验证特定位置
    print("\n验证位置 (第11行, 第1列):")
    print(f"  行索引: 直径通道11 = {diameters_center[10]:.3f} mm")
    print(f"  列索引: 速度通道1 = {measured_velocity[0]:.3f} m/s")
    print(f"  该点经典速度: {classical_velocities[10]:.3f} m/s")
    print(f"  该点实测速度: {measured_velocity[0]:.3f} m/s")
    print(f"  速度差异: {velocity_diff[10, 0]:.3f} m/s")
    print(f"  绝对值: {abs(velocity_diff[10, 0]):.3f} m/s")
    print(f"  是否应剔除: {velocity_mask[10, 0]} (1=是, 0=否)")

    return velocity_mask, diameters_center, measured_velocity, excel_path


# 运行函数
if __name__ == "__main__":
    velocity_mask, diameters, measured_vel, excel_file = calculate_velocity_mask_and_save_excel()