#pragma once

#include "FluidSim/Engine.hpp"
#include <algorithm>
#include <cmath>
#include <deque>
#include <limits>
#include <stdexcept>
#include <utility>
#include <vector>

namespace fluidsim::port {

class Volume0D : public Component {
public:
    explicit Volume0D(float volume_m3, float pressure = 101325.0f, float temperature = 300.0f,
                      float gas_constant_value = 287.0f, float gamma_value = 1.4f)
        : gamma(gamma_value), gas_constant(gas_constant_value), volume_(volume_m3),
          air_mass_(pressure * volume_m3 / (gas_constant_value * temperature)),
          energy_(pressure * volume_m3 / (gamma_value - 1.0f)) {
        if (!std::isfinite(volume_m3) || volume_m3 <= 0 || !std::isfinite(pressure) || pressure <= 0 ||
            !std::isfinite(temperature) || temperature <= 0 || !std::isfinite(gas_constant_value) ||
            gas_constant_value <= 0 || !std::isfinite(gamma_value) || gamma_value <= 1)
            throw std::invalid_argument("volume thermodynamic state must be finite and positive, gamma > 1");
    }

    float gamma{1.4f};
    float gas_constant{287.0f};
    float ambient_pressure{101325.0f};
    float ambient_temperature{300.0f};
    float energy_relaxation_rate{0.0f};
    float volume_change_rate{};

    Port &create_port() {
        ports_.emplace_back();
        auto &p = ports_.back();
        p.owner_energy = &energy_;
        p.owner_air_mass = &air_mass_;
        p.owner_exhaust_mass = &exhaust_mass_;
        sync_port(p);
        return p;
    }
    float pressure() const { return (gamma - 1.0f) * energy_ / std::max(volume_, 1e-12f); }
    float density() const { return mass() / std::max(volume_, 1e-12f); }
    float temperature() const { return pressure() / std::max(density() * gas_constant, 1e-9f); }
    float volume() const { return volume_; }
    float mass() const { return air_mass_ + exhaust_mass_; }
    float internal_energy() const { return energy_; }
    float air_fraction() const { return air_mass_ / std::max(mass(), 1e-12f); }
    void set_pressure(float value, float temp = -1.0f) {
        if (!std::isfinite(value) || value <= 0 || !std::isfinite(temp))
            throw std::invalid_argument("invalid volume pressure/temperature");
        const float target_t = temp > 0.0f ? temp : temperature();
        const float total = value * volume_ / (gas_constant * target_t);
        const float fraction = air_fraction();
        air_mass_ = total * fraction;
        exhaust_mass_ = total * (1.0f - fraction);
        energy_ = value * volume_ / (gamma - 1.0f);
        air_roundoff_ = exhaust_roundoff_ = energy_roundoff_ = 0;
        sync_ports();
    }
    void update(float dt) {
        if (!std::isfinite(dt) || dt < 0)
            throw std::invalid_argument("invalid volume timestep");
        if (dt == 0)
            return;
        double d_air = 0, d_exhaust = 0, d_energy = 0;
        const float fraction = air_fraction();
        for (const Port &p : ports_) {
            const float flow_fraction = p.mass_flow >= 0.0f ? p.air_fraction : fraction;
            d_air += double(p.mass_flow) * flow_fraction;
            d_exhaust += double(p.mass_flow) * (1.0f - flow_fraction);
            d_energy += double(p.mass_flow) * p.enthalpy;
        }
        accumulate(air_mass_, air_roundoff_, d_air * dt);
        accumulate(exhaust_mass_, exhaust_roundoff_, d_exhaust * dt);
        accumulate(energy_, energy_roundoff_, d_energy * dt);
        const float new_volume = volume_ + volume_change_rate * dt;
        if (!(new_volume > 0) || !std::isfinite(new_volume))
            throw std::runtime_error("invalid moving volume");
        if (new_volume != volume_) {
            const double scale = std::pow(double(volume_) / new_volume, gamma - 1.0f);
            const double precise = (double(energy_)+energy_roundoff_)*scale;
            energy_ = float(precise);
            energy_roundoff_ = precise-energy_;
        }
        volume_ = new_volume;
        if (energy_relaxation_rate > 0.0f && mass() > 1e-12f) {
            const float target = mass() * gas_constant * ambient_temperature / (gamma - 1.0f);
            energy_ = target + (energy_ - target) * std::exp(-energy_relaxation_rate * dt);
            energy_roundoff_ = 0;
        }
        if (!(mass() > 0) || !(energy_ > 0) || !std::isfinite(mass()) || !std::isfinite(energy_))
            throw std::runtime_error(
                "nonphysical volume state; reduce timestep or check external port fluxes");
        if (air_mass_ < -1e-12f || exhaust_mass_ < -1e-12f)
            throw std::runtime_error("negative volume species mass");
        air_mass_ = std::max(air_mass_, 0.0f);
        exhaust_mass_ = std::max(exhaust_mass_, 0.0f);
        sync_ports();
    }
    void update_state(float dt) override { update(dt); }
    std::vector<Port *> ports() override {
        std::vector<Port *> result;
        result.reserve(ports_.size());
        for (Port &port : ports_)
            result.push_back(&port);
        return result;
    }
    const std::deque<Port> &ports() const { return ports_; }

private:
    float volume_, air_mass_{}, exhaust_mass_{}, energy_;
    // Preserve sub-ulp transfers without changing the float port/budget ABI.
    double air_roundoff_{}, exhaust_roundoff_{}, energy_roundoff_{};
    static void accumulate(float& value, double& roundoff, double increment) {
        const double precise = double(value)+roundoff+increment;
        value = float(precise);
        roundoff = precise-value;
    }
    std::deque<Port> ports_;
    void sync_port(Port &p) const {
        p.gamma_value = gamma;
        p.gas_constant_value = gas_constant;
        p.pressure = pressure();
        p.density = density();
        p.temperature = temperature();
        p.enthalpy = gamma / (gamma - 1.0f) * p.pressure / std::max(p.density, 1e-12f);
        p.air_fraction = air_fraction();
    }
    void sync_ports() {
        for (Port &p : ports_)
            sync_port(p);
    }
};

struct GasCell {
    float rho{1.225f}, momentum{};
    float energy{1.225f * (101325.0f / (pipe_gm1 * 1.225f))};
    float air{1.0f};
};
struct GasState {
    float rho{1.225f}, velocity{}, pressure{101325.0f}, air{1.0f};
};

class Pipe1D {
public:
    // Integrated face fluxes: kg/s, N, W, kg/s of fresh air, positive along +x.
    struct Flux {
        float m{}, p{}, e{}, y{};
    };
    float area{}, dx{};
    float damping_multiplier{0.0f}, energy_relaxation_rate{0.0f}, ambient_pressure{101325.0f};
    std::vector<float> cell_areas, cell_lengths;
    std::vector<GasCell> cells;
    GasState left_ghost{}, right_ghost{};
    std::uint64_t positivity_repairs{};
    Pipe1D(float a, float length, int count)
        : Pipe1D(uniform(a, count), uniform(length / std::max(count, 1), count)) {}
    Pipe1D(std::vector<float> areas, std::vector<float> lengths)
        : area(areas.empty() ? 0 : areas.front()), dx(lengths.empty() ? 0 : lengths.front()),
          cell_areas(std::move(areas)), cell_lengths(std::move(lengths)), cells(cell_areas.size()),
          next_cells_(cells.size()), states_(cells.size()), waves_(cells.size()), faces_(cells.size() + 1) {
        if (cells.empty() || cell_lengths.size() != cells.size())
            throw std::invalid_argument("pipe profiles must be nonempty and have equal length");
        for (std::size_t i = 0; i < cells.size(); ++i)
            if (!std::isfinite(cell_areas[i]) || !std::isfinite(cell_lengths[i]) || cell_areas[i] <= 0 ||
                cell_lengths[i] <= 0)
                throw std::invalid_argument("pipe geometry must be finite and positive");
    }
    GasState state(int i) const { return primitive(cells.at(static_cast<std::size_t>(i))); }
    float cell_mass(int i) const { return cells.at(i).rho * cell_areas.at(i) * cell_lengths.at(i); }
    float cell_area(int i) const { return cell_areas.at(i); }
    float cell_length(int i) const { return cell_lengths.at(i); }
    float minimum_cell_length() const { return *std::min_element(cell_lengths.begin(), cell_lengths.end()); }
    float stable_timestep(float cfl = 0.8f) const {
        float dt = std::numeric_limits<float>::max();
        for (std::size_t i = 0; i < cells.size(); ++i) {
            auto s = primitive(cells[i]);
            if (!(s.rho > 0) || !(s.pressure > 0) || !std::isfinite(s.velocity) || !std::isfinite(s.pressure))
                throw std::runtime_error("nonphysical pipe state in CFL calculation");
            dt = std::min(dt, cfl * cell_lengths[i] / wave(s));
        }
        return dt;
    }
    // Solver preparation is valid until advance; boundary resolution changes
    // fluxes/ghosts, not cells. Standalone advance still prepares its own cache.
    float prepare_step(float cfl = 0.8f) {
        float dt = std::numeric_limits<float>::max();
        for (std::size_t i = 0; i < cells.size(); ++i) {
            states_[i] = primitive(cells[i]);
            waves_[i] = wave(states_[i]);
        }
        // Keep the arithmetic loop free of exceptional control flow so SIMD
        // can evaluate independent cells, then validate before any flux update.
        for (std::size_t i = 0; i < cells.size(); ++i) {
            const auto &s = states_[i];
            if (!(s.rho > 0) || !(s.pressure > 0) || !std::isfinite(s.velocity) || !std::isfinite(s.pressure))
                throw std::runtime_error("nonphysical pipe state in preparation");
            dt = std::min(dt, cfl * cell_lengths[i] / waves_[i]);
        }
        prepared_ = true;
        return dt;
    }
    int required_substeps(float dt, float cfl = 0.8f) const {
        if (!std::isfinite(dt) || dt < 0 || !std::isfinite(cfl) || cfl <= 0 || cfl > 1)
            throw std::invalid_argument("invalid timestep/CFL");
        const double n = std::ceil(double(dt) / stable_timestep(cfl));
        if (n > 1000000)
            throw std::runtime_error("excessive pipe substeps");
        return std::max(1, static_cast<int>(n));
    }
    void set_boundary_flux(bool left, Flux f) {
        if (left) {
            left_flux_ = f;
            left_override_ = true;
        } else {
            right_flux_ = f;
            right_override_ = true;
        }
    }
    static GasState primitive(const GasCell &c) {
        const float u = c.momentum / std::max(c.rho, 1e-12f);
        return {c.rho, u, pipe_gm1 * (c.energy - .5f * c.momentum * u), std::clamp(c.air, 0.0f, 1.0f)};
    }
    static float wave(const GasState &s) {
        return std::abs(s.velocity) + std::sqrt(gamma * s.pressure / s.rho);
    }
    static Flux riemann(const GasState &a, const GasState &b) {
        return riemann(a, b, std::max(wave(a), wave(b)));
    }
    void advance(float dt) {
        if (dt <= 0)
            return;
        if (!std::isfinite(dt))
            throw std::invalid_argument("nonfinite pipe timestep");
        const std::size_t n = cells.size();
        if (!prepared_)
            prepare_step();
        faces_[0] = left_override_ ? left_flux_ : scaled(riemann(left_ghost, states_[0]), cell_areas[0]);
        faces_[n] =
            right_override_ ? right_flux_ : scaled(riemann(states_[n - 1], right_ghost), cell_areas[n - 1]);
        for (std::size_t i = 1; i < n; ++i)
            faces_[i] = scaled(riemann(states_[i - 1], states_[i], std::max(waves_[i - 1], waves_[i])),
                               .5f * (cell_areas[i - 1] + cell_areas[i]));
        const float relaxation = energy_relaxation_rate > 0 ? std::exp(-energy_relaxation_rate * dt) : 1.0f;
        for (std::size_t i = 0; i < n; ++i) {
            const auto fl = faces_[i], fr = faces_[i + 1];
            const float factor = dt / (cell_areas[i] * cell_lengths[i]);
            const float al = i ? .5f * (cell_areas[i - 1] + cell_areas[i]) : cell_areas[i];
            const float ar = i + 1 < n ? .5f * (cell_areas[i] + cell_areas[i + 1]) : cell_areas[i];
            float rho = cells[i].rho - factor * (fr.m - fl.m);
            float mom = cells[i].momentum - factor * (fr.p - fl.p - states_[i].pressure * (ar - al));
            float en = cells[i].energy - factor * (fr.e - fl.e);
            const float air = cells[i].rho * cells[i].air - factor * (fr.y - fl.y);
            if (!std::isfinite(rho) || !std::isfinite(mom) || !std::isfinite(en) || !std::isfinite(air))
                throw std::runtime_error("nonfinite finite-volume update");
            if (rho < 1e-12f) {
                rho = 1e-12f;
                ++positivity_repairs;
            }
            if (damping_multiplier > 0)
                mom *= std::exp(-.1f * damping_multiplier * dt / cells[i].rho);
            if (energy_relaxation_rate > 0)
                en = ambient_pressure * pipe_gm1_inverse +
                     (en - ambient_pressure * pipe_gm1_inverse) * relaxation;
            const float minimum = 1e-6f * pipe_gm1_inverse + .5f * mom * mom / rho;
            if (en < minimum) {
                en = minimum;
                ++positivity_repairs;
            }
            // Purged tracer tails below the smallest normal float carry less
            // than 1.18e-38 of the gas mass. Avoid repeatedly transporting
            // subnormal concentrations; total gas mass and energy are unchanged.
            const float fraction = std::clamp(air / rho, 0.0f, 1.0f);
            next_cells_[i] = {rho, mom, en,
                             fraction < std::numeric_limits<float>::min() ? 0.0f : fraction};
        }
        cells.swap(next_cells_);
        prepared_ = false;
        left_override_ = right_override_ = false;
    }
    static Flux scaled(Flux f, float a) { return {f.m * a, f.p * a, f.e * a, f.y * a}; }

private:
    std::vector<GasCell> next_cells_;
    std::vector<GasState> states_;
    std::vector<float> waves_;
    std::vector<Flux> faces_;
    Flux left_flux_{}, right_flux_{};
    bool left_override_{}, right_override_{}, prepared_{};
    static std::vector<float> uniform(float v, int n) {
        if (n <= 0)
            throw std::invalid_argument("pipe count must be positive");
        return std::vector<float>(n, v);
    }
    static Flux riemann(const GasState &a, const GasState &b, float c) {
        const float ma = a.rho * a.velocity, mb = b.rho * b.velocity;
        const float ea = a.pressure * pipe_gm1_inverse + .5f * ma * a.velocity;
        const float eb = b.pressure * pipe_gm1_inverse + .5f * mb * b.velocity;
        return {.5f * (ma + mb - c * (b.rho - a.rho)),
                .5f * (ma * a.velocity + a.pressure + mb * b.velocity + b.pressure - c * (mb - ma)),
                .5f * ((ea + a.pressure) * a.velocity + (eb + b.pressure) * b.velocity - c * (eb - ea)),
                .5f * (ma * a.air + mb * b.air - c * (b.rho * b.air - a.rho * a.air))};
    }
};

class BoundarySystem {
public:
    explicit BoundarySystem(std::vector<Pipe1D> &pipes) : pipes_(pipes) {}
    void add_orifice(Port &port, int pipe, bool left, float &area, float cd = 1.0f) {
        if (std::abs(port.gamma_value - gamma) > 1e-6f ||
            std::abs(port.gas_constant_value - gas_constant) > .001f)
            throw std::invalid_argument("pipe coupling requires the same calorically perfect gas");
        if (!std::isfinite(cd) || cd < 0 || cd > 1)
            throw std::invalid_argument("orifice Cd must be in [0,1]");
        for (auto &o : orifices_)
            if (o.port == &port)
                throw std::invalid_argument("each connection requires a distinct component port");
        claim(pipe, left);
        Orifice o;
        o.port = &port;
        o.pipe = pipe;
        o.left = left;
        o.area = &area;
        o.cd = cd;
        if (port.owner_air_mass) {
            auto it = std::find_if(budgets_.begin(), budgets_.end(),
                                   [&](const Budget &b) { return b.air == port.owner_air_mass; });
            if (it == budgets_.end()) {
                o.budget = budgets_.size();
                budgets_.push_back({port.owner_air_mass, port.owner_exhaust_mass, port.owner_energy});
            } else
                o.budget = static_cast<std::size_t>(it - budgets_.begin());
        }
        orifices_.push_back(o);
    }
    void add_orifice_with_inertance(Port &p, int pipe, bool left, float &a, float cd, float length,
                                    float loss = 0) {
        if (!std::isfinite(length) || length <= 0 || !std::isfinite(loss) || loss < 0)
            throw std::invalid_argument("invalid inertance");
        add_orifice(p, pipe, left, a, cd);
        auto &o = orifices_.back();
        o.length = length;
        o.loss = loss;
    }
    void add_closed_end(int pipe, bool left) {
        claim(pipe, left);
        Orifice o;
        o.pipe = pipe;
        o.left = left;
        orifices_.push_back(o);
    }
    void add_open_end(int pipe, bool left, float p, float a) {
        if (!std::isfinite(p) || p <= 0 || !std::isfinite(a) || a <= 0)
            throw std::invalid_argument("invalid open boundary");
        claim(pipe, left);
        opens_.push_back({pipe, left, p, a});
    }
    // An ideal, zero-volume junction between two duct faces. Orientation is explicit.
    void add_connection(int a, bool left_a, int b, bool left_b) {
        if (a < 0 || b < 0 || static_cast<std::size_t>(a) >= pipes_.size() ||
            static_cast<std::size_t>(b) >= pipes_.size())
            throw std::invalid_argument("invalid connection pipe index");
        const float aa = pipes_[a].cell_area(end(pipes_[a], left_a));
        const float ab = pipes_[b].cell_area(end(pipes_[b], left_b));
        if (std::abs(aa - ab) > 1e-6f * std::max(aa, ab))
            throw std::invalid_argument("unequal-area junction needs a profiled pipe transition");
        if (a == b && left_a == left_b)
            throw std::invalid_argument("cannot join a face to itself");
        check_face(a, left_a);
        check_face(b, left_b);
        claim(a, left_a);
        claim(b, left_b);
        connections_.push_back({a, b, left_a, left_b});
    }
    float open_end_mass_flow(std::size_t i) const { return opens_.at(i).last_mass_flow; }
    float open_end_pressure(std::size_t i) const { return opens_.at(i).last_face_pressure; }
    std::uint64_t limited_fluxes{}, nonlinear_iterations{};
    void validate_topology() const {
        if (claims_.size() != 2 * pipes_.size())
            throw std::runtime_error("every pipe endpoint needs an explicit boundary or connection");
    }
    void resolve(float dt) {
        if (!std::isfinite(dt) || dt <= 0)
            throw std::invalid_argument("boundary timestep must be positive");
        for (auto &b : budgets_) {
            b.mass = .25f * (*b.air + (b.exhaust ? *b.exhaust : 0));
            b.heat = b.energy ? .25f * (*b.energy) : std::numeric_limits<float>::max();
        }
        for (auto &o : orifices_)
            resolve_orifice(o, dt);
        for (auto &o : opens_)
            resolve_open(o);
        for (auto &o : connections_)
            resolve_connection(o);
    }

private:
    struct Budget {
        float *air{}, *exhaust{}, *energy{};
        float mass{}, heat{};
    };
    struct Orifice {
        Port *port{};
        int pipe{};
        bool left{};
        float *area{};
        float cd{1}, length{}, loss{}, mdot{};
        std::size_t budget{std::numeric_limits<std::size_t>::max()};
        double previous_pressure{};
    };
    struct Open {
        int pipe;
        bool left;
        float pressure, area, last_mass_flow{}, last_face_pressure{};
    };
    struct Connection {
        int a, b;
        bool left_a, left_b;
    };
    struct Claim {
        int pipe;
        bool left;
    };
    std::vector<Pipe1D> &pipes_;
    std::vector<Orifice> orifices_;
    std::vector<Open> opens_;
    std::vector<Connection> connections_;
    std::vector<Budget> budgets_;
    std::vector<Claim> claims_;
    void check_face(int pipe, bool left) const {
        if (pipe < 0 || static_cast<std::size_t>(pipe) >= pipes_.size())
            throw std::invalid_argument("invalid boundary pipe index");
        for (auto c : claims_)
            if (c.pipe == pipe && c.left == left)
                throw std::invalid_argument("pipe face already connected");
    }
    void claim(int pipe, bool left) {
        check_face(pipe, left);
        claims_.push_back({pipe, left});
    }
    static int end(const Pipe1D &p, bool left) { return left ? 0 : static_cast<int>(p.cells.size()) - 1; }
    static void ghost(Pipe1D &p, bool left, GasState s) {
        if (left)
            p.left_ghost = s;
        else
            p.right_ghost = s;
    }
    void resolve_orifice(Orifice &o, float dt) {
        auto &p = pipes_[o.pipe];
        int i = end(p, o.left);
        auto s = p.state(i);
        float face_area = p.cell_area(i);
        if (o.area && (!std::isfinite(*o.area) || *o.area < 0))
            throw std::invalid_argument("invalid dynamic orifice area");
        if (o.port && (!(o.port->pressure > 0) || !(o.port->density > 0) ||
                       !std::isfinite(o.port->pressure) || !std::isfinite(o.port->density)))
            throw std::runtime_error("nonphysical component port");
        const float area = o.area ? std::min(std::max(*o.area, 0.f), face_area) : 0.f;
        if (!o.port || area <= 0 || o.cd == 0) {
            auto g = s;
            g.velocity = -s.velocity;
            ghost(p, o.left, g);
            auto wall_flux =
                Pipe1D::scaled(o.left ? Pipe1D::riemann(g, s) : Pipe1D::riemann(s, g), face_area);
            // A stationary adiabatic wall transfers momentum, never mass or energy.
            p.set_boundary_flux(o.left, {0, wall_flux.p, 0, 0});
            if (o.port)
                o.port->mass_flow = 0;
            o.mdot = 0;
            return;
        }
        auto &v = *o.port;
        if (std::abs(v.gamma_value - gamma) > 1e-6f || std::abs(v.gas_constant_value - gas_constant) > .001f)
            throw std::invalid_argument("pipe coupling requires the same calorically perfect gas");
        const float dp = v.pressure - s.pressure;
        const double g = gamma;
        const double critical = std::pow(2.0 / (g + 1), g / (g - 1));
        const double c = std::sqrt(g * s.pressure / s.rho);
        const double ui = s.velocity * (o.left ? 1 : -1);
        const auto face_density = [&](double pressure) {
            return s.rho * std::pow(pressure / s.pressure, 1 / g);
        };
        const auto nozzle_flow = [&](double pressure) {
            const bool from = v.pressure >= pressure;
            const double up = from ? v.pressure : pressure;
            const double down = from ? pressure : v.pressure;
            const double density = from ? v.density : face_density(pressure);
            const double ratio = std::max(down / up, critical);
            const double r = std::pow(ratio, 1 / g);
            const double speed2 = 2 * g / (g - 1) * up / density * (1 - std::pow(ratio, (g - 1) / g));
            return (from ? 1 : -1) * density * r * std::sqrt(std::max(0.0, speed2)) * area * o.cd;
        };
        // Match nozzle discharge to the outgoing duct characteristic. Imposing
        // sqrt(dp) against the *cell-centre* pressure is an unstable explicit
        // boundary near equilibrium and ignores the duct's acoustic impedance.
        const double zero_sound = c - .5 * (g - 1) * ui;
        if (zero_sound <= 0)
            throw std::runtime_error("orifice characteristic outside subsonic model");
        const double zero_pressure = s.pressure * std::pow(zero_sound / c, 2 * g / (g - 1));
        double low = std::min(zero_pressure, double(v.pressure));
        double high = std::max(zero_pressure, double(v.pressure));
        double duct_pressure = o.length == 0 && o.previous_pressure > low && o.previous_pressure < high
                                   ? o.previous_pressure : .5 * (low + high);
        const double flow_sign = v.pressure >= zero_pressure ? 1 : -1;
        const double flow_scale = area * o.cd;
        const double density_ratio_scale = std::pow(double(s.pressure) / v.pressure, 1 / g);
        const double critical_density_ratio = std::pow(critical, 1 / g);
        // Squared-flow residual avoids the nozzle sqrt singularity at dp=0.
        // Safeguarded analytic Newton; bisection remains the fallback.
        if (o.length == 0 && high > low)
            for (int iteration = 0; iteration < 40; ++iteration) {
                ++nonlinear_iterations;
                const double density = face_density(duct_pressure);
                const double sound = std::sqrt(g * duct_pressure / density);
                const double velocity = ui + 2 * (sound - c) / (g - 1);
                const double q = density * p.cell_area(i) * velocity;
                const double dq = p.cell_area(i) * density / (g * duct_pressure) * (velocity + sound);
                const bool from = v.pressure >= duct_pressure;
                const double up = from ? v.pressure : duct_pressure;
                const double down = from ? duct_pressure : v.pressure;
                const double upstream_density = from ? v.density : density;
                const double raw_ratio = down / up;
                const double ratio = std::max(raw_ratio, critical);
                const double root_ratio = raw_ratio <= critical
                                              ? critical_density_ratio
                                              : (from ? density / s.rho * density_ratio_scale
                                                      : s.rho / (density * density_ratio_scale));
                const double r1 = root_ratio * root_ratio, r2 = root_ratio * ratio;
                const double coefficient = flow_scale * flow_scale * 2 * g / (g - 1) * upstream_density * up;
                const double q2 = coefficient * (r1 - r2);
                const double dr = raw_ratio > critical ? (from ? 1 / up : -ratio / up) : 0;
                const double derivative = coefficient * ((2 / g * r1 - (g + 1) / g * r2) / ratio) * dr +
                                          (from ? 0 : q2 * (1 + 1 / g) / duct_pressure);
                const double residual = flow_sign * (q * q - q2);
                if (residual > 0)
                    high = duct_pressure;
                else
                    low = duct_pressure;
                if (high - low < std::max(1e-9, duct_pressure * 1e-12))
                    break;
                const double slope = flow_sign * (2 * q * dq - derivative);
                const double candidate = duct_pressure - residual / slope;
                if (std::abs(residual) < 4e-7 * std::max(q * q + q2, 1e-24))
                    break;
                duct_pressure = (std::isfinite(candidate) && candidate > low && candidate < high)
                                    ? candidate
                                    : .5 * (low + high);
            }
        o.previous_pressure = duct_pressure;
        float mdot = static_cast<float>(nozzle_flow(duct_pressure));
        const bool upstream_port = mdot >= 0;
        const float pup = upstream_port ? v.pressure : static_cast<float>(duct_pressure);
        const float ratio =
            static_cast<float>(std::max((upstream_port ? duct_pressure : v.pressure) / pup, critical));
        const float throat_rho =
            (upstream_port ? v.density : static_cast<float>(face_density(duct_pressure))) *
            std::pow(ratio, 1 / gamma);
        if (o.length > 0) {
            // d(mdot)/dt = A/L dp - K mdot |mdot|/(2 rho A L).
            const float donor_rho = o.mdot >= 0 ? v.density : s.rho;
            // Semi-implicit resistance remains dissipative for large dt.
            mdot = (o.mdot + area / o.length * dp * dt) /
                   (1 + dt * o.loss * std::abs(o.mdot) / (2 * donor_rho * area * o.length));
        }
        const bool from = mdot >= 0;
        const float h = from ? gamma / (gamma - 1) * v.pressure / v.density
                             : gamma / (gamma - 1) * s.pressure / s.rho + .5f * s.velocity * s.velocity;
        const float requested_mdot = mdot;
        if (from && o.budget < budgets_.size()) {
            auto &b = budgets_[o.budget];
            mdot = std::min(mdot, std::max(0.f, std::min(b.mass, b.heat / std::max(h, 1.f))) / dt);
            b.mass -= mdot * dt;
            b.heat -= mdot * h * dt;
        } else if (!from) {
            const float cell_energy = p.cells[i].energy * face_area * p.cell_length(i);
            mdot = std::max(mdot, -.25f * std::min(p.cell_mass(i), cell_energy / std::max(h, 1.f)) / dt);
        }
        if (mdot != requested_mdot)
            ++limited_fluxes;
        o.mdot = mdot;
        const float donor_rho = from ? v.density : s.rho;
        const float density = o.length > 0 ? donor_rho : throat_rho;
        const float velocity = mdot / (density * area);
        const float face_pressure =
            o.length > 0 ? .5f * (v.pressure + s.pressure) : static_cast<float>(duct_pressure);
        const float direction = o.left ? 1.f : -1.f;
        const float species = from ? v.air_fraction : s.air;
        const float signed_mass = direction * mdot;
        p.set_boundary_flux(o.left, {signed_mass,
                                     mdot * velocity + face_pressure * area + s.pressure * (face_area - area),
                                     signed_mass * h, signed_mass * species});
        ghost(p, o.left, {density, direction * velocity, face_pressure, species});
        v.mass_flow = -mdot;
        v.enthalpy = h;
        if (!from)
            v.air_fraction = s.air;
    }
    void resolve_open(Open &o) {
        auto &p = pipes_[o.pipe];
        int i = end(p, o.left);
        auto s = p.state(i);
        const float c = std::sqrt(gamma * s.pressure / s.rho), gm1 = gamma - 1;
        if (o.left ? s.velocity >= c : s.velocity <= -c)
            throw std::runtime_error(
                "supersonic inlet needs full upstream state; pressure-only boundary is insufficient");
        const float entropy = s.pressure / std::pow(s.rho, gamma);
        const float rho = std::pow(o.pressure / entropy, 1 / gamma);
        const float cb = std::sqrt(gamma * o.pressure / rho);
        const float u =
            o.left ? s.velocity - 2 * c / gm1 + 2 * cb / gm1 : s.velocity + 2 * c / gm1 - 2 * cb / gm1;
        const bool out = o.left ? s.velocity <= -c : s.velocity >= c;
        GasState g = out ? s : GasState{rho, u, o.pressure, (o.left ? u >= 0 : u <= 0) ? 1.f : s.air};
        ghost(p, o.left, g);
        const float a = std::min(o.area, p.cell_area(i));
        auto f = Pipe1D::scaled(o.left ? Pipe1D::riemann(g, s) : Pipe1D::riemann(s, g), a);
        f.p += s.pressure * (p.cell_area(i) - a);
        p.set_boundary_flux(o.left, f);
        o.last_mass_flow = f.m * (o.left ? -1 : 1);
        o.last_face_pressure = f.p / p.cell_area(i);
    }
    void resolve_connection(const Connection &o) {
        auto &a = pipes_[o.a];
        auto &b = pipes_[o.b];
        int ia = end(a, o.left_a), ib = end(b, o.left_b);
        auto sa = a.state(ia), sb = b.state(ib);
        const float da = o.left_a ? -1.f : 1.f, db = o.left_b ? 1.f : -1.f;
        sa.velocity *= da;
        sb.velocity *= db;
        const float area = std::min(a.cell_area(ia), b.cell_area(ib));
        auto f = Pipe1D::scaled(Pipe1D::riemann(sa, sb), area);
        a.set_boundary_flux(o.left_a,
                            {da * f.m, f.p + sa.pressure * (a.cell_area(ia) - area), da * f.e, da * f.y});
        b.set_boundary_flux(o.left_b,
                            {db * f.m, f.p + sb.pressure * (b.cell_area(ib) - area), db * f.e, db * f.y});
    }
};

} // namespace fluidsim::port