/// MoonMet - MoonBit 气象计算库
/// 基于 MetPy 实现的主要气象计算功能
/// 提供热力学、动力学和基础气象计算函数

// ============================================================================
// 1. 湿空气属性函数
// ============================================================================

///|
/// 计算湿空气的比气体常数 R_m
/// 公式: R_m = (1 - q) R_d + q R_v
/// q: 比湿 (kg/kg)
pub fn moist_air_gas_constant(specific_humidity : Double) -> Double {
  dry_air_gas_constant +
  specific_humidity * (water_gas_constant - dry_air_gas_constant)
}

///|
/// 计算湿空气的定压比热 C_pm
/// 公式: C_pm = (1 - q) C_pd + q C_pv
/// q: 比湿 (kg/kg)
pub fn moist_air_specific_heat_pressure(specific_humidity : Double) -> Double {
  dry_air_spec_heat_press +
  specific_humidity * (wv_specific_heat_press - dry_air_spec_heat_press)
}

///|
/// 计算湿空气的 Poisson 指数 κ_m
/// 公式: κ_m = R_m / C_pm
pub fn moist_air_poisson_exponent(specific_humidity : Double) -> Double {
  moist_air_gas_constant(specific_humidity) /
  moist_air_specific_heat_pressure(specific_humidity)
}

///|
/// 计算水的汽化潜热
/// 公式: L = L_0 - (C_pl - C_pv) * (T - T_0)
/// T: 温度 (K)
pub fn latent_heat_vaporization(temperature : Double) -> Double {
  water_heat_vaporization -
  (water_specific_heat - wv_specific_heat_press) *
  (temperature - water_triple_point_temperature)
}

///|
/// 计算水的升华潜热
/// 公式: L_s = L_v + L_f
pub fn latent_heat_sublimation(temperature : Double) -> Double {
  latent_heat_vaporization(temperature) + water_heat_fusion
}

///|
/// 计算水的融化潜热
pub fn latent_heat_melting(_temperature : Double) -> Double {
  water_heat_fusion
}

// ============================================================================
// 2. 饱和水汽压函数
// ============================================================================

///|
/// 计算液态水的饱和水汽压（内部函数）
/// 公式: e = e_{s0} * (T_0/T)^((C_pl - C_pv)/R_v) * exp((L_0/(R_v*T_0) - L/(R_v*T)))
fn _saturation_vapor_pressure_liquid(temperature : Double) -> Double {
  let latent_heat = latent_heat_vaporization(temperature)
  let heat_power = (water_specific_heat - wv_specific_heat_press) /
    water_gas_constant
  let exp_term = water_heat_vaporization /
    (water_gas_constant * water_triple_point_temperature) -
    latent_heat / (water_gas_constant * temperature)
  sat_pressure_0c *
  Double::pow(water_triple_point_temperature / temperature, heat_power) *
  @math.exp(exp_term)
}

///|
/// 计算冰的饱和水汽压（内部函数）
/// 公式: e_i = e_{i0} * (T_0/T)^((C_pi - C_pv)/R_v) * exp((L_s0/(R_v*T_0) - L_s/(R_v*T)))
fn _saturation_vapor_pressure_solid(temperature : Double) -> Double {
  let latent_heat = latent_heat_sublimation(temperature)
  let ls0 = water_heat_vaporization + water_heat_fusion
  let heat_power = (ice_specific_heat - wv_specific_heat_press) /
    water_gas_constant
  let exp_term = ls0 / (water_gas_constant * water_triple_point_temperature) -
    latent_heat / (water_gas_constant * temperature)
  sat_pressure_0c *
  Double::pow(water_triple_point_temperature / temperature, heat_power) *
  @math.exp(exp_term)
}

///|
/// 计算饱和水汽压
/// temperature: 温度 (K)
/// phase: "liquid" (默认), "solid", 或 "auto"
pub fn saturation_vapor_pressure(
  temperature : Double,
  phase : String,
) -> Double {
  match phase {
    "liquid" => _saturation_vapor_pressure_liquid(temperature)
    "solid" => _saturation_vapor_pressure_solid(temperature)
    "auto" =>
      if temperature > water_triple_point_temperature {
        _saturation_vapor_pressure_liquid(temperature)
      } else {
        _saturation_vapor_pressure_solid(temperature)
      }
    _ => _saturation_vapor_pressure_liquid(temperature) // 默认使用液态
  }
}

// ============================================================================
// 3. Exner 函数和位温
// ============================================================================

///|
/// 计算 Exner 函数
/// 公式: Π = (p/p_0)^κ
/// pressure: 气压 (Pa)
/// reference_pressure: 参考气压 (Pa)，默认 100000 Pa
pub fn exner_function(
  pressure : Double,
  reference_pressure? : Double,
) -> Double {
  let ref_p = match reference_pressure {
    Some(p) => p
    None => pot_temp_ref_press
  }
  Double::pow(pressure / ref_p, poisson_exponent)
}

///|
/// 计算位温
/// 公式: θ = T * (p_0/p)^κ
/// pressure: 气压 (Pa)
/// temperature: 温度 (K)
pub fn potential_temperature(pressure : Double, temperature : Double) -> Double {
  let ref_p : Double = pot_temp_ref_press
  temperature / exner_function(pressure, reference_pressure=ref_p)
}

///|
/// 从位温计算温度
/// 公式: T = θ * (p/p_0)^κ
pub fn temperature_from_potential_temperature(
  pressure : Double,
  potential_temperature : Double,
) -> Double {
  let ref_p : Double = pot_temp_ref_press
  potential_temperature * exner_function(pressure, reference_pressure=ref_p)
}

// ============================================================================
// 4. 混合比和比湿
// ============================================================================

///|
/// 计算混合比
/// 公式: r = ε * e / (p - e)
/// partial_press: 水汽分压 (Pa)
/// total_press: 总气压 (Pa)
/// molecular_weight_ratio: 分子量比，默认使用 ε (水汽/干空气)
pub fn mixing_ratio(
  partial_press : Double,
  total_press : Double,
  eps_param? : Double,
) -> Double {
  let eps = match eps_param {
    Some(e) => e
    None => molecular_weight_ratio
  }
  eps * partial_press / (total_press - partial_press)
}

///|
/// 计算饱和混合比
/// 公式: r_s = ε * e_s / (p - e_s)
pub fn saturation_mixing_ratio(
  total_press : Double,
  temperature : Double,
  phase : String,
) -> Double {
  let es = saturation_vapor_pressure(temperature, phase)
  let eps_val : Double = molecular_weight_ratio
  mixing_ratio(es, total_press, eps_param=eps_val)
}

///|
/// 从比湿计算混合比
/// 公式: w = q / (1 - q)
pub fn mixing_ratio_from_specific_humidity(
  specific_humidity : Double,
) -> Double {
  specific_humidity / (1.0 - specific_humidity)
}

///|
/// 从混合比计算比湿
/// 公式: q = w / (1 + w)
pub fn specific_humidity_from_mixing_ratio(mixing_ratio : Double) -> Double {
  mixing_ratio / (1.0 + mixing_ratio)
}

///|
/// 从露点和气压计算比湿
/// pressure: 气压 (Pa)
/// dewpoint: 露点温度 (K)
/// phase: "liquid" (默认), "solid", 或 "auto"
pub fn specific_humidity_from_dewpoint(
  pressure : Double,
  dewpoint : Double,
  phase : String,
) -> Double {
  let mixing_ratio_val = saturation_mixing_ratio(pressure, dewpoint, phase)
  specific_humidity_from_mixing_ratio(mixing_ratio_val)
}

///|
/// 计算水汽压
/// 公式: e = p * r / (r + ε)
pub fn vapor_pressure(pressure : Double, mixing_ratio_val : Double) -> Double {
  pressure * mixing_ratio_val / (molecular_weight_ratio + mixing_ratio_val)
}

// ============================================================================
// 5. 露点温度
// ============================================================================

///|
/// 从水汽压计算露点温度
/// 使用 Bolton (1980) 公式的逆函数
/// vapor_pressure: 水汽压 (Pa)
/// 返回: 露点温度 (K)
pub fn dewpoint(vapor_pressure : Double) -> Double {
  // Bolton 公式返回摄氏度，需要加上 273.15 转换为开尔文
  // 公式: T(°C) = 243.5 * ln(e / 6.112) / (17.67 - ln(e / 6.112))
  // 其中 e 是水汽压（毫巴），6.112 是 0°C 的饱和水汽压（毫巴）
  let e_mbar = vapor_pressure / 100.0 // 转换为毫巴
  let e0_mbar = sat_pressure_0c / 100.0 // 0°C 的饱和水汽压（毫巴）
  let val = @math.ln(e_mbar / e0_mbar)
  let td_celsius = 243.5 * val / (17.67 - val)
  273.15 + td_celsius // 转换为开尔文
}

///|
/// 从相对湿度计算露点温度
pub fn dewpoint_from_relative_humidity(
  temperature : Double,
  relative_humidity : Double,
  _phase : String,
) -> Double {
  let es = saturation_vapor_pressure(temperature, "liquid")
  let e = es * relative_humidity
  dewpoint(e)
}

///|
/// 从比湿和气压计算露点温度
pub fn dewpoint_from_specific_humidity(
  pressure : Double,
  specific_humidity : Double,
) -> Double {
  let mixing_ratio_val = mixing_ratio_from_specific_humidity(specific_humidity)
  let e = vapor_pressure(pressure, mixing_ratio_val)
  dewpoint(e)
}

// ============================================================================
// 6. 相对湿度
// ============================================================================

///|
/// 从混合比计算相对湿度
pub fn relative_humidity_from_mixing_ratio(
  pressure : Double,
  temperature : Double,
  mixing_ratio_val : Double,
  phase : String,
) -> Double {
  let e = vapor_pressure(pressure, mixing_ratio_val)
  let es = saturation_vapor_pressure(temperature, phase)
  e / es
}

///|
/// 从比湿计算相对湿度
pub fn relative_humidity_from_specific_humidity(
  pressure : Double,
  temperature : Double,
  specific_humidity : Double,
  phase : String,
) -> Double {
  let mixing_ratio_val = mixing_ratio_from_specific_humidity(specific_humidity)
  relative_humidity_from_mixing_ratio(
    pressure, temperature, mixing_ratio_val, phase,
  )
}

///|
/// 从露点计算相对湿度
pub fn relative_humidity_from_dewpoint(
  temperature : Double,
  dewpoint : Double,
  phase : String,
) -> Double {
  let e = saturation_vapor_pressure(dewpoint, phase)
  let es = saturation_vapor_pressure(temperature, phase)
  e / es
}

///|
/// 从相对湿度计算混合比
pub fn mixing_ratio_from_relative_humidity(
  pressure : Double,
  temperature : Double,
  relative_humidity : Double,
  phase : String,
) -> Double {
  let es = saturation_vapor_pressure(temperature, phase)
  let e = es * relative_humidity
  let eps_val : Double = molecular_weight_ratio
  mixing_ratio(e, pressure, eps_param=eps_val)
}

// ============================================================================
// 7. 虚温和密度
// ============================================================================

///|
/// 计算虚温
/// 公式: T_v = T * (1 + r/ε) / (1 + r)
/// temperature: 温度 (K)
/// mixing_ratio_val: 混合比
pub fn virtual_temperature(
  temperature : Double,
  mixing_ratio_val : Double,
) -> Double {
  temperature *
  (1.0 + mixing_ratio_val / molecular_weight_ratio) /
  (1.0 + mixing_ratio_val)
}

///|
/// 从露点计算虚温
pub fn virtual_temperature_from_dewpoint(
  pressure : Double,
  temperature : Double,
  dewpoint : Double,
) -> Double {
  let mixing_ratio_val = saturation_mixing_ratio(pressure, dewpoint, "liquid")
  virtual_temperature(temperature, mixing_ratio_val)
}

///|
/// 计算虚位温
pub fn virtual_potential_temperature(
  pressure : Double,
  temperature : Double,
  mixing_ratio_val : Double,
) -> Double {
  let tv = virtual_temperature(temperature, mixing_ratio_val)
  potential_temperature(pressure, tv)
}

///|
/// 计算空气密度
/// 公式: ρ = p / (R_d * T_v)
/// pressure: 气压 (Pa)
/// temperature: 温度 (K)
/// mixing_ratio_val: 混合比，默认 0（干空气）
pub fn density(
  pressure : Double,
  temperature : Double,
  mixing_ratio_val? : Double,
) -> Double {
  let mr = match mixing_ratio_val {
    Some(m) => m
    None => 0.0
  }
  let tv = virtual_temperature(temperature, mr)
  pressure / (dry_air_gas_constant * tv)
}

// ============================================================================
// 8. 相当位温和湿球温度
// ============================================================================

///|
/// 计算相当位温
/// 使用 Bolton (1980) 简化公式
/// pressure: 气压 (Pa)
/// temperature: 温度 (K)
/// dewpoint: 露点温度 (K)
pub fn equivalent_potential_temperature(
  pressure : Double,
  temperature : Double,
  dewpoint : Double,
) -> Double {
  let theta = potential_temperature(pressure, temperature)
  let mixing_ratio_val = saturation_mixing_ratio(pressure, dewpoint, "liquid")

  // Bolton (1980) 公式
  let tl = temperature - 55.0 // 简化处理
  let theta_e = theta *
    @math.exp(
      (3036.0 / tl - 1.78) * mixing_ratio_val * (1.0 + 0.448 * mixing_ratio_val),
    )
  theta_e
}

///|
/// 计算饱和相当位温
pub fn saturation_equivalent_potential_temperature(
  pressure : Double,
  temperature : Double,
) -> Double {
  equivalent_potential_temperature(pressure, temperature, temperature)
}

///|
/// 计算湿球温度（简化迭代方法）
/// pressure_pa: 气压 (Pa)
/// temperature_k: 温度 (K)
/// dewpoint_k: 露点温度 (K)
pub fn wet_bulb_temperature(
  pressure_pa : Double,
  temperature_k : Double,
  dewpoint_k : Double,
) -> Double {
  // 边界检查：如果温度低于露点，返回露点
  if temperature_k <= dewpoint_k {
    return dewpoint_k
  }

  // 初始猜测：湿球温度介于干球温度和露点温度之间
  let mut tw = (temperature_k + dewpoint_k) / 2.0

  // 确保初始值在合理范围内
  if tw < dewpoint_k {
    tw = dewpoint_k + 0.1
  }
  if tw > temperature_k {
    tw = temperature_k - 0.1
  }

  // 计算当前比湿（在循环外计算一次即可）
  let q = specific_humidity_from_dewpoint(pressure_pa, dewpoint_k, "liquid")

  // 迭代求解（最多 100 次）
  for _iter = 0; _iter < 100; _iter = _iter + 1 {
    // 计算湿球温度对应的比湿（应该是饱和的）
    let qs_w = specific_humidity_from_dewpoint(pressure_pa, tw, "liquid")

    // 检查是否出现 NaN 或无效值
    if qs_w != qs_w || qs_w <= 0.0 {
      // 如果计算失败，返回保守估计
      return (temperature_k + dewpoint_k) / 2.0
    }

    // 使用简化的湿球温度公式（基于能量平衡）
    // Tw = T - (Lv / cp) * (q - qs(Tw)) / (1 + q)
    let lv = latent_heat_vaporization(tw) // J/kg
    let cp = moist_air_specific_heat_pressure(q) // J/(K*kg)

    // 检查 cp 是否为零或无效
    if cp <= 0.0 || cp != cp {
      return (temperature_k + dewpoint_k) / 2.0
    }

    // 计算新的湿球温度
    let q_diff = q - qs_w
    let denominator = 1.0 + q
    if denominator <= 0.0 {
      return (temperature_k + dewpoint_k) / 2.0
    }
    let new_tw = temperature_k - lv / cp * q_diff / denominator

    // 检查新值是否有效
    if new_tw != new_tw ||
      new_tw < dewpoint_k - 10.0 ||
      new_tw > temperature_k + 10.0 {
      // 如果新值无效，使用阻尼回退
      tw = (tw + dewpoint_k) / 2.0
      continue
    }

    // 限制在合理范围内
    let clamped_tw = if new_tw < dewpoint_k {
      dewpoint_k + 0.01
    } else if new_tw > temperature_k {
      temperature_k - 0.01
    } else {
      new_tw
    }

    // 检查收敛（使用较小的容差）
    if Double::abs(clamped_tw - tw) < 0.001 {
      return clamped_tw
    }

    // 使用阻尼迭代（混合新旧值）以避免振荡
    tw = 0.7 * tw + 0.3 * clamped_tw
  }

  // 如果未收敛，返回最后一次迭代的结果（确保在合理范围内）
  if tw < dewpoint_k {
    dewpoint_k + 0.1
  } else if tw > temperature_k {
    temperature_k - 0.1
  } else {
    tw
  }
}

// ============================================================================
// 9. 基础气象函数（风速、风向等）
// ============================================================================

///|
/// 计算风速
/// u: 东西方向风速分量 (m/s)
/// v: 南北方向风速分量 (m/s)
pub fn wind_speed(u : Double, v : Double) -> Double {
  Double::sqrt(u * u + v * v)
}

///|
/// 计算风向
/// u: 东西方向风速分量 (m/s)
/// v: 南北方向风速分量 (m/s)
/// convention: "from" (默认，气象学约定) 或 "to" (海洋学约定)
/// 返回: 风向角度 (度)，0-360，0度为北
pub fn wind_direction(u : Double, v : Double, convention? : String) -> Double {
  let conv = match convention {
    Some(c) => c
    None => "from"
  }

  // 计算反正切，得到 -180 到 180 度的角度
  let angle_rad = @math.atan2(-v, -u)
  let angle_deg = angle_rad * 180.0 / 3.141592653589793

  // 转换为 0-360 度范围
  let mut wdir = 90.0 - angle_deg

  // 处理海洋学约定
  if conv == "to" {
    wdir = wdir - 180.0
  }

  // 规范化到 0-360 度
  if wdir <= 0.0 {
    wdir = wdir + 360.0
  }

  // 处理静风情况
  if Double::abs(u) < 0.001 && Double::abs(v) < 0.001 {
    0.0
  } else {
    wdir
  }
}

///|
/// 从风速和风向计算 u, v 分量
/// speed: 风速 (m/s)
/// wind_direction: 风向 (度)，0度为北，顺时针增加
pub fn wind_components(
  speed : Double,
  wind_direction_deg : Double,
) -> (Double, Double) {
  let angle_rad = (270.0 - wind_direction_deg) * 3.141592653589793 / 180.0
  let u = speed * @math.cos(angle_rad)
  let v = speed * @math.sin(angle_rad)
  (u, v)
}

///|
/// 计算科里奥利参数
/// latitude: 纬度 (度)
/// 返回: 科里奥利参数 (1/s)
pub fn coriolis_parameter(latitude : Double) -> Double {
  2.0 * earth_avg_angular_vel * @math.sin(latitude * 3.141592653589793 / 180.0)
}

// ============================================================================
// 10. 稳定度分析函数
// ============================================================================

///|
/// 计算 Brunt-Väisälä 频率的平方
/// height: 高度数组 (m)
/// potential_temperature: 位温数组 (K)
/// 返回: N² 数组 (1/s²)
pub fn brunt_vaisala_frequency_squared(
  height : Array[Double],
  potential_temperature : Array[Double],
) -> Array[Double]? {
  if height.length() < 2 || potential_temperature.length() != height.length() {
    None
  } else {
    let n = height.length()
    let g = earth_gravity
    let n2 = Array::make(n, 0.0)

    // 第一个点：使用前向差分
    if n > 1 {
      let dz = height[1] - height[0]
      let dtheta = potential_temperature[1] - potential_temperature[0]
      let n2_val = if dz > 0.0 && potential_temperature[0] > 0.0 {
        g / potential_temperature[0] * (dtheta / dz)
      } else {
        0.0
      }
      n2[0] = n2_val
    }

    // 中间点：使用中心差分
    for i = 1; i < n - 1; i = i + 1 {
      let dz = height[i + 1] - height[i - 1]
      let dtheta = potential_temperature[i + 1] - potential_temperature[i - 1]
      let n2_val = if dz > 0.0 && potential_temperature[i] > 0.0 {
        g / potential_temperature[i] * (dtheta / dz)
      } else {
        0.0
      }
      n2[i] = n2_val
    }

    // 最后一个点：使用后向差分
    if n > 1 {
      let dz = height[n - 1] - height[n - 2]
      let dtheta = potential_temperature[n - 1] - potential_temperature[n - 2]
      let n2_val = if dz > 0.0 && potential_temperature[n - 1] > 0.0 {
        g / potential_temperature[n - 1] * (dtheta / dz)
      } else {
        0.0
      }
      n2[n - 1] = n2_val
    }
    Some(n2)
  }
}

///|
/// 计算 Brunt-Väisälä 频率
pub fn brunt_vaisala_frequency(
  height : Array[Double],
  potential_temperature : Array[Double],
) -> Array[Double]? {
  match brunt_vaisala_frequency_squared(height, potential_temperature) {
    Some(n2) => {
      let n = Array::make(n2.length(), 0.0)
      for i = 0; i < n2.length(); i = i + 1 {
        n[i] = if n2[i] > 0.0 { Double::sqrt(n2[i]) } else { 0.0 }
      }
      Some(n)
    }
    None => None
  }
}

///|
/// 计算 Brunt-Väisälä 周期
pub fn brunt_vaisala_period(
  height : Array[Double],
  potential_temperature : Array[Double],
) -> Array[Double]? {
  match brunt_vaisala_frequency(height, potential_temperature) {
    Some(n) => {
      let period = Array::make(n.length(), 0.0)
      for i = 0; i < n.length(); i = i + 1 {
        period[i] = if n[i] > 0.0 {
          2.0 * 3.141592653589793 / n[i]
        } else {
          0.0
        }
      }
      Some(period)
    }
    None => None
  }
}

///|
/// 计算静力稳定度 (Static Stability)
/// 公式: σ = - (1/θ) * (∂θ/∂p)
/// 返回静力稳定度 (1/Pa)
/// 正值表示稳定，负值表示不稳定
pub fn static_stability(
  pressure : Array[Double],
  temperature : Array[Double],
) -> Array[Double]? {
  if pressure.length() != temperature.length() || pressure.length() < 2 {
    None
  } else {
    let n = pressure.length()
    let stability = Array::make(n, 0.0)

    // 计算位温
    let theta = Array::make(n, 0.0)
    for i = 0; i < n; i = i + 1 {
      theta[i] = potential_temperature(pressure[i], temperature[i])
    }

    // 第一个点：使用前向差分
    if n > 1 {
      let dp = pressure[0] - pressure[1] // 注意：压力向上递减
      let dtheta = theta[1] - theta[0]
      stability[0] = if dp != 0.0 && theta[0] > 0.0 {
        -1.0 / theta[0] * (dtheta / dp)
      } else {
        0.0
      }
    }

    // 中间点：使用中心差分
    for i = 1; i < n - 1; i = i + 1 {
      let dp = pressure[i - 1] - pressure[i + 1]
      let dtheta = theta[i + 1] - theta[i - 1]
      stability[i] = if dp != 0.0 && theta[i] > 0.0 {
        -1.0 / theta[i] * (dtheta / dp)
      } else {
        0.0
      }
    }

    // 最后一个点：使用后向差分
    if n > 1 {
      let dp = pressure[n - 2] - pressure[n - 1]
      let dtheta = theta[n - 1] - theta[n - 2]
      stability[n - 1] = if dp != 0.0 && theta[n - 1] > 0.0 {
        -1.0 / theta[n - 1] * (dtheta / dp)
      } else {
        0.0
      }
    }
    Some(stability)
  }
}

// ============================================================================
// 11. 厚度计算
// ============================================================================

///|
/// 计算静力学厚度
/// pressure_bottom: 底层气压 (Pa)
/// pressure_top: 顶层气压 (Pa)
/// temperature: 平均温度 (K)
/// mixing_ratio_val: 混合比（可选）
pub fn thickness_hydrostatic(
  pressure_bottom : Double,
  pressure_top : Double,
  temperature : Double,
  mixing_ratio_val? : Double,
) -> Double {
  let mr = match mixing_ratio_val {
    Some(m) => m
    None => 0.0
  }
  let tv = virtual_temperature(temperature, mr)
  let r = dry_air_gas_constant
  let g = earth_gravity

  // 使用虚温的静力学方程
  // h = (R_d * T_v / g) * ln(p_bottom / p_top)
  r * tv / g * @math.ln(pressure_bottom / pressure_top)
}

// ============================================================================
// 12. 静态能量
// ============================================================================

///|
/// 计算干静力能
/// height: 高度 (m)
/// temperature: 温度 (K)
pub fn dry_static_energy(height : Double, temperature : Double) -> Double {
  let cp = dry_air_spec_heat_press
  let g = earth_gravity
  cp * temperature + g * height
}

///|
/// 计算湿静力能
/// height: 高度 (m)
/// temperature: 温度 (K)
/// specific_humidity: 比湿 (kg/kg)
pub fn moist_static_energy(
  height : Double,
  temperature : Double,
  specific_humidity : Double,
) -> Double {
  let dse = dry_static_energy(height, temperature)
  let lv = latent_heat_vaporization(temperature)
  dse + lv * specific_humidity
}

// ============================================================================
// 13. 湿绝热过程和对流分析
// ============================================================================

///|
/// 计算湿绝热过程中的温度
/// 使用迭代方法计算气块沿湿绝热线移动时的温度
/// pressure: 目标气压 (Pa)
/// temperature: 起始温度 (K)
/// reference_pressure: 参考气压 (Pa)，默认使用起始气压
/// max_iters: 最大迭代次数，默认 50
/// eps: 收敛容差，默认 1e-6
/// 返回: 目标气压下的温度 (K)
pub fn moist_lapse(
  pressure : Double,
  temperature : Double,
  reference_pressure? : Double,
  max_iters? : Int,
  eps? : Double,
) -> Double {
  let ref_p = match reference_pressure {
    Some(p) => p
    None => pressure
  }
  let _max_iter = match max_iters {
    Some(n) => n
    None => 50
  }
  let tolerance = match eps {
    Some(e) => e
    None => 0.000001
  }

  // 如果压力相同，直接返回温度
  if (pressure - ref_p).abs() < tolerance {
    temperature
  } else {
    // 计算初始混合比（假设饱和）
    let initial_mr = saturation_mixing_ratio(ref_p, temperature, "liquid")

    // 使用迭代方法计算湿绝热温度
    // 基于能量守恒：C_p * dT = R * T * dp/p - L * dr
    let mut current_temp = temperature
    let mut current_p = ref_p
    let mut current_mr = initial_mr

    // 使用数值积分方法，将压力变化分成小步
    let dp = pressure - ref_p
    let n_steps = 20 // 分成20步以提高精度
    let dp_step = dp / Double::from_int(n_steps)
    for _step = 0; _step < n_steps; _step = _step + 1 {
      current_p = current_p + dp_step

      // 计算当前温度下的饱和混合比
      let mr_sat = saturation_mixing_ratio(current_p, current_temp, "liquid")

      // 如果当前混合比超过饱和，则按湿绝热过程
      if current_mr >= mr_sat {
        // 湿绝热过程：潜热释放
        let lv = latent_heat_vaporization(current_temp)
        let cp = moist_air_specific_heat_pressure(
          specific_humidity_from_mixing_ratio(current_mr),
        )
        let r = moist_air_gas_constant(
          specific_humidity_from_mixing_ratio(current_mr),
        )

        // 计算温度变化
        // dT = (R*T/p - L*dr_sat/dT) * dp / (C_p + L*dr_sat/dT)
        // 简化：使用近似公式
        let gamma_m = earth_gravity /
          cp *
          (1.0 + lv * mr_sat / (r * current_temp)) /
          (1.0 + lv * lv * mr_sat / (cp * r * current_temp * current_temp))
        let dt = -gamma_m * dp_step / (earth_gravity * current_temp / r)
        current_temp = current_temp + dt

        // 更新混合比为饱和值
        current_mr = mr_sat
      } else {
        // 干绝热过程
        let cp = dry_air_spec_heat_press
        let r = dry_air_gas_constant
        let gamma_d = earth_gravity / cp
        let dt = -gamma_d * dp_step / (earth_gravity * current_temp / r)
        current_temp = current_temp + dt
      }

      // 防止温度过低或过高
      if current_temp < 150.0 {
        current_temp = 150.0
      }
      if current_temp > 400.0 {
        current_temp = 400.0
      }
    }
    current_temp
  }
}

///|
/// 计算抬升凝结高度 (LCL - Lifting Condensation Level)
/// 气块从起始位置干绝热抬升至饱和的高度
/// pressure: 起始气压 (Pa)
/// temperature: 起始温度 (K)
/// dewpoint: 起始露点温度 (K)
/// max_iters: 最大迭代次数，默认 50
/// eps: 收敛容差，默认 1e-6
/// 返回: (LCL 气压 (Pa), LCL 温度 (K))
pub fn lcl(
  pressure : Double,
  temperature : Double,
  dewpoint : Double,
  max_iters? : Int,
  eps? : Double,
) -> (Double, Double) {
  let max_iter = match max_iters {
    Some(n) => n
    None => 50
  }
  let tolerance = match eps {
    Some(e) => e
    None => 0.000001
  }

  // 如果温度已经低于等于露点，说明已经饱和
  if temperature <= dewpoint {
    (pressure, temperature)
  } else {
    // 计算起始混合比（从露点）
    let initial_mr = saturation_mixing_ratio(pressure, dewpoint, "liquid")

    // 使用迭代方法找到 LCL
    // LCL 是气块干绝热抬升至饱和时的位置
    // 使用 Bolton (1980) 公式的改进版本
    let mut current_p = pressure
    let mut current_t = temperature
    let mut iter = 0
    while iter < max_iter {
      // 计算当前温度下的饱和混合比
      let mr_sat = saturation_mixing_ratio(current_p, current_t, "liquid")

      // 检查是否达到饱和
      if (mr_sat - initial_mr).abs() < tolerance * initial_mr {
        break
      }

      // 如果还没饱和，继续干绝热抬升
      // 使用二分法或牛顿法找到 LCL 气压
      // 简化：使用 Bolton 公式的近似
      let t_ratio = current_t / dewpoint
      let p_lcl = pressure * Double::pow(t_ratio, -3.5)

      // 计算在 LCL 处的温度（干绝热过程）
      let t_lcl = dewpoint * Double::pow(p_lcl / pressure, poisson_exponent)

      // 更新
      let dp = (p_lcl - current_p).abs()
      if dp < tolerance * pressure {
        current_p = p_lcl
        current_t = t_lcl
        break
      }
      current_p = p_lcl
      current_t = t_lcl
      iter = iter + 1
    }
    (current_p, current_t)
  }
}

///|
/// 计算气块剖面（温度随压力的变化）
/// 从起始位置干绝热抬升到 LCL，然后湿绝热抬升
/// pressure: 起始气压 (Pa)
/// temperature: 起始温度 (K)
/// dewpoint: 起始露点温度 (K)
/// 返回: 温度剖面数组（对应压力数组）
pub fn parcel_profile(
  pressure : Array[Double],
  temperature : Double,
  dewpoint : Double,
) -> Array[Double]? {
  if pressure.length() == 0 {
    None
  } else {
    let start_pressure = pressure[0]
    let start_temp = temperature
    let start_dewpoint = dewpoint

    // 计算 LCL
    let (lcl_p, lcl_t) = lcl(start_pressure, start_temp, start_dewpoint)
    let profile = Array::make(pressure.length(), 0.0)
    for i = 0; i < pressure.length(); i = i + 1 {
      let p = pressure[i]
      if p >= lcl_p {
        // LCL 以下：干绝热过程
        let theta = potential_temperature(start_pressure, start_temp)
        profile[i] = temperature_from_potential_temperature(p, theta)
      } else {
        // LCL 以上：湿绝热过程
        profile[i] = moist_lapse(p, lcl_t, reference_pressure=lcl_p)
      }
    }
    Some(profile)
  }
}

///|
/// 计算自由对流高度 (LFC - Level of Free Convection)
/// 气块温度首次高于环境温度的高度
/// pressure: 气压数组 (Pa)，从地面向上递减
/// temperature: 环境温度数组 (K)
/// parcel_profile: 气块温度剖面 (K)
/// 返回: (LFC 气压 (Pa), LFC 温度 (K))，如果不存在则返回 None
pub fn lfc(
  pressure : Array[Double],
  temperature : Array[Double],
  parcel_profile : Array[Double],
) -> (Double, Double)? {
  if pressure.length() != temperature.length() ||
    pressure.length() != parcel_profile.length() ||
    pressure.length() < 2 {
    None
  } else {
    // 找到气块温度首次高于环境温度的位置
    let mut found = false
    let mut lfc_p = 0.0
    let mut lfc_t = 0.0

    // 从下往上查找（数组索引从大到小，因为压力向上递减）
    for i = 0; i < pressure.length() - 1; i = i + 1 {
      let p_parcel = parcel_profile[i]
      let t_env = temperature[i]

      // 检查气块温度是否高于环境温度
      if p_parcel > t_env && !found {
        // 找到 LFC
        found = true
        lfc_p = pressure[i]
        lfc_t = p_parcel
        break
      }
    }
    if found {
      Some((lfc_p, lfc_t))
    } else {
      None
    }
  }
}

///|
/// 计算平衡高度 (EL - Equilibrium Level)
/// 气块温度再次等于环境温度的高度（LFC 之上）
/// pressure: 气压数组 (Pa)，从地面向上递减
/// temperature: 环境温度数组 (K)
/// parcel_profile: 气块温度剖面 (K)
/// which: "top" 或 "bottom"，返回哪个 EL（默认 "top"）
/// 返回: (EL 气压 (Pa), EL 温度 (K))，如果不存在则返回 None
pub fn el(
  pressure : Array[Double],
  temperature : Array[Double],
  parcel_profile : Array[Double],
  which? : String,
) -> (Double, Double)? {
  let which_el = match which {
    Some(w) => w
    None => "top"
  }
  if pressure.length() != temperature.length() ||
    pressure.length() != parcel_profile.length() ||
    pressure.length() < 2 {
    None
  } else {
    // 先找到 LFC
    let lfc_result = lfc(pressure, temperature, parcel_profile)
    match lfc_result {
      Some((lfc_p, _lfc_t)) => {
        // 找到 LFC 后，在 LFC 之上查找 EL
        let mut found = false
        let mut el_p = 0.0
        let mut el_t = 0.0
        let mut lfc_index = 0

        // 找到 LFC 的索引
        for i = 0; i < pressure.length(); i = i + 1 {
          if (pressure[i] - lfc_p).abs() < 1.0 {
            lfc_index = i
            break
          }
        }

        // 从 LFC 向上查找 EL
        if which_el == "top" {
          // 找最高的 EL
          for i = lfc_index + 1; i < pressure.length(); i = i + 1 {
            let p_parcel = parcel_profile[i]
            let t_env = temperature[i]

            // 如果气块温度等于或低于环境温度，找到 EL
            if p_parcel <= t_env {
              found = true
              el_p = pressure[i]
              el_t = p_parcel
              // 继续找最高的
            }
          }
        } else {
          // 找最低的 EL（第一个）
          for i = lfc_index + 1; i < pressure.length(); i = i + 1 {
            let p_parcel = parcel_profile[i]
            let t_env = temperature[i]
            if p_parcel <= t_env {
              found = true
              el_p = pressure[i]
              el_t = p_parcel
              break
            }
          }
        }
        if found {
          Some((el_p, el_t))
        } else {
          None
        }
      }
      None => None
    }
  }
}

///|
/// 计算对流可用位能 (CAPE) 和对流抑制能量 (CIN)
/// pressure: 气压数组 (Pa)，从地面向上递减
/// temperature: 环境温度数组 (K)
/// dewpoint: 环境露点数组 (K)
/// parcel_profile: 气块温度剖面 (K)
/// 返回: (CAPE (J/kg), CIN (J/kg))
pub fn cape_cin(
  pressure : Array[Double],
  temperature : Array[Double],
  dewpoint : Array[Double],
  parcel_profile : Array[Double],
) -> (Double, Double) {
  if pressure.length() != temperature.length() ||
    pressure.length() != dewpoint.length() ||
    pressure.length() != parcel_profile.length() ||
    pressure.length() < 2 {
    (0.0, 0.0)
  } else {
    let r = dry_air_gas_constant
    let mut cape = 0.0
    let mut cin = 0.0

    // 找到 LFC
    let lfc_result = lfc(pressure, temperature, parcel_profile)
    let mut lfc_index = 0
    let mut el_index = pressure.length() - 1
    match lfc_result {
      Some((lfc_p, _lfc_t)) => {
        // 找到 LFC 的索引
        for i = 0; i < pressure.length(); i = i + 1 {
          if (pressure[i] - lfc_p).abs() < 1.0 {
            lfc_index = i
            break
          }
        }

        // 找到 EL
        let el_result = el(pressure, temperature, parcel_profile, which="top")
        match el_result {
          Some((el_p, _el_t)) =>
            for i = 0; i < pressure.length(); i = i + 1 {
              if (pressure[i] - el_p).abs() < 1.0 {
                el_index = i
                break
              }
            }
          None => ignore(0)
        }
      }
      None =>
        // 没有 LFC，返回 0
        return (0.0, 0.0)
    }

    // 计算 CIN（从地面到 LFC）
    for i = 0; i < lfc_index; i = i + 1 {
      let t_parcel = parcel_profile[i]
      let t_env = temperature[i]
      let p = pressure[i]
      if t_parcel < t_env {
        // 负浮力区域，计算 CIN
        let tv_parcel = virtual_temperature(t_parcel, 0.0)
        let tv_env = virtual_temperature(t_env, 0.0)
        let buoyancy = r * (tv_parcel - tv_env) / p
        cin = cin + buoyancy * (pressure[i] - pressure[i + 1])
      }
    }

    // 计算 CAPE（从 LFC 到 EL）
    for i = lfc_index; i < el_index; i = i + 1 {
      let t_parcel = parcel_profile[i]
      let t_env = temperature[i]
      let p = pressure[i]
      if t_parcel > t_env {
        // 正浮力区域，计算 CAPE
        let tv_parcel = virtual_temperature(t_parcel, 0.0)
        let tv_env = virtual_temperature(t_env, 0.0)
        let buoyancy = r * (tv_parcel - tv_env) / p
        cape = cape + buoyancy * (pressure[i] - pressure[i + 1])
      }
    }
    (cape, cin)
  }
}

// ============================================================================
// 14. 运动学函数
// ============================================================================

///|
/// 计算垂直涡度
/// 公式: ζ = ∂v/∂x - ∂u/∂y
/// u: 东西方向风速 (m/s)，二维数组
/// v: 南北方向风速 (m/s)，二维数组
/// dx: x 方向网格间距 (m)
/// dy: y 方向网格间距 (m)
/// 返回: 垂直涡度数组 (1/s)
pub fn vorticity(
  u : Array[Array[Double]],
  v : Array[Array[Double]],
  dx : Double,
  dy : Double,
) -> Array[Array[Double]]? {
  if u.length() != v.length() || u.length() == 0 {
    None
  } else {
    let ny = u.length()
    let nx = if u[0].length() > 0 { u[0].length() } else { 0 }
    if nx == 0 || v.length() == 0 || v[0].length() != nx {
      None
    } else {
      let vort = Array::make(ny, Array::make(nx, 0.0))

      // 使用中心差分计算涡度
      for j = 0; j < ny; j = j + 1 {
        for i = 0; i < nx; i = i + 1 {
          // ∂v/∂x
          let dv_dx = if i > 0 && i < nx - 1 {
            (v[j][i + 1] - v[j][i - 1]) / (2.0 * dx)
          } else if i == 0 {
            (v[j][i + 1] - v[j][i]) / dx
          } else {
            (v[j][i] - v[j][i - 1]) / dx
          }

          // ∂u/∂y
          let du_dy = if j > 0 && j < ny - 1 {
            (u[j + 1][i] - u[j - 1][i]) / (2.0 * dy)
          } else if j == 0 {
            (u[j + 1][i] - u[j][i]) / dy
          } else {
            (u[j][i] - u[j - 1][i]) / dy
          }

          // 涡度 = ∂v/∂x - ∂u/∂y
          vort[j][i] = dv_dx - du_dy
        }
      }
      Some(vort)
    }
  }
}

///|
/// 计算水平散度
/// 公式: D = ∂u/∂x + ∂v/∂y
/// u: 东西方向风速 (m/s)，二维数组
/// v: 南北方向风速 (m/s)，二维数组
/// dx: x 方向网格间距 (m)
/// dy: y 方向网格间距 (m)
/// 返回: 散度数组 (1/s)
pub fn divergence(
  u : Array[Array[Double]],
  v : Array[Array[Double]],
  dx : Double,
  dy : Double,
) -> Array[Array[Double]]? {
  if u.length() != v.length() || u.length() == 0 {
    None
  } else {
    let ny = u.length()
    let nx = if u[0].length() > 0 { u[0].length() } else { 0 }
    if nx == 0 || v.length() == 0 || v[0].length() != nx {
      None
    } else {
      let div = Array::make(ny, Array::make(nx, 0.0))

      // 使用中心差分计算散度
      for j = 0; j < ny; j = j + 1 {
        for i = 0; i < nx; i = i + 1 {
          // ∂u/∂x
          let du_dx = if i > 0 && i < nx - 1 {
            (u[j][i + 1] - u[j][i - 1]) / (2.0 * dx)
          } else if i == 0 {
            (u[j][i + 1] - u[j][i]) / dx
          } else {
            (u[j][i] - u[j][i - 1]) / dx
          }

          // ∂v/∂y
          let dv_dy = if j > 0 && j < ny - 1 {
            (v[j + 1][i] - v[j - 1][i]) / (2.0 * dy)
          } else if j == 0 {
            (v[j + 1][i] - v[j][i]) / dy
          } else {
            (v[j][i] - v[j - 1][i]) / dy
          }

          // 散度 = ∂u/∂x + ∂v/∂y
          div[j][i] = du_dx + dv_dy
        }
      }
      Some(div)
    }
  }
}

///|
/// 计算地转风
/// 从位势高度计算地转风速分量
/// height: 位势高度 (m)，二维数组
/// dx: x 方向网格间距 (m)
/// dy: y 方向网格间距 (m)
/// latitude_deg: 纬度 (度)
/// 返回: (u_g (m/s), v_g (m/s))，u_g 为东西分量，v_g 为南北分量
pub fn geostrophic_wind(
  height : Array[Array[Double]],
  dx : Double,
  dy : Double,
  latitude_deg : Double,
) -> (Array[Array[Double]], Array[Array[Double]])? {
  if height.length() == 0 {
    None
  } else {
    let ny = height.length()
    let nx = if height[0].length() > 0 { height[0].length() } else { 0 }
    if nx == 0 {
      None
    } else {
      // 计算科里奥利参数
      let f = coriolis_parameter(latitude_deg)
      let u_g = Array::make(ny, Array::make(nx, 0.0))
      let v_g = Array::make(ny, Array::make(nx, 0.0))

      // 地转风公式:
      // u_g = - (g/f) * (∂h/∂y)
      // v_g = (g/f) * (∂h/∂x)
      let g = earth_gravity
      let coeff = g / f
      for j = 0; j < ny; j = j + 1 {
        for i = 0; i < nx; i = i + 1 {
          // ∂h/∂x (用于计算 v_g)
          let dh_dx = if i > 0 && i < nx - 1 {
            (height[j][i + 1] - height[j][i - 1]) / (2.0 * dx)
          } else if i == 0 {
            (height[j][i + 1] - height[j][i]) / dx
          } else {
            (height[j][i] - height[j][i - 1]) / dx
          }

          // ∂h/∂y (用于计算 u_g)
          let dh_dy = if j > 0 && j < ny - 1 {
            (height[j + 1][i] - height[j - 1][i]) / (2.0 * dy)
          } else if j == 0 {
            (height[j + 1][i] - height[j][i]) / dy
          } else {
            (height[j][i] - height[j - 1][i]) / dy
          }

          // u_g = - (g/f) * (∂h/∂y)
          u_g[j][i] = -coeff * dh_dy

          // v_g = (g/f) * (∂h/∂x)
          v_g[j][i] = coeff * dh_dx
        }
      }
      Some((u_g, v_g))
    }
  }
}

///|
/// 计算位涡（正压）
/// 公式: PV = (ζ + f) / H
/// u: 东西方向风速 (m/s)，二维数组
/// v: 南北方向风速 (m/s)，二维数组
/// dx: x 方向网格间距 (m)
/// dy: y 方向网格间距 (m)
/// latitude_deg: 纬度 (度)
/// height: 位势高度或厚度 (m)，二维数组
/// 返回: 位涡数组 (PVU, 1e-6 m²·K·s⁻¹·kg⁻¹)
pub fn potential_vorticity_barotropic(
  u : Array[Array[Double]],
  v : Array[Array[Double]],
  dx : Double,
  dy : Double,
  latitude_deg : Double,
  height : Array[Array[Double]],
) -> Array[Array[Double]]? {
  // 计算相对涡度
  let vort_result = vorticity(u, v, dx, dy)
  match vort_result {
    Some(vort) =>
      if height.length() != vort.length() || height.length() == 0 {
        None
      } else {
        let ny = vort.length()
        let nx = if vort[0].length() > 0 { vort[0].length() } else { 0 }
        if nx == 0 || height[0].length() != nx {
          None
        } else {
          let f = coriolis_parameter(latitude_deg)
          let pv = Array::make(ny, Array::make(nx, 0.0))
          for j = 0; j < ny; j = j + 1 {
            for i = 0; i < nx; i = i + 1 {
              // 绝对涡度 = 相对涡度 + 科里奥利参数
              let abs_vort = vort[j][i] + f

              // 位涡 = 绝对涡度 / 高度
              if height[j][i] > 0.0 {
                pv[j][i] = abs_vort / height[j][i]
              } else {
                pv[j][i] = 0.0
              }
            }
          }
          Some(pv)
        }
      }
    None => None
  }
}

///|
/// 计算位涡（斜压）
/// 公式: PV = -g * (ζ + f) * (∂θ/∂p)
/// u: 东西方向风速 (m/s)，二维数组
/// v: 南北方向风速 (m/s)，二维数组
/// temperature: 温度 (K)，二维数组
/// pressure: 气压 (Pa)，一维数组（垂直方向）
/// dx: x 方向网格间距 (m)
/// dy: y 方向网格间距 (m)
/// latitude_deg: 纬度 (度)
/// 返回: 位涡数组 (PVU, 1e-6 m²·K·s⁻¹·kg⁻¹)
pub fn potential_vorticity_baroclinic(
  u : Array[Array[Double]],
  v : Array[Array[Double]],
  temperature : Array[Array[Double]],
  pressure : Array[Double],
  dx : Double,
  dy : Double,
  latitude_deg : Double,
) -> Array[Array[Double]]? {
  // 计算相对涡度
  let vort_result = vorticity(u, v, dx, dy)
  match vort_result {
    Some(vort) =>
      if temperature.length() != vort.length() || temperature.length() == 0 {
        None
      } else {
        let ny = vort.length()
        let nx = if vort[0].length() > 0 { vort[0].length() } else { 0 }
        if nx == 0 || temperature[0].length() != nx || pressure.length() < 2 {
          None
        } else {
          let f = coriolis_parameter(latitude_deg)
          let g = earth_gravity
          let pv = Array::make(ny, Array::make(nx, 0.0))

          // 计算位温梯度
          for j = 0; j < ny; j = j + 1 {
            for i = 0; i < nx; i = i + 1 {
              // 计算位温（使用平均气压）
              let avg_p = if pressure.length() > 0 {
                let mut sum = 0.0
                for k = 0; k < pressure.length(); k = k + 1 {
                  sum = sum + pressure[k]
                }
                sum / Double::from_int(pressure.length())
              } else {
                100000.0
              }
              let theta = potential_temperature(avg_p, temperature[j][i])

              // 计算 ∂θ/∂p（使用中心差分）
              let dtheta_dp = if pressure.length() > 1 {
                // 简化为使用压力差的倒数
                let dp = pressure[0] - pressure[pressure.length() - 1]
                if dp > 0.0 {
                  // 近似：使用平均温度梯度
                  -theta / dp
                } else {
                  0.0
                }
              } else {
                0.0
              }

              // 绝对涡度
              let abs_vort = vort[j][i] + f

              // 位涡 = -g * (ζ + f) * (∂θ/∂p)
              pv[j][i] = -g * abs_vort * dtheta_dp
            }
          }
          Some(pv)
        }
      }
    None => None
  }
}

// ============================================================================
// 11. 天气指数函数
// ============================================================================

///|
/// 在指定气压层进行线性插值
/// pressure: 气压数组 (Pa)，从高到低排列
/// values: 对应的值数组
/// target_pressure: 目标气压 (Pa)
/// 返回: 插值得到的值，如果无法插值则返回 None
fn interpolate_to_pressure(
  pressure : Array[Double],
  values : Array[Double],
  target_pressure : Double,
) -> Double? {
  if pressure.length() != values.length() || pressure.length() < 2 {
    None
  } else {
    // 找到包含目标气压的区间
    let mut found = false
    let mut idx = 0

    // 如果目标气压在数组范围内
    if target_pressure <= pressure[0] &&
      target_pressure >= pressure[pressure.length() - 1] {
      // 找到目标气压所在的区间
      for i = 0; i < pressure.length() - 1; i = i + 1 {
        if (
            target_pressure <= pressure[i] && target_pressure >= pressure[i + 1]
          ) ||
          (target_pressure >= pressure[i] && target_pressure <= pressure[i + 1]) {
          idx = i
          found = true
          break
        }
      }
    }
    if !found {
      // 如果目标气压在数组外，使用最近的点
      if target_pressure > pressure[0] {
        Some(values[0])
      } else if target_pressure < pressure[pressure.length() - 1] {
        Some(values[pressure.length() - 1])
      } else {
        None
      }
    } else {
      // 线性插值
      let p1 = pressure[idx]
      let p2 = pressure[idx + 1]
      let v1 = values[idx]
      let v2 = values[idx + 1]
      if Double::abs(p2 - p1) < 1.0 {
        // 如果压力差太小，直接返回平均值
        Some((v1 + v2) / 2.0)
      } else {
        let ratio = (target_pressure - p1) / (p2 - p1)
        Some(v1 + ratio * (v2 - v1))
      }
    }
  }
}

///|
/// 计算抬升指数 (Lifted Index)
/// LI = T500 - Tp500
/// 其中 T500 是 500 hPa 的观测温度，Tp500 是气块在 500 hPa 的温度
/// pressure: 气压数组 (Pa)，从高到低排列
/// temperature: 环境温度数组 (K)
/// parcel_profile: 气块温度剖面 (K)
/// 返回: 抬升指数 (K)，负值表示不稳定
pub fn lifted_index(
  pressure : Array[Double],
  temperature : Array[Double],
  parcel_profile : Array[Double],
) -> Double? {
  if pressure.length() != temperature.length() ||
    pressure.length() != parcel_profile.length() {
    None
  } else {
    // 500 hPa = 50000 Pa
    let p500 = 50000.0

    // 插值得到 500 hPa 的观测温度
    let t500_result = interpolate_to_pressure(pressure, temperature, p500)

    // 插值得到 500 hPa 的气块温度
    let tp500_result = interpolate_to_pressure(pressure, parcel_profile, p500)
    match (t500_result, tp500_result) {
      (Some(t500), Some(tp500)) => Some(t500 - tp500)
      _ => None
    }
  }
}

///|
/// 计算 K 指数
/// K = (T850 - T500) + Td850 - (T700 - Td700)
/// pressure: 气压数组 (Pa)，从高到低排列
/// temperature: 温度数组 (K)
/// dewpoint: 露点温度数组 (K)
/// 返回: K 指数 (K)，正值越大表示对流潜力越大
pub fn k_index(
  pressure : Array[Double],
  temperature : Array[Double],
  dewpoint : Array[Double],
) -> Double? {
  if pressure.length() != temperature.length() ||
    pressure.length() != dewpoint.length() {
    None
  } else {
    // 标准气压层：850, 700, 500 hPa
    let p850 = 85000.0
    let p700 = 70000.0
    let p500 = 50000.0

    // 插值得到各层的温度和露点
    let t850_result = interpolate_to_pressure(pressure, temperature, p850)
    let t700_result = interpolate_to_pressure(pressure, temperature, p700)
    let t500_result = interpolate_to_pressure(pressure, temperature, p500)
    let td850_result = interpolate_to_pressure(pressure, dewpoint, p850)
    let td700_result = interpolate_to_pressure(pressure, dewpoint, p700)
    match (t850_result, t700_result, t500_result, td850_result, td700_result) {
      (Some(t850), Some(t700), Some(t500), Some(td850), Some(td700)) =>
        // K = (T850 - T500) + Td850 - (T700 - Td700)
        Some(t850 - t500 + td850 - (t700 - td700))
      _ => None
    }
  }
}

///|
/// 计算可降水量 (Precipitable Water)
/// 通过对比湿进行气压积分得到
/// pressure: 气压数组 (Pa)，从高到低排列
/// dewpoint: 露点温度数组 (K)
/// bottom: 底层气压 (Pa)，可选，默认使用数组最小值
/// top: 顶层气压 (Pa)，可选，默认使用数组最大值
/// 返回: 可降水量 (mm)
pub fn precipitable_water(
  pressure : Array[Double],
  dewpoint : Array[Double],
  bottom? : Double,
  top? : Double,
) -> Double? {
  if pressure.length() != dewpoint.length() || pressure.length() < 2 {
    None
  } else {
    let p_bottom = match bottom {
      Some(b) => b
      None => pressure[pressure.length() - 1] // 最低气压（数组末尾）
    }
    let p_top = match top {
      Some(t) => t
      None => pressure[0] // 最高气压（数组开头）
    }

    // 找到积分区间
    let mut start_idx = 0
    let mut end_idx = pressure.length() - 1

    // 找到顶部和底部对应的索引
    for i = 0; i < pressure.length(); i = i + 1 {
      if pressure[i] <= p_top && start_idx == 0 {
        start_idx = i
      }
      if pressure[i] >= p_bottom {
        end_idx = i
      }
    }
    if start_idx >= end_idx {
      None
    } else {
      // 计算可降水量
      // PW = (1/g) * ∫ q dp
      // 其中 q 是比湿，g 是重力加速度
      let g = earth_gravity
      let mut pw = 0.0

      // 使用梯形法则进行积分
      for i = start_idx; i < end_idx; i = i + 1 {
        let p1 = pressure[i]
        let p2 = pressure[i + 1]
        let td1 = dewpoint[i]
        let td2 = dewpoint[i + 1]

        // 计算比湿
        let q1 = specific_humidity_from_dewpoint(p1, td1, "liquid")
        let q2 = specific_humidity_from_dewpoint(p2, td2, "liquid")

        // 梯形积分
        let dp = Double::abs(p2 - p1)
        let q_avg = (q1 + q2) / 2.0

        // 积分：∫ q dp / g
        // 转换为 mm: (kg/kg * Pa) / (m/s²) * (1000 mm/m) = mm
        pw = pw + q_avg * dp / g * 1000.0
      }
      Some(pw)
    }
  }
}

// ============================================================================
// 12. 运动学扩展函数
// ============================================================================

///|
/// 计算风暴相对螺旋度 (Storm Relative Helicity)
/// 公式: H = ∫ (v - vs) · (∇ × v) dz
/// height: 高度数组 (m)，从低到高排列
/// u: 东西方向风速分量数组 (m/s)
/// v: 南北方向风速分量数组 (m/s)
/// depth: 积分深度 (m)，默认 3000 m
/// bottom: 底层高度 (m)，可选，默认使用数组最小值
/// storm_u: 风暴移动速度的 u 分量 (m/s)，可选，默认 0
/// storm_v: 风暴移动速度的 v 分量 (m/s)，可选，默认 0
/// 返回: 风暴相对螺旋度 (m²/s²)
pub fn storm_relative_helicity(
  height : Array[Double],
  u : Array[Double],
  v : Array[Double],
  depth? : Double,
  bottom? : Double,
  storm_u? : Double,
  storm_v? : Double,
) -> Double? {
  if height.length() != u.length() ||
    height.length() != v.length() ||
    height.length() < 2 {
    None
  } else {
    let depth_val = match depth {
      Some(d) => d
      None => 3000.0 // 默认 3000 m
    }
    let bottom_val = match bottom {
      Some(b) => b
      None => height[0] // 最低高度（数组开头）
    }
    let top_val = bottom_val + depth_val
    let storm_u_val = match storm_u {
      Some(u) => u
      None => 0.0
    }
    let storm_v_val = match storm_v {
      Some(v) => v
      None => 0.0
    }

    // 找到积分区间
    let mut start_idx = 0
    let mut end_idx = height.length() - 1
    for i = 0; i < height.length(); i = i + 1 {
      if height[i] >= bottom_val && start_idx == 0 {
        start_idx = i
      }
      if height[i] <= top_val {
        end_idx = i
      }
    }
    if start_idx >= end_idx {
      None
    } else {
      // 计算相对风速
      let rel_u = Array::make(end_idx - start_idx + 1, 0.0)
      let rel_v = Array::make(end_idx - start_idx + 1, 0.0)
      let h_arr = Array::make(end_idx - start_idx + 1, 0.0)
      for i = start_idx; i <= end_idx; i = i + 1 {
        let idx = i - start_idx
        rel_u[idx] = u[i] - storm_u_val
        rel_v[idx] = v[i] - storm_v_val
        h_arr[idx] = height[i]
      }

      // 计算螺旋度
      // H = ∫ (v_rel · (∇ × v_rel)) dz
      // 对于一维垂直剖面，简化计算
      // H ≈ Σ (v_rel[i] · (u_rel[i+1] - u_rel[i]) - u_rel[i] · (v_rel[i+1] - v_rel[i])) * Δz
      let mut helicity = 0.0
      for i = 0; i < rel_u.length() - 1; i = i + 1 {
        let dz = h_arr[i + 1] - h_arr[i]
        if dz > 0.0 {
          // 计算相对涡度项
          let du_dz = (rel_u[i + 1] - rel_u[i]) / dz
          let dv_dz = (rel_v[i + 1] - rel_v[i]) / dz

          // 螺旋度贡献: v · (du/dz) - u · (dv/dz)
          let contrib = rel_v[i] * du_dz - rel_u[i] * dv_dz
          helicity = helicity + contrib * dz
        }
      }
      Some(helicity)
    }
  }
}

// ============================================================================
// 13. 体感温度函数
// ============================================================================

///|
/// 计算风寒指数 (Wind Chill)
/// 基于 NWS 风寒指数公式（2001年版本）
/// temperature: 环境温度 (K)
/// speed: 风速 (m/s)
/// face_level_winds: 是否使用面部高度风速（默认 false，使用 10 m 高度）
/// mask_undefined: 是否对超出有效范围的值返回 None（默认 true）
/// 返回: 风寒指数 (K)
/// 注意: 仅当温度 < 10°C (283.15 K) 且风速 > 4.8 km/h (1.33 m/s) 时有效
pub fn windchill(
  temperature : Double,
  speed : Double,
  _face_level_winds? : Bool,
  mask_undefined? : Bool,
) -> Double? {
  let mask = match mask_undefined {
    Some(m) => m
    None => true
  }

  // 转换为摄氏度用于计算
  let temp_c = temperature - 273.15

  // 风速转换为 km/h
  let speed_kmh = speed * 3.6

  // 有效范围检查
  let is_valid = temp_c < 10.0 && speed_kmh > 4.8
  if !is_valid && mask {
    None
  } else {
    // NWS 风寒指数公式（2001年）
    // WC = 13.12 + 0.6215*T - 11.37*V^0.16 + 0.3965*T*V^0.16
    // 其中 T 是温度（°C），V 是风速（km/h）
    let v_pow = Double::pow(speed_kmh, 0.16)
    let wc_c = 13.12 + 0.6215 * temp_c - 11.37 * v_pow + 0.3965 * temp_c * v_pow

    // 转换回开尔文
    Some(wc_c + 273.15)
  }
}

///|
/// 计算热指数 (Heat Index)
/// 基于 Rothfusz 公式（NWS 使用）
/// temperature: 环境温度 (K)
/// relative_humidity: 相对湿度 (0-1 或 0-100)
/// mask_undefined: 是否对超出有效范围的值返回 None（默认 true）
/// 返回: 热指数 (K)
/// 注意: 仅当温度 > 27°C (300.15 K) 且相对湿度 > 40% 时有效
pub fn heat_index(
  temperature : Double,
  relative_humidity : Double,
  mask_undefined? : Bool,
) -> Double? {
  let mask = match mask_undefined {
    Some(m) => m
    None => true
  }

  // 转换为摄氏度用于计算
  let temp_c = temperature - 273.15

  // 处理相对湿度：如果 > 1，假设是百分比形式
  let rh = if relative_humidity > 1.0 {
    relative_humidity / 100.0
  } else {
    relative_humidity
  }

  // 有效范围检查
  let is_valid = temp_c > 27.0 && rh > 0.4
  if !is_valid && mask {
    None
  } else {
    // Rothfusz 热指数公式
    // HI = c1 + c2*T + c3*RH + c4*T*RH + c5*T² + c6*RH² + c7*T²*RH + c8*T*RH² + c9*T²*RH²
    let c1 : Double = -8.78469475556
    let c2 : Double = 1.61139411
    let c3 : Double = 2.33854883889
    let c4 : Double = -0.14611605
    let c5 : Double = -0.012308094
    let c6 : Double = -0.0164248277778
    let c7 : Double = 0.002211732
    let c8 : Double = 0.00072546
    let c9 : Double = -0.000003582
    let t2 = temp_c * temp_c
    let rh2 = rh * rh
    let hi_c = c1 +
      c2 * temp_c +
      c3 * rh +
      c4 * temp_c * rh +
      c5 * t2 +
      c6 * rh2 +
      c7 * t2 * rh +
      c8 * temp_c * rh2 +
      c9 * t2 * rh2

    // 转换回开尔文
    Some(hi_c + 273.15)
  }
}

// ============================================================================
// 14. 新增功能函数
// ============================================================================

///|
/// 计算表观温度 (Apparent Temperature)
/// 根据条件自动选择风寒指数或热指数
/// temperature: 环境温度 (K)
/// relative_humidity: 相对湿度 (0-1 或 0-100)
/// speed: 风速 (m/s)
/// face_level_winds: 是否使用面部高度风速（默认 false）
/// mask_undefined: 是否对超出有效范围的值返回 None（默认 true）
pub fn apparent_temperature(
  temperature : Double,
  relative_humidity : Double,
  speed : Double,
  face_level_winds? : Bool,
  mask_undefined? : Bool,
) -> Double? {
  let mask = match mask_undefined {
    Some(m) => m
    None => true
  }

  // 计算风寒指数和热指数
  let face_level_winds_val = match face_level_winds {
    Some(flw) => flw
    None => false
  }
  let wc_result = windchill(
    temperature,
    speed,
    _face_level_winds=face_level_winds_val,
    mask_undefined=mask,
  )
  let hi_result = heat_index(
    temperature,
    relative_humidity,
    mask_undefined=mask,
  )

  // 如果风寒指数有效，使用风寒指数；否则使用热指数
  match (wc_result, hi_result) {
    (Some(wc), _) => Some(wc)
    (None, Some(hi)) => Some(hi)
    (None, None) =>
      if mask {
        None
      } else {
        Some(temperature) // 如果都无效且不mask，返回原温度
      }
  }
}

///|
/// 计算干绝热过程中的温度
/// 气块干绝热抬升时的温度变化
/// pressure: 目标气压 (Pa)
/// temperature: 起始温度 (K)
/// reference_pressure: 参考气压 (Pa)，默认使用 pressure
/// 返回: 目标气压下的温度 (K)
pub fn dry_lapse(
  pressure : Double,
  temperature : Double,
  reference_pressure? : Double,
) -> Double {
  let ref_p = match reference_pressure {
    Some(p) => p
    None => pressure
  }
  // 干绝热过程：T = T_ref * (p / p_ref)^κ
  temperature * Double::pow(pressure / ref_p, poisson_exponent)
}

///|
/// 计算对流凝结高度 (CCL - Convective Condensation Level)
/// 使用地面混合比线与环境温度廓线的交点
/// pressure: 气压数组 (Pa)，从高到低排列
/// temperature: 环境温度数组 (K)
/// dewpoint: 露点温度数组 (K)
/// which: "top" (默认) 或 "bottom"，返回哪个 CCL
/// 返回: (CCL 气压 (Pa), CCL 温度 (K), 对流温度 (K))
pub fn ccl(
  pressure : Array[Double],
  temperature : Array[Double],
  dewpoint_arr : Array[Double],
  which? : String,
) -> (Double, Double, Double)? {
  if pressure.length() != temperature.length() ||
    pressure.length() != dewpoint_arr.length() ||
    pressure.length() < 2 {
    None
  } else {
    let which_val = match which {
      Some(w) => w
      None => "top"
    }

    // 使用地面露点计算起始混合比
    let p_start = pressure[0]
    let dewpoint_start = dewpoint_arr[0]
    let es_start = saturation_vapor_pressure(dewpoint_start, "liquid")
    let r_start = mixing_ratio(es_start, p_start)

    // 计算沿等混合比线的温度（即露点温度）
    fn calc_rt(i : Int) -> Double {
      let e = vapor_pressure(pressure[i], r_start)
      dewpoint(e)
    }

    let rt_profile = Array::makei(pressure.length(), calc_rt)

    // 找到等混合比线与环境温度线的交点
    // 简化实现：找到第一个交点
    let mut found_intersection = false
    let mut p_ccl = 0.0
    let mut t_ccl = 0.0
    for i = 0; i < pressure.length() - 1; i = i + 1 {
      if !found_intersection {
        let rt1 = rt_profile[i]
        let rt2 = rt_profile[i + 1]
        let t1 = temperature[i]
        let t2 = temperature[i + 1]

        // 如果等混合比线从一侧穿过环境温度线到另一侧
        if (rt1 > t1 && rt2 <= t2) || (rt1 <= t1 && rt2 > t2) {
          // 线性插值找到交点
          // 解方程: rt1 + ratio * (rt2 - rt1) = t1 + ratio * (t2 - t1)
          // ratio = (t1 - rt1) / ((rt2 - rt1) - (t2 - t1))
          let p1 = pressure[i]
          let p2 = pressure[i + 1]
          let denominator = rt2 - rt1 - (t2 - t1)
          let ratio = if denominator != 0.0 {
            (t1 - rt1) / denominator
          } else {
            0.5
          }
          // 限制 ratio 在 [0, 1] 范围内
          let ratio_clamped = if ratio < 0.0 {
            0.0
          } else if ratio > 1.0 {
            1.0
          } else {
            ratio
          }
          p_ccl = p1 + ratio_clamped * (p2 - p1)
          t_ccl = t1 + ratio_clamped * (t2 - t1)
          found_intersection = true
          // 如果 which 是 "top"，继续查找最后一个交点
          if which_val != "top" {
            break
          }
        }
      } else if which_val == "top" {
        // 继续查找最后一个交点
        let rt1 = rt_profile[i]
        let rt2 = rt_profile[i + 1]
        let t1 = temperature[i]
        let t2 = temperature[i + 1]
        if (rt1 > t1 && rt2 <= t2) || (rt1 <= t1 && rt2 > t2) {
          let p1 = pressure[i]
          let p2 = pressure[i + 1]
          let denominator = rt2 - rt1 - (t2 - t1)
          let ratio = if denominator != 0.0 {
            (t1 - rt1) / denominator
          } else {
            0.5
          }
          let ratio_clamped = if ratio < 0.0 {
            0.0
          } else if ratio > 1.0 {
            1.0
          } else {
            ratio
          }
          p_ccl = p1 + ratio_clamped * (p2 - p1)
          t_ccl = t1 + ratio_clamped * (t2 - t1)
        }
      }
    }
    if !found_intersection {
      None
    } else {
      // 计算对流温度（从地面干绝热抬升到 CCL 的温度）
      let t_convective = dry_lapse(p_start, t_ccl, reference_pressure=p_ccl)
      Some((p_ccl, t_ccl, t_convective))
    }
  }
}

///|
/// 标准大气压转高度
/// 使用 U.S. 标准大气公式
/// pressure: 气压 (Pa)
/// 返回: 高度 (m)
pub fn pressure_to_height_std(pressure : Double) -> Double {
  // 标准大气常数
  let t0 : Double = 288.0 // K
  let p0 : Double = 101325.0 // Pa (1013.25 hPa)
  let gamma : Double = 0.0065 // K/m (6.5 K/km)

  // 公式: Z = (T0 / γ) * [1 - (p/p0)^(R*γ/g)]
  let r = dry_air_gas_constant
  let g = earth_gravity
  let exponent = r * gamma / g
  t0 / gamma * (1.0 - Double::pow(pressure / p0, exponent))
}

///|
/// 高度转标准大气压
/// 使用 U.S. 标准大气公式的逆函数
/// height: 高度 (m)
/// 返回: 气压 (Pa)
pub fn height_to_pressure_std(height : Double) -> Double {
  // 标准大气常数
  let t0 : Double = 288.0 // K
  let p0 : Double = 101325.0 // Pa
  let gamma : Double = 0.0065 // K/m

  // 公式: p = p0 * [1 - (γ*Z/T0)]^(g/(R*γ))
  let r = dry_air_gas_constant
  let g = earth_gravity
  let exponent = g / (r * gamma)
  p0 * Double::pow(1.0 - gamma * height / t0, exponent)
}

///|
/// 计算体积切变 (Bulk Shear)
/// 计算指定层内的风切变
/// pressure: 气压数组 (Pa)，从高到低排列
/// u: u 分量数组 (m/s)
/// v: v 分量数组 (m/s)
/// depth: 层厚度 (Pa)，默认 10000 Pa (100 hPa)
/// bottom: 底层气压 (Pa)，默认使用数组最大值
/// 返回: (u 切变 (m/s), v 切变 (m/s))
pub fn bulk_shear(
  pressure : Array[Double],
  u : Array[Double],
  v : Array[Double],
  depth? : Double,
  bottom? : Double,
) -> (Double, Double)? {
  if pressure.length() != u.length() ||
    pressure.length() != v.length() ||
    pressure.length() < 2 {
    None
  } else {
    let depth_val = match depth {
      Some(d) => d
      None => 10000.0 // 100 hPa
    }
    let bottom_val = match bottom {
      Some(b) => b
      None => {
        // 找到最大气压（数组末尾）
        let mut max_p = pressure[0]
        for i = 0; i < pressure.length(); i = i + 1 {
          if pressure[i] > max_p {
            max_p = pressure[i]
          }
        }
        max_p
      }
    }
    let top_val = bottom_val - depth_val

    // 找到底层和顶层的索引
    let mut bottom_idx = 0
    let mut top_idx = pressure.length() - 1
    for i = 0; i < pressure.length(); i = i + 1 {
      if Double::abs(pressure[i] - bottom_val) <
        Double::abs(pressure[bottom_idx] - bottom_val) {
        bottom_idx = i
      }
      if Double::abs(pressure[i] - top_val) <
        Double::abs(pressure[top_idx] - top_val) {
        top_idx = i
      }
    }

    // 计算切变 = 顶层风 - 底层风
    let u_shear = u[top_idx] - u[bottom_idx]
    let v_shear = v[top_idx] - v[bottom_idx]
    Some((u_shear, v_shear))
  }
}

///|
/// 计算 Bunkers 风暴移动向量
/// 计算右移和左移超级单体风暴的移动方向
/// pressure: 气压数组 (Pa)
/// u: u 分量数组 (m/s)
/// v: v 分量数组 (m/s)
/// height: 高度数组 (m)
/// 返回: ((右移 u, 右移 v), (左移 u, 左移 v), (平均 u, 平均 v))
pub fn bunkers_storm_motion(
  pressure : Array[Double],
  u : Array[Double],
  v : Array[Double],
  height : Array[Double],
) -> ((Double, Double), (Double, Double), (Double, Double))? {
  if pressure.length() != u.length() ||
    pressure.length() != v.length() ||
    pressure.length() != height.length() ||
    pressure.length() < 2 {
    None
  } else {
    // 找到 0-6km 层的平均风
    let mut u_sum = 0.0
    let mut v_sum = 0.0
    let mut count = 0
    let h6km = 6000.0
    for i = 0; i < height.length(); i = i + 1 {
      if height[i] <= h6km {
        u_sum = u_sum + u[i]
        v_sum = v_sum + v[i]
        count = count + 1
      }
    }
    if count == 0 {
      None
    } else {
      let u_mean = u_sum / Double::from_int(count)
      let v_mean = v_sum / Double::from_int(count)

      // 找到 0-0.5km 层的平均风
      let mut u_500m = 0.0
      let mut v_500m = 0.0
      let mut count_500m = 0
      let h500m = 500.0
      for i = 0; i < height.length(); i = i + 1 {
        if height[i] <= h500m {
          u_500m = u_500m + u[i]
          v_500m = v_500m + v[i]
          count_500m = count_500m + 1
        }
      }
      if count_500m == 0 {
        None
      } else {
        u_500m = u_500m / Double::from_int(count_500m)
        v_500m = v_500m / Double::from_int(count_500m)

        // 找到 5.5-6km 层的平均风
        let mut u_5500m = 0.0
        let mut v_5500m = 0.0
        let mut count_5500m = 0
        let h5500m = 5500.0
        for i = 0; i < height.length(); i = i + 1 {
          if height[i] >= h5500m && height[i] <= h6km {
            u_5500m = u_5500m + u[i]
            v_5500m = v_5500m + v[i]
            count_5500m = count_5500m + 1
          }
        }
        if count_5500m == 0 {
          None
        } else {
          u_5500m = u_5500m / Double::from_int(count_5500m)
          v_5500m = v_5500m / Double::from_int(count_5500m)

          // 计算切变向量
          let shear_u = u_5500m - u_500m
          let shear_v = v_5500m - v_500m
          let shear_mag = Double::sqrt(shear_u * shear_u + shear_v * shear_v)
          if shear_mag < 0.001 {
            None
          } else {
            // 计算偏差向量（垂直于切变，幅度 7.5 m/s）
            let dev_mag = 7.5
            let rdev_u = shear_v * (dev_mag / shear_mag)
            let rdev_v = -shear_u * (dev_mag / shear_mag)

            // 右移：平均风 + 偏差
            let rm_u = u_mean + rdev_u
            let rm_v = v_mean + rdev_v

            // 左移：平均风 - 偏差
            let lm_u = u_mean - rdev_u
            let lm_v = v_mean - rdev_v
            Some(((rm_u, rm_v), (lm_u, lm_v), (u_mean, v_mean)))
          }
        }
      }
    }
  }
}

///|
/// 计算超级单体复合参数
/// mucape: 最不稳定 CAPE (J/kg)
/// effective_storm_helicity: 有效层风暴相对螺旋度 (m²/s²)
/// effective_shear: 有效切变 (m/s)
/// 返回: 超级单体复合参数（无量纲）
pub fn supercell_composite(
  mucape : Double,
  effective_storm_helicity : Double,
  effective_shear : Double,
) -> Double {
  // 公式: SCP = (MUCAPE/1000) * (ESRH/50) * (EShear/20)
  // 其中 EShear 在 10 m/s 以下设为 0，20 m/s 以上设为 1
  let eshear_term = if effective_shear < 10.0 {
    0.0
  } else if effective_shear > 20.0 {
    1.0
  } else {
    effective_shear / 20.0
  }
  mucape / 1000.0 * (effective_storm_helicity / 50.0) * eshear_term
}

///|
/// 计算显著龙卷风参数
/// sbcape: 地面基础 CAPE (J/kg)
/// surface_based_lcl_height: 地面基础 LCL 高度 (m)
/// storm_helicity_1km: 0-1km 风暴相对螺旋度 (m²/s²)
/// shear_6km: 0-6km 切变 (m/s)
/// 返回: 显著龙卷风参数（无量纲）
pub fn significant_tornado(
  sbcape : Double,
  surface_based_lcl_height : Double,
  storm_helicity_1km : Double,
  shear_6km : Double,
) -> Double {
  // 公式: SIGTOR = (SBCAPE/1500) * ((2000-LCL)/1000) * (SRH1km/150) * (Shear6km/20)
  // LCL 高度项：>2000m 设为 0，<1000m 设为 1
  let lcl_term = if surface_based_lcl_height > 2000.0 {
    0.0
  } else if surface_based_lcl_height < 1000.0 {
    1.0
  } else {
    (2000.0 - surface_based_lcl_height) / 1000.0
  }

  // 切变项：<12.5 m/s 设为 0，>30 m/s 设为 1.5
  let shear_term = if shear_6km < 12.5 {
    0.0
  } else if shear_6km > 30.0 {
    1.5
  } else {
    shear_6km / 20.0
  }
  sbcape / 1500.0 * lcl_term * (storm_helicity_1km / 150.0) * shear_term
}

///|
/// 计算切变涡度
/// 公式: ζ_shear = (∂v/∂x - ∂u/∂y) 的切变部分
/// u: 东西方向风速 (m/s)，二维数组
/// v: 南北方向风速 (m/s)，二维数组
/// dx: x 方向网格间距 (m)
/// dy: y 方向网格间距 (m)
/// 返回: 切变涡度数组 (1/s)
pub fn shear_vorticity(
  u : Array[Array[Double]],
  v : Array[Array[Double]],
  dx : Double,
  dy : Double,
) -> Array[Array[Double]]? {
  // 切变涡度与总涡度相同，因为这是二维水平涡度
  vorticity(u, v, dx, dy)
}

///|
/// 计算曲率涡度
/// 公式: ζ_curvature = (v/r) - (u/r) * (∂u/∂x + ∂v/∂y)
/// 其中 r 是曲率半径（简化计算）
/// u: 东西方向风速 (m/s)，二维数组
/// v: 南北方向风速 (m/s)，二维数组
/// dx: x 方向网格间距 (m)
/// dy: y 方向网格间距 (m)
/// 返回: 曲率涡度数组 (1/s)
pub fn curvature_vorticity(
  u : Array[Array[Double]],
  v : Array[Array[Double]],
  dx : Double,
  dy : Double,
) -> Array[Array[Double]]? {
  // 简化实现：曲率涡度通常需要流线曲率信息
  // 这里使用近似方法
  if u.length() != v.length() || u.length() == 0 {
    None
  } else {
    let ny = u.length()
    let nx = if u[0].length() > 0 { u[0].length() } else { 0 }
    if nx == 0 || v.length() == 0 || v[0].length() != nx {
      None
    } else {
      let curv_vort = Array::make(ny, Array::make(nx, 0.0))

      // 计算风速
      for j = 0; j < ny; j = j + 1 {
        for i = 0; i < nx; i = i + 1 {
          let speed = Double::sqrt(u[j][i] * u[j][i] + v[j][i] * v[j][i])

          // 计算风向变化率作为曲率涡度的近似
          // 这是一个简化实现
          let du_dx = if i > 0 && i < nx - 1 {
            (u[j][i + 1] - u[j][i - 1]) / (2.0 * dx)
          } else if i == 0 {
            (u[j][i + 1] - u[j][i]) / dx
          } else {
            (u[j][i] - u[j][i - 1]) / dx
          }
          let dv_dy = if j > 0 && j < ny - 1 {
            (v[j + 1][i] - v[j - 1][i]) / (2.0 * dy)
          } else if j == 0 {
            (v[j + 1][i] - v[j][i]) / dy
          } else {
            (v[j][i] - v[j - 1][i]) / dy
          }

          // 曲率涡度的简化计算（基于速度梯度的组合）
          if speed > 0.001 {
            curv_vort[j][i] = (v[j][i] / dx - u[j][i] / dy) *
              (du_dx + dv_dy) /
              speed
          } else {
            curv_vort[j][i] = 0.0
          }
        }
      }
      Some(curv_vort)
    }
  }
}

///|
/// 计算平流
/// 计算标量场沿着风场的平流
/// field: 标量场（二维数组）
/// u: 东西方向风速 (m/s)，二维数组
/// v: 南北方向风速 (m/s)，二维数组
/// dx: x 方向网格间距 (m)
/// dy: y 方向网格间距 (m)
/// 返回: 平流数组（与 field 相同单位/s）
pub fn advection(
  field : Array[Array[Double]],
  u : Array[Array[Double]],
  v : Array[Array[Double]],
  dx : Double,
  dy : Double,
) -> Array[Array[Double]]? {
  if field.length() != u.length() ||
    field.length() != v.length() ||
    field.length() == 0 {
    None
  } else {
    let ny = field.length()
    let nx = if field[0].length() > 0 { field[0].length() } else { 0 }
    if nx == 0 || u[0].length() != nx || v[0].length() != nx {
      None
    } else {
      let adv = Array::make(ny, Array::make(nx, 0.0))

      // 平流公式: -u * (∂field/∂x) - v * (∂field/∂y)
      for j = 0; j < ny; j = j + 1 {
        for i = 0; i < nx; i = i + 1 {
          // ∂field/∂x
          let df_dx = if i > 0 && i < nx - 1 {
            (field[j][i + 1] - field[j][i - 1]) / (2.0 * dx)
          } else if i == 0 {
            (field[j][i + 1] - field[j][i]) / dx
          } else {
            (field[j][i] - field[j][i - 1]) / dx
          }

          // ∂field/∂y
          let df_dy = if j > 0 && j < ny - 1 {
            (field[j + 1][i] - field[j - 1][i]) / (2.0 * dy)
          } else if j == 0 {
            (field[j + 1][i] - field[j][i]) / dy
          } else {
            (field[j][i] - field[j - 1][i]) / dy
          }

          // 平流 = -u * df/dx - v * df/dy
          adv[j][i] = -u[j][i] * df_dx - v[j][i] * df_dy
        }
      }
      Some(adv)
    }
  }
}

///|
/// 计算锋生函数
/// 计算温度梯度的时间变化率
/// temperature: 温度场 (K)，二维数组
/// u: 东西方向风速 (m/s)，二维数组
/// v: 南北方向风速 (m/s)，二维数组
/// dx: x 方向网格间距 (m)
/// dy: y 方向网格间距 (m)
/// 返回: 锋生函数数组 (K/(m·s))
pub fn frontogenesis(
  temperature : Array[Array[Double]],
  u : Array[Array[Double]],
  v : Array[Array[Double]],
  dx : Double,
  dy : Double,
) -> Array[Array[Double]]? {
  if temperature.length() != u.length() ||
    temperature.length() != v.length() ||
    temperature.length() == 0 {
    None
  } else {
    let ny = temperature.length()
    let nx = if temperature[0].length() > 0 {
      temperature[0].length()
    } else {
      0
    }
    if nx == 0 || u[0].length() != nx || v[0].length() != nx {
      None
    } else {
      let fronto = Array::make(ny, Array::make(nx, 0.0))

      // 锋生函数: F = -0.5 * |∇T|^(-1) * [∇T · ∇(dV/dt)]
      // 简化实现：计算温度梯度的平流
      for j = 0; j < ny; j = j + 1 {
        for i = 0; i < nx; i = i + 1 {
          // 计算温度梯度
          let dT_dx = if i > 0 && i < nx - 1 {
            (temperature[j][i + 1] - temperature[j][i - 1]) / (2.0 * dx)
          } else if i == 0 {
            (temperature[j][i + 1] - temperature[j][i]) / dx
          } else {
            (temperature[j][i] - temperature[j][i - 1]) / dx
          }
          let dT_dy = if j > 0 && j < ny - 1 {
            (temperature[j + 1][i] - temperature[j - 1][i]) / (2.0 * dy)
          } else if j == 0 {
            (temperature[j + 1][i] - temperature[j][i]) / dy
          } else {
            (temperature[j][i] - temperature[j - 1][i]) / dy
          }
          let grad_T_mag = Double::sqrt(dT_dx * dT_dx + dT_dy * dT_dy)
          if grad_T_mag > 0.001 {
            // 计算速度梯度
            let du_dx = if i > 0 && i < nx - 1 {
              (u[j][i + 1] - u[j][i - 1]) / (2.0 * dx)
            } else if i == 0 {
              (u[j][i + 1] - u[j][i]) / dx
            } else {
              (u[j][i] - u[j][i - 1]) / dx
            }
            let du_dy = if j > 0 && j < ny - 1 {
              (u[j + 1][i] - u[j - 1][i]) / (2.0 * dy)
            } else if j == 0 {
              (u[j + 1][i] - u[j][i]) / dy
            } else {
              (u[j][i] - u[j - 1][i]) / dy
            }
            let dv_dx = if i > 0 && i < nx - 1 {
              (v[j][i + 1] - v[j][i - 1]) / (2.0 * dx)
            } else if i == 0 {
              (v[j][i + 1] - v[j][i]) / dx
            } else {
              (v[j][i] - v[j][i - 1]) / dx
            }
            let dv_dy = if j > 0 && j < ny - 1 {
              (v[j + 1][i] - v[j - 1][i]) / (2.0 * dy)
            } else if j == 0 {
              (v[j + 1][i] - v[j][i]) / dy
            } else {
              (v[j][i] - v[j - 1][i]) / dy
            }

            // 锋生函数: F = -0.5/|∇T| * [dT/dx * (du/dx * dT/dx + dv/dx * dT/dy) + 
            //                            dT/dy * (du/dy * dT/dx + dv/dy * dT/dy)]
            let term1 = dT_dx * (du_dx * dT_dx + dv_dx * dT_dy)
            let term2 = dT_dy * (du_dy * dT_dx + dv_dy * dT_dy)
            fronto[j][i] = -0.5 / grad_T_mag * (term1 + term2)
          } else {
            fronto[j][i] = 0.0
          }
        }
      }
      Some(fronto)
    }
  }
}