-- Numerical and contract regression tests for the public LuaCoolProp API.
--
-- The assertions use physical invariants and cross-check two CoolProp entry
-- points instead of freezing version-specific reference tables.  This keeps
-- the suite useful when scripts/update-coolprop.sh selects another release.

local lcp = require("luacoolprop")

local function finite(value)
  return type(value) == "number"
    and value == value
    and value ~= math.huge
    and value ~= -math.huge
end

local function assert_close(actual, expected, relative_tolerance, message)
  assert(finite(actual), message .. ": actual value is not finite")
  assert(finite(expected), message .. ": expected value is not finite")
  local scale = math.max(1, math.abs(actual), math.abs(expected))
  local difference = math.abs(actual - expected)
  assert(difference <= relative_tolerance * scale,
    string.format("%s: relative difference %.6g exceeds %.6g",
      message, difference / scale, relative_tolerance))
end

local function assert_error_contains(callback, needle, message)
  local ok, result = pcall(callback)
  assert(not ok, message .. ": expected an error")
  assert(tostring(result):find(needle, 1, true),
    message .. ": unexpected diagnostic: " .. tostring(result))
end

-- Measure the same midpoint defect used by the TS/HS quality samplers, but in
-- the logarithmic reduced-pressure coordinate employed near the critical
-- point.  This independently verifies that tightening `tolerance` really
-- refines the critical tail rather than merely inserting a fixed point list.
local function critical_quality_midpoint_error(kind, fluid, quality, pcrit,
    first, second)
  local epsilon_first = 1 - first.p / pcrit
  local epsilon_second = 1 - second.p / pcrit
  if first.critical or second.critical
      or epsilon_first <= 0 or epsilon_second <= 0 then
    return 0
  end
  local epsilon_midpoint = math.sqrt(epsilon_first * epsilon_second)
  local pressure = pcrit * (1 - epsilon_midpoint)
  local entropy = lcp.propsSI("S", "P", pressure, "Q", quality, fluid)
  local ordinate
  if kind == "TS" then
    ordinate = lcp.propsSI("T", "P", pressure, "Q", quality, fluid)
  else
    ordinate = lcp.propsSI("H", "P", pressure, "Q", quality, fluid) * 1e-3
  end
  local dx = entropy * 1e-3 - 0.5 * (first.x + second.x)
  local dy = (ordinate - 0.5 * (first.y + second.y)) * 0.01
  return math.sqrt(dx * dx + dy * dy)
end

-- TeX key values can contain indentation around a braced path.  The loader
-- trims only this outer whitespace and preserves spaces inside directory
-- names.  Use a syntactically distinct `/./` path so an environment fallback
-- cannot make this regression pass accidentally.
local configured_library = os.getenv("LUACOOLPROP_LIB")
  or os.getenv("COOLPROP_LIB")
local explicit_library = configured_library
if explicit_library then
  local directory, basename = explicit_library:match("^(.*)/([^/]+)$")
  if directory and basename then
    explicit_library = directory .. "/./" .. basename
  else
    directory, basename = explicit_library:match("^(.*)\\([^\\]+)$")
    if directory and basename then
      explicit_library = directory .. "\\.\\" .. basename
    end
  end
end
local loaded, library_path = lcp.load_library(
  explicit_library and (" \n\t" .. explicit_library .. "  \r\n") or nil)
assert(loaded ~= nil, "CoolProp shared library must load")
assert(type(library_path) == "string" and library_path ~= "",
  "loaded library path must be reported")
if explicit_library then
  assert(library_path == explicit_library,
    "explicit library paths must ignore surrounding whitespace")
end
assert(type(lcp.version()) == "string" and lcp.version() ~= "",
  "CoolProp version must be non-empty")

local pressure = 2.5e5
local quality = 1
local direct_h = lcp.propsSI("Hmass", "P", pressure, "Q", quality, "R134a")
assert(finite(direct_h), "PropsSI enthalpy must be finite")

local state
local state_ok, state_result = xpcall(function()
  state = lcp.AbstractState("HEOS", "R134a")
  state:update("PQ_INPUTS", pressure, quality)
  return {
    h = state:keyed_output("Hmass"),
    p = state:keyed_output("P"),
  }
end, debug.traceback)
if state ~= nil then
  state:free()
  state:free()
  state = nil
end
if not state_ok then
  error(state_result, 0)
end
assert_close(state_result.h, direct_h, 2e-10,
  "PropsSI and AbstractState enthalpy")
assert_close(state_result.p, pressure, 2e-10,
  "AbstractState pressure after PQ update")

local constants = lcp.fluid_constants("R134a")
assert(finite(constants.pcrit) and constants.pcrit > 1e6,
  "critical pressure must be physically plausible")
assert(finite(constants.ptriple) and constants.ptriple > 0,
  "triple-point pressure must be positive")
assert(constants.pcrit > constants.ptriple,
  "critical pressure must exceed triple-point pressure")
assert(finite(constants.Tcrit) and finite(constants.Ttriple)
  and constants.Tcrit > constants.Ttriple,
  "critical temperature must exceed triple-point temperature")

local ph = lcp.diagram.get_type("PH")
assert(type(ph) == "table", "PH diagram registry entry must exist")
for _, family in ipairs({"quality", "isotherm", "isentrope", "isochore"}) do
  assert(type(ph.families[family]) == "table",
    family .. " family must be registered")
  assert(type(ph.families[family].curve) == "function",
    family .. " curve generator must exist")
  assert(type(ph.families[family].plots) == "function",
    family .. " plot serializer must exist")
end

local pv = lcp.diagram.get_type("PV")
assert(type(pv) == "table", "PV diagram registry entry must exist")
for _, family in ipairs({"quality", "isotherm", "isentrope"}) do
  assert(type(pv.families[family]) == "table",
    "PV " .. family .. " family must be registered")
  assert(type(pv.families[family].curve) == "function",
    "PV " .. family .. " curve generator must exist")
  assert(type(pv.families[family].plots) == "function",
    "PV " .. family .. " plot serializer must exist")
end

local ts = lcp.diagram.get_type("TS")
assert(type(ts) == "table", "TS diagram registry entry must exist")
for _, family in ipairs({"quality", "isenthalp"}) do
  assert(type(ts.families[family]) == "table",
    "TS " .. family .. " family must be registered")
  assert(type(ts.families[family].curve) == "function",
    "TS " .. family .. " curve generator must exist")
  assert(type(ts.families[family].plots) == "function",
    "TS " .. family .. " plot serializer must exist")
end

local hs = lcp.diagram.get_type("HS")
assert(type(hs) == "table", "HS diagram registry entry must exist")
for _, family in ipairs({"quality", "isochore", "isotherm", "isobar"}) do
  assert(type(hs.families[family]) == "table",
    "HS " .. family .. " family must be registered")
  assert(type(hs.families[family].curve) == "function",
    "HS " .. family .. " curve generator must exist")
  assert(type(hs.families[family].plots) == "function",
    "HS " .. family .. " plot serializer must exist")
end

local pt = lcp.diagram.get_type("PT")
assert(type(pt) == "table", "PT diagram registry entry must exist")
for _, family in ipairs({"phase_envelope", "isentrope", "isenthalp",
    "isochore"}) do
  assert(type(pt.families[family]) == "table",
    "PT " .. family .. " family must be registered")
  assert(type(pt.families[family].curve) == "function",
    "PT " .. family .. " curve generator must exist")
  assert(type(pt.families[family].plots) == "function",
    "PT " .. family .. " plot serializer must exist")
end

local curve_options = {
  fluid = "R134a",
  pressure_min = 1e5,
  pressure_max = 2e6,
  enthalpy_scale = 1e-3,
  pressure_scale = 1e-5,
}

local curve_cases = {
  {name = "quality", points = ph.quality_curve(curve_options, 0.5)},
  {name = "isotherm", points = ph.isotherm_curve(curve_options, 293.15)},
  {name = "isentrope", points = ph.isentrope_curve(curve_options, 1800)},
  {name = "isochore", points = ph.isochore_curve(curve_options, 0.05)},
}
for _, case in ipairs(curve_cases) do
  assert(#case.points >= 2, case.name .. " curve must contain at least two points")
  for index, point in ipairs(case.points) do
    assert(finite(point.x) and finite(point.y),
      string.format("%s point %d must have finite plot coordinates",
        case.name, index))
    assert(point.y > 0, case.name .. " pressure coordinate must be positive")
  end
end

local pv_curve_options = {
  fluid = "R134a",
  pressure_min = 1e5,
  pressure_max = 2e6,
  specific_volume_scale = 1,
  pressure_scale = 1e-5,
}
local pv_curve_cases = {
  {name = "quality", points = pv.quality_curve(pv_curve_options, 0.5)},
  {name = "isotherm", points = pv.isotherm_curve(pv_curve_options, 293.15)},
  {name = "isentrope", points = pv.isentrope_curve(pv_curve_options, 1800)},
}
for _, case in ipairs(pv_curve_cases) do
  assert(#case.points >= 2,
    "PV " .. case.name .. " curve must contain at least two points")
  for index, point in ipairs(case.points) do
    assert(finite(point.x) and point.x > 0
        and finite(point.y) and point.y > 0,
      string.format("PV %s point %d must be valid on log-log axes",
        case.name, index))
    assert_close(point.x, point.v, 2e-13,
      "PV x coordinate must equal SI specific volume at scale one")
    assert_close(point.y, point.p * 1e-5, 2e-13,
      "PV y coordinate must equal scaled pressure")
  end
end

local ts_curve_options = {
  fluid = "R134a",
  pressure_min = 1e5,
  pressure_max = 2e6,
  entropy_scale = 1e-3,
  temperature_scale = 1,
}
local ts_curve_cases = {
  {name = "quality", points = ts.quality_curve(ts_curve_options, 0.5)},
  {name = "isenthalp", points = ts.isenthalp_curve(ts_curve_options, 3e5)},
}
for _, case in ipairs(ts_curve_cases) do
  assert(#case.points >= 2,
    "TS " .. case.name .. " curve must contain at least two points")
  for index, point in ipairs(case.points) do
    assert(finite(point.x) and finite(point.y) and point.y > 0,
      string.format("TS %s point %d must have finite coordinates",
        case.name, index))
    assert_close(point.x, point.s * 1e-3, 2e-13,
      "TS x coordinate must equal scaled specific entropy")
    assert_close(point.y, point.T, 2e-13,
      "TS y coordinate must equal absolute temperature")
  end
end

local hs_curve_options = {
  fluid = "Water", pressure_min = 1e4, pressure_max = 3e7,
  entropy_scale = 1e-3, enthalpy_scale = 1e-3,
  temperature_unit = "kelvin",
  isobar_temperature_min = 280, isobar_temperature_max = 700,
}
local hs_curve_cases = {
  {name = "quality", points = hs.quality_curve(hs_curve_options, 0.5)},
  {name = "isochore", points = hs.isochore_curve(hs_curve_options, 0.1)},
  {name = "isotherm", points = hs.isotherm_curve(hs_curve_options, 473.15)},
  {name = "isobar", points = hs.isobar_curve(hs_curve_options, 1e5)},
}
for _, case in ipairs(hs_curve_cases) do
  assert(#case.points >= 2,
    "HS " .. case.name .. " curve must contain at least two points")
  for _, point in ipairs(case.points) do
    assert_close(point.x, point.s * 1e-3, 2e-13,
      "HS x coordinate must equal scaled entropy")
    assert_close(point.y, point.h * 1e-3, 2e-13,
      "HS y coordinate must equal scaled enthalpy")
  end
end

local pt_options = {
  fluid = "R134a",
  pressure_min = constants.ptriple,
  pressure_max = constants.pcrit * 1.1,
  temperature_axis_unit = "celsius",
  pressure_axis_unit = "bar",
}
local phase_envelope = pt.phase_envelope_curve(pt_options)
assert(#phase_envelope >= 3,
  "PT phase envelope must contain sampled points and exact endpoints")
assert(phase_envelope[1].triple and phase_envelope[#phase_envelope].critical,
  "PT phase envelope must expose exact triple and critical endpoints")
assert_close(phase_envelope[1].T, constants.Ttriple, 2e-12,
  "PT phase envelope triple-point temperature")
assert_close(phase_envelope[1].p, constants.ptriple, 2e-12,
  "PT phase envelope triple-point pressure")
assert_close(phase_envelope[#phase_envelope].T, constants.Tcrit, 2e-12,
  "PT phase envelope critical temperature")
assert_close(phase_envelope[#phase_envelope].p, constants.pcrit, 2e-12,
  "PT phase envelope critical pressure")
assert(phase_envelope[1].quality_defined == false
    and phase_envelope[#phase_envelope].quality_defined == false,
  "quality must remain undefined on the single PT coexistence locus")
for index = 2, #phase_envelope do
  assert(phase_envelope[index].p > phase_envelope[index - 1].p,
    "PT coexistence pressure must increase strictly from triple to critical")
  assert(phase_envelope[index].T >= phase_envelope[index - 1].T - 1e-9,
    "PT coexistence temperature must not decrease")
end

-- Endpoint and monotonicity invariants must follow the selected fluid rather
-- than an R134a-specific table.  These two fluids exercise markedly different
-- triple pressures and critical temperatures.
for _, fluid in ipairs({"Water", "Propane"}) do
  local fluid_constants = lcp.fluid_constants(fluid)
  local curve = pt.phase_envelope_curve({
    fluid = fluid,
    pressure_min = fluid_constants.ptriple,
    pressure_max = fluid_constants.pcrit * 1.02,
    initial_intervals = 6,
    max_depth = 4,
  })
  assert(#curve >= 3 and curve[1].triple and curve[#curve].critical,
    "PT coexistence endpoints must be available for " .. fluid)
  assert_close(curve[1].p, fluid_constants.ptriple, 2e-12,
    fluid .. " PT triple pressure")
  assert_close(curve[#curve].p, fluid_constants.pcrit, 2e-12,
    fluid .. " PT critical pressure")
  for index = 2, #curve do
    assert(curve[index].p > curve[index - 1].p,
      fluid .. " PT coexistence pressure ordering")
    assert(curve[index].T >= curve[index - 1].T - 1e-9,
      fluid .. " PT coexistence temperature ordering")
  end
end

local pt_cases = {
  {name = "isentrope", points = pt.isentrope_curve(pt_options, 1700)},
  {name = "isenthalp", points = pt.isenthalp_curve(pt_options, 3e5)},
  {name = "isochore", points = pt.isochore_curve(pt_options, 0.02)},
}
for _, case in ipairs(pt_cases) do
  assert(#case.points >= 2, "PT " .. case.name .. " must contain points")
  for _, point in ipairs(case.points) do
    assert_close(point.x, point.T - 273.15, 2e-12,
      "PT Celsius coordinate coupling")
    assert_close(point.y, point.p * 1e-5, 2e-12,
      "PT bar coordinate coupling")
  end
  assert(type(case.points.segments) == "table",
    "PT adaptive curves must expose connected-domain metadata")
end

local pt_process, pt_process_metadata = pt.process_points({
  fluid = "R134a", temperature_axis_unit = "celsius",
  pressure_axis_unit = "bar", type = "isentropic",
  from = "pressure=2bar,quality=1", to = "pressure=10bar",
})
assert(#pt_process >= 2 and pt_process_metadata.diagram_type == "PT",
  "PT process must be sampled through the shared process contract")
assert_close(pt_process_metadata.from.x,
  pt_process_metadata.from.T - 273.15, 2e-12,
  "PT process Celsius coordinate")

local plateau = pv.isotherm_curve({
  fluid = "R134a", pressure_min = 1e5, pressure_max = 5e6,
}, 293.15)
local has_two_phase_plateau = false
for index = 2, #plateau do
  if plateau[index - 1].q == 1 and plateau[index].q == 0
      and plateau[index - 1].p == plateau[index].p
      and plateau[index - 1].v > plateau[index].v then
    has_two_phase_plateau = true
  end
end
assert(has_two_phase_plateau,
  "subcritical PV isotherms must contain their exact two-phase plateau")

-- Saturated liquid and saturated vapour are distinct at the triple pressure,
-- but every quality curve must converge to one common state at the critical
-- point.  Exercise several fluid families so this invariant cannot regress to
-- a refrigerant-specific workaround.
for _, fluid in ipairs({
    "R134a", "R717", "Water", "CO2", "Propane", "Methane", "Nitrogen",
  }) do
  local fluid_constants = lcp.fluid_constants(fluid)
  local critical_options = {
    fluid = fluid,
    pressure_min = math.max(fluid_constants.ptriple * 1.05,
      fluid_constants.pcrit * 0.02),
    pressure_max = fluid_constants.pcrit * 1.05,
    enthalpy_scale = 1e-3,
    pressure_scale = 1e-5,
    initial_intervals = 6,
    max_depth = 5,
  }
  local q0 = ph.quality_curve(critical_options, 0)
  local qhalf = ph.quality_curve(critical_options, 0.5)
  local q1 = ph.quality_curve(critical_options, 1)
  assert(#q0 >= 5 and #qhalf >= 5 and #q1 >= 5,
    fluid .. " quality curves must be refined near the critical point")
  local ends = {q0[#q0], qhalf[#qhalf], q1[#q1]}
  for _, endpoint in ipairs(ends) do
    assert(endpoint.critical == true,
      fluid .. " quality curve must expose its canonical critical endpoint")
    assert(endpoint.critical_limit == true
        and endpoint.quality_defined == false and endpoint.q == nil,
      fluid .. " critical endpoint must be a quality-independent limit")
    assert_close(endpoint.p, fluid_constants.pcrit, 2e-13,
      fluid .. " quality endpoint pressure")
  end
  assert(ends[1].x == ends[2].x and ends[2].x == ends[3].x
      and ends[1].y == ends[2].y and ends[2].y == ends[3].y,
    fluid .. " quality curves must have one bit-identical critical endpoint")
  for _, points in ipairs({q0, qhalf, q1}) do
    local penultimate = points[#points - 1]
    assert(penultimate.p < fluid_constants.pcrit
        and penultimate.p > 0.989 * fluid_constants.pcrit,
      fluid .. " quality sampling must approach the critical point from below")
  end

  local pv_q0 = pv.quality_curve(critical_options, 0)
  local pv_qhalf = pv.quality_curve(critical_options, 0.5)
  local pv_q1 = pv.quality_curve(critical_options, 1)
  local pv_ends = {pv_q0[#pv_q0], pv_qhalf[#pv_qhalf], pv_q1[#pv_q1]}
  assert(pv_ends[1].critical and pv_ends[2].critical
      and pv_ends[3].critical,
    fluid .. " PV qualities must reach the canonical critical state")
  for _, endpoint in ipairs(pv_ends) do
    assert(endpoint.critical_limit == true
        and endpoint.quality_defined == false and endpoint.q == nil,
      fluid .. " PV critical endpoint must not carry a quality")
  end
  assert(pv_ends[1].x == pv_ends[2].x
      and pv_ends[2].x == pv_ends[3].x
      and pv_ends[1].y == pv_ends[2].y
      and pv_ends[2].y == pv_ends[3].y,
    fluid .. " PV quality curves must share one critical endpoint")
  assert_close(pv_ends[1].v, 1 / fluid_constants.rhocrit, 2e-13,
    fluid .. " PV critical volume")

  local ts_q0 = ts.quality_curve(critical_options, 0)
  local ts_qhalf = ts.quality_curve(critical_options, 0.5)
  local ts_q1 = ts.quality_curve(critical_options, 1)
  local ts_ends = {ts_q0[#ts_q0], ts_qhalf[#ts_qhalf], ts_q1[#ts_q1]}
  assert(ts_ends[1].critical and ts_ends[2].critical
      and ts_ends[3].critical,
    fluid .. " TS qualities must reach the canonical critical state")
  for _, endpoint in ipairs(ts_ends) do
    assert(endpoint.critical_limit == true
        and endpoint.quality_defined == false and endpoint.q == nil,
      fluid .. " TS critical endpoint must not carry a quality")
  end
  assert(ts_ends[1].x == ts_ends[2].x
      and ts_ends[2].x == ts_ends[3].x
      and ts_ends[1].y == ts_ends[2].y
      and ts_ends[2].y == ts_ends[3].y,
    fluid .. " TS quality curves must share one critical endpoint")
  assert_close(ts_ends[1].T, fluid_constants.Tcrit, 2e-13,
    fluid .. " TS critical temperature")

  local hs_q0 = hs.quality_curve(critical_options, 0)
  local hs_qhalf = hs.quality_curve(critical_options, 0.5)
  local hs_q1 = hs.quality_curve(critical_options, 1)
  local hs_ends = {hs_q0[#hs_q0], hs_qhalf[#hs_qhalf], hs_q1[#hs_q1]}
  assert(hs_ends[1].critical and hs_ends[2].critical
      and hs_ends[3].critical,
    fluid .. " HS qualities must reach the canonical critical state")
  for _, endpoint in ipairs(hs_ends) do
    assert(endpoint.critical_limit == true
        and endpoint.quality_defined == false and endpoint.q == nil,
      fluid .. " HS critical endpoint must not carry a quality")
  end
  assert(hs_ends[1].x == hs_ends[2].x
      and hs_ends[2].x == hs_ends[3].x
      and hs_ends[1].y == hs_ends[2].y
      and hs_ends[2].y == hs_ends[3].y,
    fluid .. " HS quality curves must share one critical endpoint")

  local strict_tolerance = 0.002
  local strict_options = {
    fluid = fluid,
    pressure_min = critical_options.pressure_min,
    pressure_max = critical_options.pressure_max,
    entropy_scale = 1e-3,
    enthalpy_scale = 1e-3,
    initial_intervals = 4,
    max_depth = 12,
    tolerance = strict_tolerance,
  }
  for _, quality in ipairs({0, 0.5, 1}) do
    for _, case in ipairs({
        {name = "TS", points = ts.quality_curve(strict_options, quality)},
        {name = "HS", points = hs.quality_curve(strict_options, quality)},
      }) do
      local checked = 0
      for index = 1, #case.points - 1 do
        local first, second = case.points[index], case.points[index + 1]
        local epsilon_first = 1 - first.p / fluid_constants.pcrit
        if epsilon_first > 0 and epsilon_first <= 0.5
            and not second.critical then
          local defect = critical_quality_midpoint_error(case.name, fluid,
            quality, fluid_constants.pcrit, first, second)
          assert(defect <= strict_tolerance * (1 + 1e-8),
            string.format(
              "%s %s Q=%.1f critical-tail defect %.6g exceeds %.6g",
              fluid, case.name, quality, defect, strict_tolerance))
          checked = checked + 1
        end
      end
      assert(checked > 0,
        fluid .. " " .. case.name
          .. " quality curve must exercise adaptive critical-tail segments")
    end
  end
end

-- CoolProp's built-in fluid list changes between releases.  Audit it
-- dynamically so updating the external library also tests every newly added
-- fluid without hard-coding package data in LuaCoolProp.
local audited_fluid_count = 0
for fluid in lcp.global_param_string("FluidsList"):gmatch("[^,]+") do
  if lcp.fluid_param_string(fluid, "pure"):lower() == "true" then
    local fluid_constants = lcp.fluid_constants(fluid)
  assert(finite(fluid_constants.pcrit) and finite(fluid_constants.ptriple)
      and fluid_constants.pcrit > fluid_constants.ptriple,
    fluid .. " must expose an ordered triple-to-critical pressure interval")
  local audit_options = {
    fluid = fluid,
    pressure_min = math.max(fluid_constants.ptriple * 1.05,
      fluid_constants.pcrit * 0.02),
    pressure_max = fluid_constants.pcrit * 1.01,
    initial_intervals = 2,
    max_depth = 1,
  }
  local liquid = ph.quality_curve(audit_options, 0)
  local vapour = ph.quality_curve(audit_options, 1)
  local liquid_end, vapour_end = liquid[#liquid], vapour[#vapour]
  assert(liquid_end and vapour_end
      and liquid_end.critical == true and vapour_end.critical == true
      and liquid_end.p == vapour_end.p
      and liquid_end.x == vapour_end.x and liquid_end.y == vapour_end.y,
    fluid .. " saturation branches must share one critical endpoint")
  local pv_liquid = pv.quality_curve(audit_options, 0)
  local pv_vapour = pv.quality_curve(audit_options, 1)
  local pv_liquid_end = pv_liquid[#pv_liquid]
  local pv_vapour_end = pv_vapour[#pv_vapour]
  assert(pv_liquid_end and pv_vapour_end
      and pv_liquid_end.critical == true
      and pv_vapour_end.critical == true
      and pv_liquid_end.x == pv_vapour_end.x
      and pv_liquid_end.y == pv_vapour_end.y,
    fluid .. " PV saturation branches must share one critical endpoint")
  local ts_liquid = ts.quality_curve(audit_options, 0)
  local ts_vapour = ts.quality_curve(audit_options, 1)
  local ts_liquid_end = ts_liquid[#ts_liquid]
  local ts_vapour_end = ts_vapour[#ts_vapour]
  assert(ts_liquid_end and ts_vapour_end
      and ts_liquid_end.critical == true
      and ts_vapour_end.critical == true
      and ts_liquid_end.x == ts_vapour_end.x
      and ts_liquid_end.y == ts_vapour_end.y,
    fluid .. " TS saturation branches must share one critical endpoint")
  local hs_liquid = hs.quality_curve(audit_options, 0)
  local hs_vapour = hs.quality_curve(audit_options, 1)
  local hs_liquid_end = hs_liquid[#hs_liquid]
  local hs_vapour_end = hs_vapour[#hs_vapour]
  assert(hs_liquid_end and hs_vapour_end
      and hs_liquid_end.critical and hs_vapour_end.critical
      and hs_liquid_end.x == hs_vapour_end.x
      and hs_liquid_end.y == hs_vapour_end.y,
    fluid .. " HS saturation branches must share one critical endpoint")
    audited_fluid_count = audited_fluid_count + 1
  end
end
assert(audited_fluid_count > 0,
  "CoolProp must report at least one built-in fluid for the closure audit")

local truncated_quality = ph.quality_curve({
  fluid = "R134a",
  pressure_min = constants.ptriple * 1.05,
  pressure_max = constants.pcrit * 0.8,
}, 0)
assert(#truncated_quality >= 2
    and truncated_quality[#truncated_quality].critical ~= true,
  "a quality curve requested below pc must not invent a critical endpoint")
assert_close(truncated_quality[#truncated_quality].p,
  constants.pcrit * 0.8, 2e-13,
  "a truncated quality curve must preserve its requested upper pressure")

local plot_code = ph.plots({
  fluid = "R134a",
  quality = true,
  isotherm = true,
  quality_values = "0,0.5,1",
  temperature_values = "-20,0,20",
  temperature_unit = "celsius",
  labels = true,
})
assert(plot_code:find("\\addplot", 1, true),
  "PH serializer must emit PGFPlots plots")
assert(plot_code:find("name path=lcp-ph-q-0p5", 1, true),
  "PH serializer must emit stable quality-curve names")
assert(plot_code:find("\\pgfplotsautonode", 1, true),
  "labelled PH serialization must delegate to pgfplots-autonode")
assert(plot_code:find("\\csname LCP@format@number\\endcsname", 1, true),
  "dimensionless generated labels must use the TeX number formatter")
assert(plot_code:find("\\csname LCP@format@quantity\\endcsname", 1, true)
    and plot_code:find("\\degreeCelsius", 1, true),
  "unit-bearing generated labels must use the TeX quantity formatter")

local symbol_code = ph.plots({
  fluid = "R134a",
  quality = true,
  isotherm = true,
  isentrope = true,
  isochore = true,
  quality_values = "0.5",
  temperature_values = "20",
  entropy_values = "1.8",
  specific_volume_values = "0.05",
  quality_symbol = "x",
  temperature_symbol = "\\theta",
  entropy_symbol = "\\sigma",
  specific_volume_symbol = "\\nu",
  labels = true,
  label_every = 1,
})
for _, labelled_symbol in ipairs({"x", "\\theta", "\\sigma", "\\nu"}) do
  assert(symbol_code:find("$" .. labelled_symbol .. "=", 1, true),
    labelled_symbol .. " must replace the corresponding generated symbol")
end

local pv_plot_code = pv.plots({
  fluid = "R134a",
  quality = true,
  isotherm = true,
  isentrope = true,
  quality_values = "0,0.5,1",
  temperature_values = "0,20",
  entropy_values = "1.7,1.9",
  quality_symbol = "x",
  temperature_symbol = "\\theta",
  entropy_symbol = "\\sigma",
  labels = true,
  label_every = 1,
})
assert(pv_plot_code:find("name path=lcp-pv-q-0p5", 1, true),
  "PV serializer must emit diagram-specific stable path names")
for _, labelled_symbol in ipairs({"x", "\\theta", "\\sigma"}) do
  assert(pv_plot_code:find("$" .. labelled_symbol .. "=", 1, true),
    "PV labels must honour the configured " .. labelled_symbol .. " symbol")
end

local ts_plot_code = ts.plots({
  fluid = "R134a",
  quality = true,
  isenthalp = true,
  quality_values = "0,0.5,1",
  enthalpy_values = "200,300,400",
  enthalpy_unit = "kjkg",
  quality_symbol = "x",
  enthalpy_symbol = "\\eta",
  labels = true,
  label_every = 1,
})
assert(ts_plot_code:find("name path=lcp-ts-q-0p5", 1, true),
  "TS serializer must emit diagram-specific stable path names")
assert(ts_plot_code:find("name path=lcp-ts-h-300", 1, true),
  "TS serializer must emit stable isenthalp path names")
for _, labelled_symbol in ipairs({"x", "\\eta"}) do
  assert(ts_plot_code:find("$" .. labelled_symbol .. "=", 1, true),
    "TS labels must honour the configured " .. labelled_symbol .. " symbol")
end

local hs_plot_code = hs.plots({
  fluid = "Water", pressure_min = 1e4, pressure_max = 3e7,
  quality = true, isochore = true, isotherm = true, isobar = true,
  quality_values = "0,0.5,1", specific_volume_values = "0.01,0.1",
  temperature_values = "100,200,300", isobar_values = "1,10,100",
  labels = true, label_every = 1,
})
for _, name in ipairs({"lcp-hs-q-0p5", "lcp-hs-v-", "lcp-hs-T-",
    "lcp-hs-p-"}) do
  assert(hs_plot_code:find("name path=" .. name, 1, true),
    "HS serializer must emit stable " .. name .. " paths")
end

local process_points, process_metadata = ph.process_points({
  fluid = "R134a",
  type = "isentrope",
  from = "pressure=2bar,quality=1",
  to = "pressure=10bar",
})
assert(#process_points >= 2, "isentrope process must contain points")
assert(process_metadata.kind == "isentropic",
  "process vocabulary must normalize to the thermodynamic implementation")
assert_close(process_metadata.from.s, process_metadata.to.s, 2e-8,
  "isentrope endpoints must conserve entropy")
for _, endpoint in ipairs({process_metadata.from, process_metadata.to}) do
  assert(finite(endpoint.x) and finite(endpoint.y),
    "completed process endpoints must expose finite plot coordinates")
end

local process_code, rendered_process_metadata = ph.process_plot({
  fluid = "R134a",
  type = "isobar",
  pressure = "8bar",
  from = "quality=0",
  to = "quality=1",
  name = "test-process-path",
  coord_digits = 7,
})
assert(rendered_process_metadata.id == "test-process-path",
  "a process name must become the metadata path ID")
assert(process_code:find("name path=test-process-path", 1, true),
  "a named process must expose a TikZ path for intersections")
assert(process_code:find("(" .. tostring(rendered_process_metadata.from.x):sub(1, 5),
    1, true),
  "process serialization must contain its numeric endpoint")

local _, throttle_metadata = ph.process_points({
  fluid = "R134a",
  type = "isenthalp",
  from = "pressure=18bar,temperature=60C",
  to = "pressure=3bar",
})
assert(finite(throttle_metadata.to.q)
    and throttle_metadata.to.q > 0 and throttle_metadata.to.q < 1,
  "a two-phase process endpoint must expose its resolved quality")
assert(finite(throttle_metadata.to.s) and finite(throttle_metadata.to.v),
  "a resolved two-phase endpoint must expose entropy and specific volume")

local pv_process_points, pv_process_metadata = pv.process_points({
  fluid = "R134a",
  type = "isentropic",
  from = "pressure=2bar,quality=1",
  to = "pressure=12bar",
})
assert(#pv_process_points >= 2,
  "PV isentropic process must contain sampled points")
assert(pv_process_metadata.id:find("lcp%-pv%-process"),
  "PV processes must use diagram-specific path identifiers")
assert_close(pv_process_metadata.from.s, pv_process_metadata.to.s, 2e-8,
  "PV process endpoints must conserve entropy")
for _, endpoint in ipairs({pv_process_metadata.from, pv_process_metadata.to}) do
  assert_close(endpoint.x, endpoint.v, 2e-13,
    "PV process endpoint x must be its specific volume")
  assert(endpoint.x > 0 and endpoint.y > 0,
    "PV process endpoints must be valid on log-log axes")
end

local pv_process_code = pv.process_plot({
  fluid = "R134a", type = "isobar", pressure = "8bar",
  from = "quality=0", to = "quality=1", name = "test-pv-process",
})
assert(pv_process_code:find("name path=test-pv-process", 1, true),
  "named PV processes must expose TikZ paths for intersections")

local ts_process_points, ts_process_metadata = ts.process_points({
  fluid = "R134a",
  type = "isentropic",
  from = "pressure=2bar,quality=1",
  to = "pressure=12bar",
})
assert(#ts_process_points == 2,
  "a TS isentrope must be represented by its exact straight endpoints")
assert_close(ts_process_metadata.from.s, ts_process_metadata.to.s, 2e-8,
  "TS isentrope endpoints must conserve entropy")
assert_close(ts_process_metadata.from.x, ts_process_metadata.to.x, 2e-8,
  "TS isentrope endpoints must share one horizontal coordinate")
for _, endpoint in ipairs({ts_process_metadata.from, ts_process_metadata.to}) do
  assert_close(endpoint.x, endpoint.s * 1e-3, 2e-13,
    "TS process endpoint x must be scaled entropy")
  assert_close(endpoint.y, endpoint.T, 2e-13,
    "TS process endpoint y must be absolute temperature")
end

local ts_isobar = ts.process_points({
  fluid = "R134a", type = "isobar", pressure = "8bar",
  from = "quality=0", to = "quality=1",
})
assert(#ts_isobar >= 2,
  "a TS two-phase isobar must contain its entropy-parametrized plateau")
for _, point in ipairs(ts_isobar) do
  assert_close(point.T, ts_isobar[1].T, 2e-8,
    "TS two-phase isobar temperature")
end

local hs_process, hs_metadata = hs.process_points({
  fluid = "Water", type = "isentropic",
  from = "pressure=1bar,quality=1", to = "pressure=20bar",
})
assert(#hs_process == 2,
  "an HS isentrope must use its exact vertical endpoints")
assert_close(hs_metadata.from.x, hs_metadata.to.x, 2e-8,
  "HS isentrope plot coordinate")
for _, endpoint in ipairs({hs_metadata.from, hs_metadata.to}) do
  assert_close(endpoint.x, endpoint.s * 1e-3, 2e-13,
    "HS process x coordinate")
  assert_close(endpoint.y, endpoint.h * 1e-3, 2e-13,
    "HS process y coordinate")
end

-- A rigid two-phase tank is a demanding isochore case: its initial T-Q pair
-- determines v, while the target gives T only.  Completing pressure first must
-- not replace the valid target T-rho pair with the ambiguous saturation P-T
-- pair before enthalpy, entropy, and quality are evaluated.
local butane_isochore, butane_isochore_metadata = pv.process_points({
  fluid = "Butane", type = "isochore",
  from = "temperature=305K,quality=0.3",
  to = "temperature=335K",
})
assert(#butane_isochore == 2,
  "an isochore must be emitted as one exact vertical segment")
assert_close(butane_isochore_metadata.from.v,
  butane_isochore_metadata.to.v, 2e-13,
  "PV isochore endpoints must retain one specific volume")
assert_close(butane_isochore_metadata.from.q, 0.3, 2e-11,
  "butane isochore initial quality")
assert_close(butane_isochore_metadata.to.q, 0.665250479, 2e-8,
  "butane isochore target quality resolved from T and v")
assert(butane_isochore_metadata.to.p > butane_isochore_metadata.from.p,
  "heating the rigid two-phase butane tank must raise its pressure")

-- TeX-facing quantities use an intentionally strict, locale-independent
-- grammar.  Decimal floating-point and E notation are accepted, while a
-- decimal comma, a partial token, or an unknown unit is always an error.
local _, scientific_pressure = ph.process_points({
  fluid = "R134a", type = "isobar", pressure = "1.2E5Pa",
  from = "quality=0", to = "quality=1",
})
assert_close(scientific_pressure.constant, 1.2e5, 1e-13,
  "scientific-notation pressure with an attached SI unit")
assert_close(scientific_pressure.from.p, 1.2e5, 1e-12,
  "scientific-notation pressure propagated to a process endpoint")

local _, decimal_quality = ph.process_points({
  fluid = "R134a", type = "isoquality", quality = "5.0E-1",
  from = "pressure=1.2E5Pa", to = "pressure=8.0E5Pa",
})
assert_close(decimal_quality.constant, 0.5, 1e-13,
  "floating-point quality in scientific notation")

local decimal_bound_curve = ph.quality_curve({
  fluid = "R134a", pressure_min = "1.0E5", pressure_max = "2.0E6",
}, 0.5)
assert(#decimal_bound_curve >= 2,
  "diagram APIs must accept strict decimal scientific bounds")

assert_error_contains(function()
  ph.quality_curve({
    fluid = "R134a", pressure_min = "0x186A0", pressure_max = "2E6",
  }, 0.5)
end, "full finite decimal number", "hexadecimal Lua numeric spelling")

assert_error_contains(function()
  ph.process_points({
    fluid = "R134a", type = "isobar", pressure = "1,2bar",
    from = "quality=0", to = "quality=1",
  })
end, "use '.' as the decimal separator", "decimal-comma quantity")

assert_error_contains(function()
  ph.process_points({
    fluid = "R134a", type = "isobar", pressure = "1.2E5Paa",
    from = "quality=0", to = "quality=1",
  })
end, "unknown pressure unit", "unknown quantity unit")

assert_error_contains(function()
  ph.axis_style({
    fluid = "R134a", pressure_min = "not-a-number",
    pressure_max = 2e6,
  })
end, "invalid numeric value", "malformed numeric key")

assert_error_contains(function()
  ph.axis_style({fluid = "R134a", pressure_min = 2e6,
    pressure_max = 1e6})
end, "pressure_max must be greater", "reversed pressure bounds")

assert_error_contains(function()
  ph.plots({
    fluid = "R134a", pressure_min = 1e5, pressure_max = 2e5,
    isoquality = "treu",
  })
end, "invalid boolean value", "malformed boolean key")

assert_error_contains(function()
  ph.plots({
    fluid = "R134a", pressure_min = 1e5, pressure_max = 2e5,
    quality_values = "0,broken,1",
  })
end, "invalid numeric list element", "malformed numeric-list member")

assert_error_contains(function()
  ph.plots({
    fluid = "R134a", pressure_min = 1e5, pressure_max = 2e5,
    quality_mode = "linar",
  })
end, "invalid grid mode", "unknown grid mode")

assert_error_contains(function()
  ph.plots({
    fluid = "R134a", pressure_min = 1e5, pressure_max = 2e5,
    quality_preset = "dens",
  })
end, "invalid quality preset", "unknown grid preset")

-- A redundant property must describe the state selected by the defining
-- CoolProp pair; it is never retained merely because the user supplied it.
assert_error_contains(function()
  ph.process_points({
    fluid = "R134a", type = "isobar",
    from = "pressure=1bar,quality=0,temperature=500K",
    to = "pressure=1bar,quality=1",
  })
end, "inconsistent temperature T", "overdetermined contradictory state")

-- Likewise, both endpoints must satisfy the declared conservation law even
-- when each endpoint is a valid thermodynamic state on its own.
assert_error_contains(function()
  ph.process_points({
    fluid = "R134a", type = "isobar",
    from = "pressure=1bar,quality=0",
    to = "pressure=2bar,quality=1",
  })
end, "inconsistent pressure p", "contradictory process endpoints")

local saturation_temperature = lcp.propsSI(
  "T", "P", 1e5, "Q", 0, "R134a")
local _, consistent_redundancy = ph.process_points({
  fluid = "R134a", type = "isobar", pressure = "1E5Pa",
  from = string.format("p=1E5Pa,Q=0,T=%.12gK", saturation_temperature),
  to = "p=1E5Pa,Q=1",
})
assert_close(consistent_redundancy.from.T, saturation_temperature, 5e-7,
  "consistent redundant state property")

local unknown_ok, unknown_error = pcall(lcp.diagram.get_type, "UNKNOWN")
assert(not unknown_ok and tostring(unknown_error):find("unknown diagram type", 1, true),
  "unknown diagram types must fail with a useful error")

local dummy = {plots = function() return "dummy" end}
lcp.diagram.register_type("TEST", dummy)
assert(lcp.diagram.get_type("TEST") == dummy,
  "diagram registration must preserve uppercase identifiers")

-- Background families are capabilities, not keys inherited from whichever
-- diagram happened to be implemented first.
assert_error_contains(function()
  pv.plots({fluid = "R134a", pressure_min = 1e5, pressure_max = 2e5,
    isochore = true})
end, "does not provide the 'isochore'", "unsupported PV family")
assert_error_contains(function()
  ts.plots({fluid = "R134a", pressure_min = 1e5, pressure_max = 2e5,
    isotherm = true})
end, "does not provide the 'isotherm'", "unsupported TS family")
assert_error_contains(function()
  hs.plots({fluid = "R134a", pressure_min = 1e5, pressure_max = 2e5,
    isentrope = true})
end, "does not provide the 'isentrope'", "unsupported HS family")
assert_error_contains(function()
  ph.plots({fluid = "R134a", pressure_min = 1e5, pressure_max = 2e5,
    isenthalp = true})
end, "does not provide the 'isenthalp'", "unsupported PH family")
assert_error_contains(function()
  pt.plots({fluid = "R134a", quality = true})
end, "does not provide the 'quality'", "unsupported PT quality family")
assert_error_contains(function()
  pv.plots({fluid = "R134a", pressure_min = 1e5, pressure_max = 2e5,
    enthalpy_axis_unit = "kjkg"})
end, "is not an axis option of diagram type PV",
  "coordinate-unit capability mismatch")

-- A semantic axis unit changes the coordinate and its advertised unit as one
-- operation.  Raw SI fields remain invariant.
local unit_base = {fluid = "R134a", pressure_min = 1e5,
  pressure_max = 2e5, initial_intervals = 2, max_depth = 1}
local ph_named_default = ph.quality_curve(unit_base, 0.5)[1]
local ph_named_si = ph.quality_curve({fluid = "R134a", pressure_min = 1e5,
  pressure_max = 2e5, initial_intervals = 2, max_depth = 1,
  enthalpy_axis_unit = "jkg", pressure_axis_unit = "mpa"}, 0.5)[1]
assert_close(ph_named_si.h, ph_named_default.h, 1e-13,
  "semantic enthalpy unit preserves SI data")
assert_close(ph_named_si.x, ph_named_default.x * 1000, 1e-13,
  "J/kg coordinate scaling")
assert_close(ph_named_si.y, ph_named_default.y * 0.1, 1e-13,
  "MPa coordinate scaling")
local ts_kelvin = ts.quality_curve(unit_base, 0.5)[1]
local ts_celsius = ts.quality_curve({fluid = "R134a", pressure_min = 1e5,
  pressure_max = 2e5, initial_intervals = 2, max_depth = 1,
  temperature_axis_unit = "celsius"}, 0.5)[1]
assert_close(ts_celsius.y, ts_kelvin.y - 273.15, 1e-12,
  "Celsius coordinate offset")

-- Diagram samplers return explicit connected components.  The serializer
-- consumes `_break_before` markers and therefore cannot bridge an invalid
-- CoolProp domain with a fictitious line segment.
local segmented = ph.isentrope_curve({fluid = "R134a",
  pressure_min = 1e3, pressure_max = 1e9, initial_intervals = 4,
  max_depth = 4, domain_policy = "ignore"}, 1700)
assert(type(segmented.segments) == "table"
    and segmented.domain_gap_count == math.max(0, #segmented.segments - 1)
    and type(segmented.omitted_sample_count) == "number",
  "adaptive curves expose connected-domain metadata")

-- Diagram construction is deliberately pure-fluid-only.  This restriction
-- does not leak into the low-level CoolProp wrappers.
local mixture = "HEOS::R32[0.697615]&R125[0.302385]"
assert(finite(lcp.propsSI("H", "P", 1e5, "T", 300, mixture)),
  "low-level wrappers must remain mixture-capable")
assert_error_contains(function()
  ph.axis_style({fluid = mixture, pressure_min = 1e5,
    pressure_max = 2e6})
end, "diagrams require a pure fluid", "mixture diagram boundary")
assert_error_contains(function()
  pv.axis_style({fluid = "R407C", pressure_min = 1e5,
    pressure_max = 2e6})
end, "diagrams require a pure fluid", "predefined blend diagram boundary")

-- Enthalpy and entropy coordinates in a document share one immutable
-- CoolProp reference-state convention per fluid.
assert(lcp.set_reference_state("Propane", "DEF") == "DEF")
ph.axis_style({fluid = "Propane", pressure_min = 1e5,
  pressure_max = 2e6, reference_state = "DEF"})
assert(lcp.reference_state("Propane") == "DEF")
assert_error_contains(function()
  lcp.set_reference_state("Propane", "IIR")
end, "is locked to DEF", "locked reference-state contract")
local _, reference_metadata = ph.process_points({fluid = "R134a",
  type = "isobar", pressure = "1bar", from = "quality=0",
  to = "quality=1"})
assert(reference_metadata.reference_state == "DEF",
  "process metadata records the reference-state convention")

-- TS, HS and PT need the same sampled records both to derive automatic axis
-- limits and to serialize plots.  The TeX bridge must consume those records
-- exactly once, without changing a plotted coordinate or an autonode path.
local prepared_specs = {
  {code = "TS", generator = "ts_quality_curve",
    opts = {fluid = "R134a", pressure_min = 1e5, pressure_max = 1e6,
      quality = true, q_mode = "list", q_values = "0.5", labels = true,
      initial_intervals = 4, max_depth = 2}},
  {code = "HS", generator = "hs_quality_curve",
    opts = {fluid = "R134a", pressure_min = 1e5, pressure_max = 1e6,
      quality = true, q_mode = "list", q_values = "0.5", labels = true,
      initial_intervals = 4, max_depth = 2}},
  {code = "PT", generator = "pt_phase_envelope_curve",
    opts = {fluid = "R134a", pressure_min = 1e5, pressure_max = 1e6,
      phase_envelope = true, labels = true,
      initial_intervals = 4, max_depth = 2}},
}
local captured = {}
tex = {sprint = function(value) captured[#captured + 1] = value end}
for _, spec in ipairs(prepared_specs) do
  local implementation = lcp.diagram.get_type(spec.code)
  local expected_style = "\\pgfplotsset{every axis/.append style={"
    .. implementation.axis_style(spec.opts) .. "}}"
  local expected_plots = implementation.plots(spec.opts)
  local original = lcp.diagram[spec.generator]
  local calls = 0
  lcp.diagram[spec.generator] = function(...)
    calls = calls + 1
    return original(...)
  end
  captured = {}
  lcp.tex.print_diagram(spec.code, "axis_style", spec.opts)
  lcp.tex.print_diagram(spec.code, "plots", spec.opts)
  lcp.diagram[spec.generator] = original
  assert(calls == 1,
    spec.code .. " high-level TeX diagram must sample its curve once")
  assert(captured[1] == expected_style,
    spec.code .. " prepared axis bounds must match uncached output")
  assert(captured[2] == expected_plots,
    spec.code .. " prepared plot and autonode path must be byte-identical")
end

-- The prepared records are one-shot and keyed by every input option.  A
-- changed family value after axis preparation must trigger new sampling.
local changed_opts = {fluid = "R134a", pressure_min = 1e5,
  pressure_max = 1e6, quality = true, q_mode = "list", q_values = "0.5",
  initial_intervals = 4, max_depth = 2}
captured = {}
lcp.tex.print_diagram("TS", "axis_style", changed_opts)
changed_opts.q_values = "0.2"
local changed_expected = lcp.diagram.ts.plots(changed_opts)
lcp.tex.print_diagram("TS", "plots", changed_opts)
assert(captured[2] == changed_expected,
  "changed curve options must invalidate prepared axis records")

io.write("LuaCoolProp numerical regression tests passed\n")