-- TOFCalculator.lua
-- 完整版：助推/二级/滑翔全程数值积分 + 30° loft + 空气阻力模型
-- 性能参数表已补全

local M = {}
local g0      = 9.81          -- 重力加速度 (m/s²)
local rho0    = 1.225         -- 海平面空气密度 (kg/m³)
local H_scale = 7640          -- 大气标度高度 (m)
local Cd      = 0.15          -- 阻力系数
local A       = 0.126         -- 弹体迎风面积 (m²)

-- 标准大气密度随高度衰减
local function airDensity(h)
  if h < 0 then h = 0 end
  return rho0 * math.exp(-h / H_scale)
end

-- 解 At² + 2B t − 2R = 0 的正根：R = B·t + ½A·t²
local function solveTime(A, B, R)
  if A == 0 then return R / B end
  local disc = B*B + 2*A*R
  return (-B + math.sqrt(disc)) / A
end

-- loft sinθ 线性插值，最大 30°
local function getLoftSin(p, R0)
  if R0 <= p.loft_off_range then
    return 0
  elseif R0 < p.loft_min_range then
    local t = (R0 - p.loft_off_range)
            / (p.loft_min_range - p.loft_off_range)
    return p.loft_sin * t
  end
  return p.loft_sin
end

-- 导弹性能参数表
local perf = {
  --------------------------------------------------------------------------
  -- PL-10E（短程近防，单级助推，无 loft）
  --------------------------------------------------------------------------
  ["{PLAAF_PL-10E}"] = {
    loft_active      = false,
    loft_min_range   = 0,
    loft_off_range   = 0,
    loft_sin         = 0,

    boost_Isp        = 240.0,    -- s
    boost_thrust     = 12802.0,  -- N
    fuel_mass_boost  = 27.2,     -- kg

    march_Isp        = nil,
    march_thrust     = nil,
    fuel_mass_march  = 0.0,

    mass_total       = 84.46,    -- kg
  },

  --------------------------------------------------------------------------
  -- PL-15（中/远程，多用途，二级助推 + loft）
  --------------------------------------------------------------------------
  ["{J20A_PL15}"] = {
    loft_active      = true,
    loft_min_range   = 30000.0,
    loft_off_range   = 18000.0,
    loft_sin         = math.sin(math.rad(30)),  -- 30°

    boost_Isp        = 242.0,
    boost_thrust     = 22770.0,
    fuel_mass_boost  = 70.6,

    march_Isp        = 236.0,
    march_thrust     = 6530.0,
    fuel_mass_march  = 30.2,

    mass_total       = 210.0,    -- kg
  },

  --------------------------------------------------------------------------
  -- PL-17（远程超视距，二级助推 + loft）
  --------------------------------------------------------------------------
  ["{J20A_PL17}"] = {
    loft_active      = true,
    loft_min_range   = 30000.0,
    loft_off_range   = 18000.0,
    loft_sin         = math.sin(math.rad(30)),  -- 30°

    boost_Isp        = 472.0,
    boost_thrust     = 16000.0,
    fuel_mass_boost  = 70.6,

    march_Isp        = 336.0,
    march_thrust     = 4700.0,
    fuel_mass_march  = 30.2,

    mass_total       = 300.0,    -- kg
  },

  --------------------------------------------------------------------------
  -- PL-12（中程，多用途，单级助推 + loft）
  --------------------------------------------------------------------------
  ["DIS_PL-12"] = {
    loft_active      = true,
    loft_min_range   = 30000.0,
    loft_off_range   = 18000.0,
    loft_sin         = math.sin(math.rad(30)),  -- 30°

    boost_Isp        = 242.0,
    boost_thrust     = 22770.0,
    fuel_mass_boost  = 70.6,

    march_Isp        = nil,
    march_thrust     = nil,
    fuel_mass_march  = 0.0,

    mass_total       = 175.0,    -- kg
  },
}

-- 兼容/复用条目
perf["{BayArm_PLAAF_PL-10E_x4}"] = perf["{PLAAF_PL-10E}"]
perf["{PLAAF_PL-10E_DUAL}"]     = perf["{PLAAF_PL-10E}"]
perf["J20_PL-15_DUAL"]          = perf["{J20A_PL15}"]
perf["BayArm_J20_PL-15_DUAL"]   = perf["{J20A_PL15}"]
perf["{BayArm_J20_PL-17_x4}"]   = perf["{J20A_PL17}"]
perf["DIS_PL-12_DUAL_L"]        = perf["DIS_PL-12"]

-- TOF 计算主函数
-- host: { vel = m/s, alt = m }
-- tgt : { range = m, rdot = m/s }
function M.calculateTOF(clsid, host, tgt)
    -- MOD1: 如果 clsid 没在 perf 表中，打印提示并返回 0
    if not perf[clsid] then
        print(string.format("TOF 无配置：clsid='%s'", tostring(clsid)))
        return 0
    end
    local p = perf[clsid]

    -- 初始化参数
    local R_rem = tgt.range or 0
    local Vt    = tgt.rdot  or 0
    local v0    = host.vel  or 0
    local h     = host.alt  or 0
    local tof   = 0

    -- 1) 计算抛射角（loft）
    local sin_t = 0
    if p.loft_active and R_rem >= p.loft_min_range then
        sin_t = getLoftSin(p, R_rem)
    end
    local cos_t = math.sqrt(1 - sin_t*sin_t)

    -- 2) Boost 段 数值积分
    local t_b        = p.boost_Isp
                       and (p.boost_Isp * g0 * p.fuel_mass_boost / p.boost_thrust)
                       or 0
    local steps      = 40
    local dt         = t_b / steps
    local mass0      = p.mass_total
    local mf_b       = p.fuel_mass_boost
    local thrust_b   = p.boost_thrust

    local v  = v0
    local s_h = 0

    for i = 1, steps do
        local frac = (i - 0.5) / steps
        local m_i  = mass0 - mf_b * frac
        local rho  = airDensity(h)
        local Fd   = 0.5 * rho * Cd * A * v * v
        local a    = (thrust_b - Fd) / m_i

        local v_new = v + a * dt
        local v_avg = 0.5 * (v + v_new)
        local ds    = v_avg * dt
        local ds_h  = ds * cos_t
        local ds_v  = ds * sin_t

        if s_h + ds_h >= R_rem then
            local Rp    = R_rem - s_h
            local A_p   = a * cos_t
            local B_p   = v * cos_t - Vt
            local t_hit = solveTime(A_p, B_p, Rp)
            return tof + (i - 1) * dt + t_hit
        end

        s_h = s_h + ds_h
        v   = v_new
        h   = h + ds_v
    end

    tof   = tof + t_b
    R_rem = R_rem - s_h
    local v1_h = v * cos_t

    -- 3) Coast-1 （此处可插 loft 重点推进段后处理，原示例未实现）

    -- 4) March 数值积分
    -- MOD2: 如果没有巡航推力/推进段，就跳过 March 阶段
    local v2_h = v1_h
    if p.march_Isp and p.march_thrust and p.fuel_mass_march > 0 then
        local t_m       = p.march_Isp * g0 * p.fuel_mass_march / p.march_thrust
        local mf_m      = p.fuel_mass_march
        local thrust_m  = p.march_thrust
        local steps2    = 30
        local dt2       = t_m / steps2

        local v_m = v1_h
        local s_m = 0

        for i = 1, steps2 do
            local frac2 = (i - 0.5) / steps2
            local m2    = mass0 - mf_b - mf_m * frac2

            local rho2 = airDensity(h)
            local Fd2  = 0.5 * rho2 * Cd * A * v_m * v_m
            local a_m  = (thrust_m - Fd2) / m2

            local v2_new = v_m + a_m * dt2
            local v2_avg = 0.5 * (v_m + v2_new)
            local ds2    = v2_avg * dt2

            if s_m + ds2 >= R_rem then
                local Rp2    = R_rem - s_m
                local A2     = a_m
                local B2     = v_m - Vt
                local t_hit2 = solveTime(A2, B2, Rp2)
                return tof + t_m * ((i - 1) / steps2) + t_hit2
            end

            s_m = s_m + ds2
            v_m = v2_new
            h   = h + ds2 * sin_t
        end

        tof   = tof + t_m
        R_rem = R_rem - s_m
        v2_h  = v_m
    end

    -- 5) Coast-2 匀速滑翔 + 阻力减速
    local dt3 = 0.1
    local s3  = 0
    local v3  = v2_h

    while s3 < R_rem do
        local rho3 = airDensity(h)
        local Fd3  = 0.5 * rho3 * Cd * A * v3 * v3
        local a3   = -Fd3 / mass0

        local v3_new = v3 + a3 * dt3
        local v3_avg = 0.5 * (v3 + v3_new)
        local ds3    = v3_avg * dt3

        if s3 + ds3 >= R_rem then
            local Rp3 = R_rem - s3
            return tof + Rp3 / v3_avg
        end

        s3  = s3 + ds3
        v3  = v3_new
        h   = h    -- 高度微小变化忽略
        tof = tof + dt3
    end

    return tof
end

--[[function M.calculateTOF(clsid, host, tgt)
  local p = perf[clsid]
  if not p then return 0 end

  local R_rem = tgt.range or 0
  local Vt    = tgt.rdot  or 0
  local v0    = host.vel  or 0
  local h     = host.alt  or 0
  local tof   = 0

  -- 1) 计算抛射角
  local sin_t = 0
  if p.loft_active and R_rem >= p.loft_min_range then
    sin_t = getLoftSin(p, R_rem)
  end
  local cos_t = math.sqrt(1 - sin_t*sin_t)

  ----------------------------------------------------------------
  -- 2) Boost 段 数值积分 (考虑质量变化 + 阻力)
  ----------------------------------------------------------------
  local t_b   = p.boost_Isp
             and (p.boost_Isp * g0 * p.fuel_mass_boost / p.boost_thrust)
             or 0
  local steps = 40
  local dt    = t_b / steps
  local mass0 = p.mass_total
  local mf_b  = p.fuel_mass_boost
  local thrust_b = p.boost_thrust

  local v     = v0
  local s_h   = 0

  for i = 1, steps do
    local frac = (i - 0.5) / steps
    local m_i  = mass0 - mf_b * frac
    local rho  = airDensity(h)
    local Fd   = 0.5 * rho * Cd * A * v * v
    local a    = (thrust_b - Fd) / m_i

    local v_new = v + a * dt
    local v_avg = 0.5 * (v + v_new)
    local ds    = v_avg * dt
    local ds_h  = ds * cos_t
    local ds_v  = ds * sin_t

    if s_h + ds_h >= R_rem then
      local Rp   = R_rem - s_h
      local A_p  = a * cos_t
      local B_p  = v * cos_t - Vt
      local t_hit = solveTime(A_p, B_p, Rp)
      return tof + (i - 1) * dt + t_hit
    end

    s_h = s_h + ds_h
    v   = v_new
    h   = h   + ds_v
  end

  tof   = tof + t_b
  R_rem = R_rem - s_h
  local v1_h = v * cos_t

  ----------------------------------------------------------------
  -- 3) Coast-1 (紧接二级点火，t=0)
  ----------------------------------------------------------------

  ----------------------------------------------------------------
  -- 4) March 数值积分
  ----------------------------------------------------------------
  local t_m   = p.march_Isp
             and (p.march_Isp * g0 * p.fuel_mass_march / p.march_thrust)
             or 0
  local mf_m  = p.fuel_mass_march
  local thrust_m = p.march_thrust
  local steps2 = 30
  local dt2   = t_m / steps2

  local v_m   = v1_h
  local s_m   = 0

  for i = 1, steps2 do
    local frac2 = (i - 0.5) / steps2
    local m2    = mass0 - mf_b - mf_m * frac2

    local rho2  = airDensity(h)
    local Fd2   = 0.5 * rho2 * Cd * A * v_m * v_m
    local a_m   = (thrust_m - Fd2) / m2

    local v2_new = v_m + a_m * dt2
    local v2_avg = 0.5 * (v_m + v2_new)
    local ds2    = v2_avg * dt2

    if s_m + ds2 >= R_rem then
      local Rp2   = R_rem - s_m
      local A2    = a_m
      local B2    = v_m - Vt
      local t_hit2 = solveTime(A2, B2, Rp2)
      return tof + t_m * ((i - 1)/steps2) + t_hit2
    end

    s_m = s_m + ds2
    v_m = v2_new
    h   = h   + ds2 * sin_t
  end

  tof   = tof + t_m
  R_rem = R_rem - s_m
  local v2_h = v_m

  ----------------------------------------------------------------
  -- 5) Coast-2 匀速滑翔 + 阻力减速 (简单积分)
  ----------------------------------------------------------------
  local dt3 = 0.1
  local s3  = 0
  local v3  = v2_h

  while s3 < R_rem do
    local rho3 = airDensity(h)
    local Fd3  = 0.5 * rho3 * Cd * A * v3 * v3
    local a3   = -Fd3 / mass0

    local v3_new = v3 + a3 * dt3
    local v3_avg = 0.5*(v3 + v3_new)
    local ds3    = v3_avg * dt3

    if s3 + ds3 >= R_rem then
      local Rp3 = R_rem - s3
      return tof + Rp3 / v3_avg
    end

    s3  = s3 + ds3
    v3  = v3_new
    h   = h     -- 高度微小变化可忽略
    tof = tof + dt3
  end

  return tof
end]]

return M
