#pragma once

#include "FluidSim/Engine.hpp"
#include "FluidSim/OpenEndRadiation.hpp"
#include "FluidSim/PipeSolver.hpp"
#include <algorithm>
#include <array>
#include <span>
#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_;
        p.owner_fuel_mass = &fuel_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_ + fuel_mass_; }
    float internal_energy() const { return energy_; }
    float air_fraction() const { return air_mass_ / std::max(mass(), 1e-12f); }
    float fuel_fraction() const { return fuel_mass_ / std::max(mass(), 1e-12f); }
    float fuel_mass() const { return fuel_mass_; }
    void add_fuel_mass(float value, float temperature = 300.0f) {
        if (!std::isfinite(value) || value < 0.0f || !std::isfinite(temperature) || temperature <= 0.0f)
            throw std::invalid_argument("invalid fuel addition");
        fuel_mass_ += value;
        energy_ += value * gas_constant / (gamma - 1.0f) * temperature;
        sync_ports();
    }
    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();
        const float fuel = fuel_fraction();
        air_mass_ = total * fraction;
        fuel_mass_ = total * fuel;
        exhaust_mass_ = total * std::max(1.0f - fraction - fuel, 0.0f);
        energy_ = value * volume_ / (gamma - 1.0f);
        air_roundoff_ = exhaust_roundoff_ = fuel_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_fuel = 0, d_energy = 0;
        const float fraction = air_fraction();
        const float fuel = fuel_fraction();
        for (const Port &p : ports_) {
            const float flow_fraction = p.mass_flow >= 0.0f ? p.air_fraction : fraction;
            const float flow_fuel = p.mass_flow >= 0.0f ? p.fuel_fraction : fuel;
            d_air += double(p.mass_flow) * flow_fraction;
            d_fuel += double(p.mass_flow) * flow_fuel;
            d_exhaust += double(p.mass_flow) * std::max(1.0f - flow_fraction - flow_fuel, 0.0f);
            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(fuel_mass_, fuel_roundoff_, d_fuel * 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 || fuel_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);
        fuel_mass_ = std::max(fuel_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_{}, fuel_mass_{}, energy_;
    // Preserve sub-ulp transfers without changing the float port/budget ABI.
    double air_roundoff_{}, exhaust_roundoff_{}, fuel_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();
        p.fuel_fraction = fuel_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};
    float fuel{};
};
struct GasState {
    float rho{1.225f}, velocity{}, pressure{101325.0f}, air{1.0f};
    float fuel{};
};

class Pipe1D {
public:
    std::vector<float> linear_resistance, quadratic_resistance;
    // Integrated face fluxes: kg/s, N, W, kg/s of fresh air and fuel, positive along +x.
    struct Flux {
        float m{}, p{}, e{}, y{}, fuel{};
    };
    struct RoeState { double root_density{}, enthalpy{}; };
    PipeSolver solver{PipeSolver::Rusanov};
    std::uint64_t hllc_fallbacks{};
    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{};
    std::size_t scratch_bytes() const {
        return next_cells_.capacity()*sizeof(GasCell)+states_.capacity()*sizeof(GasState)+
            reconstructed_left_.capacity()*sizeof(GasState)+reconstructed_right_.capacity()*sizeof(GasState)+
            waves_.capacity()*sizeof(float)+faces_.capacity()*sizeof(Flux)+roe_states_.capacity()*sizeof(RoeState);
    }
    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) {
        if(!scratch_solver_initialized_ || scratch_solver_!=solver) {
            if(solver==PipeSolver::Muscl) {
                if(reconstructed_left_.size()!=cells.size()) reconstructed_left_.resize(cells.size());
                if(reconstructed_right_.size()!=cells.size()) reconstructed_right_.resize(cells.size());
            } else if(solver==PipeSolver::Hlle || solver==PipeSolver::Hllc) {
                if(roe_states_.size()!=cells.size()) roe_states_.resize(cells.size());
            }
            scratch_solver_=solver;
            scratch_solver_initialized_=true;
        }
        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]);
        }
        if (solver == PipeSolver::Hlle || solver == PipeSolver::Hllc)
            for (std::size_t i=0; i<cells.size(); ++i) {
                const auto& s=states_[i];
                const double energy=s.pressure*pipe_gm1_inverse+.5*double(s.rho)*s.velocity*s.velocity;
                roe_states_[i]={std::sqrt(s.rho),(energy+s.pressure)/s.rho};
            }
        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);
        const float air = std::clamp(c.air, 0.0f, 1.0f);
        return {c.rho, u, pipe_gm1 * (c.energy - .5f * c.momentum * u), air,
                std::clamp(c.fuel, 0.0f, 1.0f - air)};
    }
    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)));
    }
    // HLLE/HLLC use bounds enclosing endpoint and Roe acoustic speeds.
    template<PipeSolver method>
    static Flux flux_kernel(const GasState& a, const GasState& b,
                               float ca, float cb, bool* fallback = nullptr,
                               const RoeState* left_roe = nullptr, const RoeState* right_roe = nullptr) {
        if constexpr (method == PipeSolver::Rusanov)
            return riemann(a,b,std::max(std::abs(a.velocity)+ca,std::abs(b.velocity)+cb));
        const double ra=left_roe?left_roe->root_density:std::sqrt(a.rho);
        const double rb=right_roe?right_roe->root_density:std::sqrt(b.rho), sum=ra+rb;
        const double ea=a.pressure*pipe_gm1_inverse+.5*double(a.rho)*a.velocity*a.velocity;
        const double eb=b.pressure*pipe_gm1_inverse+.5*double(b.rho)*b.velocity*b.velocity;
        const double ua=a.velocity, ub=b.velocity;
        const double uroe=(ra*ua+rb*ub)/sum;
        const double ha=left_roe?left_roe->enthalpy:(ea+a.pressure)/a.rho;
        const double hb=right_roe?right_roe->enthalpy:(eb+b.pressure)/b.rho;
        const double hroe=(ra*ha+rb*hb)/sum;
        const double croe=std::sqrt(std::max(0.,double(gamma-1)*(hroe-.5*uroe*uroe)));
        const double sl=std::min({ua-ca,ub-cb,uroe-croe});
        const double sr=std::max({ua+ca,ub+cb,uroe+croe});
        const auto physical=[](const GasState& q,double energy) {
            const float m=q.rho*q.velocity;
            return Flux{m,m*q.velocity+q.pressure,float((energy+q.pressure)*q.velocity),m*q.air,m*q.fuel};
        };
        const Flux fa=physical(a,ea), fb=physical(b,eb);
        if(sl>=0) return fa;
        if(sr<=0) return fb;
        const auto hlle=[&] {
            const double inv=1/(sr-sl), diff=sl*sr;
            return Flux{float((sr*fa.m-sl*fb.m+diff*(b.rho-a.rho))*inv),
                        float((sr*fa.p-sl*fb.p+diff*(b.rho*ub-a.rho*ua))*inv),
                        float((sr*fa.e-sl*fb.e+diff*(eb-ea))*inv),
                        float((sr*fa.y-sl*fb.y+diff*(b.rho*b.air-a.rho*a.air))*inv),
                        float((sr*fa.fuel-sl*fb.fuel+diff*(b.rho*b.fuel-a.rho*a.fuel))*inv)};
        };
        if constexpr(method==PipeSolver::Hlle) return hlle();
        const double dl=a.rho*(sl-ua), dr=b.rho*(sr-ub);
        const double sm=(b.pressure-a.pressure+dl*ua-dr*ub)/(dl-dr);
        const double ps=a.pressure+dl*(sm-ua);
        // Reject nonphysical HLLC intermediate states without modifying conserved cells.
        const bool left=sm>=0;
        const auto& q=left?a:b;
        const auto f=left?fa:fb;
        const double speed=left?sl:sr;
        const double inv=1/(speed-sm);
        const double density=(left?dl:dr)*inv;
        const double energy=((speed-q.velocity)*(left?ea:eb)-q.pressure*q.velocity+ps*sm)*inv;
        if(!std::isfinite(sm) || !(sl<sm && sm<sr) || !(ps>0) || !(density>0) ||
           !(energy>.5*density*sm*sm)) {
            if(fallback) *fallback=true;
            return hlle();
        }
        const double mass=f.m+speed*(density-q.rho);
        return {float(mass),float(f.p+speed*(density*sm-q.rho*q.velocity)),
                float(f.e+speed*(energy-(left?ea:eb))),float(mass*q.air),float(mass*q.fuel)};
    }
    static Flux numerical_flux(const GasState& a,const GasState& b,PipeSolver solver,
                               float ca,float cb,bool* fallback=nullptr) {
        if(solver==PipeSolver::Rusanov) return flux_kernel<PipeSolver::Rusanov>(a,b,ca,cb,fallback);
        if(solver==PipeSolver::Hlle) return flux_kernel<PipeSolver::Hlle>(a,b,ca,cb,fallback);
        return flux_kernel<PipeSolver::Hllc>(a,b,ca,cb,fallback);
    }
    static Flux numerical_flux(const GasState& a,const GasState& b,PipeSolver solver) {
        return numerical_flux(a,b,solver,std::sqrt(gamma*a.pressure/a.rho),std::sqrt(gamma*b.pressure/b.rho));
    }
    void set_acoustic_neighbor(bool left,GasState state) {
        if(left) { left_acoustic_neighbor_=state; has_left_acoustic_neighbor_=true; }
        else { right_acoustic_neighbor_=state; has_right_acoustic_neighbor_=true; }
    }
    void prepare_acoustic_reconstruction(float dt) {
        if(solver!=PipeSolver::Muscl) return;
        if(!prepared_) prepare_step();
        reconstruct(dt); reconstructed_prepared_=true;
    }
    GasState acoustic_edge(bool left) const {
        return left?reconstructed_left_.front():reconstructed_right_.back();
    }
    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(numerical_flux(left_ghost, states_[0], solver), cell_areas[0]);
        faces_[n] =
            right_override_ ? right_flux_ : scaled(numerical_flux(states_[n - 1], right_ghost, solver), cell_areas[n - 1]);
        if(solver==PipeSolver::Muscl) {
            if(!reconstructed_prepared_) reconstruct(dt);
            for(std::size_t i=1;i<n;++i) {
                bool fallback=false;
                const auto& a=reconstructed_right_[i-1];
                const auto& b=reconstructed_left_[i];
                faces_[i]=scaled(numerical_flux(a,b,PipeSolver::Hllc,
                    std::sqrt(gamma*a.pressure/a.rho),std::sqrt(gamma*b.pressure/b.rho),&fallback),
                    .5f*(cell_areas[i-1]+cell_areas[i]));
                hllc_fallbacks+=fallback;
            }
        } else if (solver == PipeSolver::Rusanov) {
            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]));
        } else {
            const auto fill_faces=[&]<PipeSolver method>() {
                for (std::size_t i = 1; i < n; ++i) {
                    bool fallback=false;
                    faces_[i] = scaled(flux_kernel<method>(states_[i-1], states_[i],
                        waves_[i-1]-std::abs(states_[i-1].velocity), waves_[i]-std::abs(states_[i].velocity),
                        &fallback, &roe_states_[i-1], &roe_states_[i]),
                        .5f*(cell_areas[i-1]+cell_areas[i]));
                    hllc_fallbacks += fallback;
                }
            };
            if(solver==PipeSolver::Hlle) fill_faces.template operator()<PipeSolver::Hlle>();
            else fill_faces.template operator()<PipeSolver::Hllc>();
        }
        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);
            const float fuel = cells[i].rho * cells[i].fuel - factor * (fr.fuel - fl.fuel);
            if (!std::isfinite(rho) || !std::isfinite(mom) || !std::isfinite(en) ||
                !std::isfinite(air) || !std::isfinite(fuel))
                throw std::runtime_error("nonfinite finite-volume update");
            if (rho < 1e-12f) {
                rho = 1e-12f;
                ++positivity_repairs;
            }
            if (!linear_resistance.empty()) {
                const float speed = std::abs(mom/rho);
                mom /= 1 + dt*(linear_resistance[i] + .5f*quadratic_resistance[i]*speed);
                // Retain total energy: dissipated kinetic energy becomes heat.
            }
            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);
            const float fuel_fraction = std::clamp(fuel / rho, 0.0f, 1.0f - fraction);
            next_cells_[i] = {rho, mom, en,
                             fraction < std::numeric_limits<float>::min() ? 0.0f : fraction,
                             fuel_fraction < std::numeric_limits<float>::min() ? 0.0f : fuel_fraction};
        }
        cells.swap(next_cells_);
        prepared_ = false;
        reconstructed_prepared_=has_left_acoustic_neighbor_=has_right_acoustic_neighbor_=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, f.fuel * a};
    }

private:
    std::vector<GasCell> next_cells_;
    std::vector<GasState> states_;
    std::vector<GasState> reconstructed_left_, reconstructed_right_;
    std::vector<float> waves_;
    std::vector<Flux> faces_;
    std::vector<RoeState> roe_states_;
    PipeSolver scratch_solver_{PipeSolver::Rusanov};
    bool scratch_solver_initialized_{};
    Flux left_flux_{}, right_flux_{};
    bool left_override_{}, right_override_{}, prepared_{};
    GasState left_acoustic_neighbor_{}, right_acoustic_neighbor_{};
    bool has_left_acoustic_neighbor_{}, has_right_acoustic_neighbor_{}, reconstructed_prepared_{};
    // Primitive minmod MUSCL-Hancock on uniform sections. Component interfaces
    // and area/profile transitions retain the existing first-order coupling.
    void reconstruct(float dt) {
        const auto minmod=[](float a,float b) {
            return a*b<=0 ? 0.f : std::copysign(std::min(std::abs(a),std::abs(b)),a);
        };
        for(std::size_t i=0;i<states_.size();++i) {
            const auto q=states_[i];
            reconstructed_left_[i]=reconstructed_right_[i]=q;
            if((i==0 && !has_left_acoustic_neighbor_) ||
               (i+1==states_.size() && !has_right_acoustic_neighbor_) ||
               (i>0 && (cell_areas[i-1]!=cell_areas[i] || cell_lengths[i-1]!=cell_lengths[i])) ||
               (i+1<states_.size() && (cell_areas[i+1]!=cell_areas[i] || cell_lengths[i+1]!=cell_lengths[i]))) continue;
            const auto a=i==0?left_acoustic_neighbor_:states_[i-1];
            const auto b=i+1==states_.size()?right_acoustic_neighbor_:states_[i+1];
            const float dr=minmod(q.rho-a.rho,b.rho-q.rho);
            const float du=minmod(q.velocity-a.velocity,b.velocity-q.velocity);
            const float dp=minmod(q.pressure-a.pressure,b.pressure-q.pressure);
            const float da=minmod(q.air-a.air,b.air-q.air);
            const float df=minmod(q.fuel-a.fuel,b.fuel-q.fuel);
            const float h=.5f*dt/cell_lengths[i];
            const GasState half{q.rho-h*(q.velocity*dr+q.rho*du),
                q.velocity-h*(q.velocity*du+dp/q.rho),
                q.pressure-h*(q.velocity*dp+gamma*q.pressure*du),
                q.air-h*q.velocity*da,q.fuel-h*q.velocity*df};
            const auto edge=[&](float sign) { return GasState{half.rho+sign*dr,
                half.velocity+sign*du,half.pressure+sign*dp,half.air+sign*da,half.fuel+sign*df}; };
            const auto left=edge(-.5f), right=edge(.5f);
            const auto physical=[](const GasState& v) {
                return v.rho>1e-12f && v.pressure>1e-6f && std::isfinite(v.pressure) &&
                    std::isfinite(v.rho) && std::isfinite(v.velocity) && v.air>=0 && v.fuel>=0 && v.air+v.fuel<=1.000001f;
            };
            if(physical(left) && physical(right)) {
                reconstructed_left_[i]=left; reconstructed_right_[i]=right;
            }
        }
    }
    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)),
                .5f * (ma * a.fuel + mb * b.fuel - c * (b.rho * b.fuel - a.rho * a.fuel))};
    }
};

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_fuel_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});
    }
    // Full circular mouth radiating into free space; no flange or mean-flow correction.
    void add_radiating_open_end(int pipe, bool left, float pressure, float area, float temperature = 300) {
        if (!std::isfinite(temperature) || temperature <= 0)
            throw std::invalid_argument("invalid ambient temperature");
        if (pipe < 0 || static_cast<std::size_t>(pipe) >= pipes_.size())
            throw std::invalid_argument("invalid boundary pipe index");
        const float face_area = pipes_[pipe].cell_area(end(pipes_[pipe], left));
        if (std::abs(area-face_area) > 1e-6f*face_area)
            throw std::invalid_argument("radiating mouth must match pipe face area");
        add_open_end(pipe, left, pressure, area);
        opens_.back().radiating = true;
        opens_.back().temperature = temperature;
    }
    // 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});
    }
    struct JunctionFace {
        int pipe; bool left;
        // Unit vector pointing OUT of the node into this branch, in world axes.
        std::array<float,3> direction{};
    };
    enum class JunctionMethod { Characteristic, LegacyRusanov };
    void add_junction(std::initializer_list<JunctionFace> ends,
                      JunctionMethod method=JunctionMethod::Characteristic,float turning_loss=0) {
        add_junction(std::span(ends.begin(),ends.size()),method,turning_loss);
    }
    void add_junction(std::span<const JunctionFace> ends,
                      JunctionMethod method=JunctionMethod::Characteristic,float turning_loss=0) {
        if(ends.size()<2 || !std::isfinite(turning_loss) || turning_loss<0 || turning_loss>10)
            throw std::invalid_argument("invalid junction size or turning loss");
        JunctionNode node; node.method=method; node.turning_loss=turning_loss;
        for(auto e:ends) {
            check_face(e.pipe,e.left);
            for(auto previous:node.ends)
                if(previous.pipe==e.pipe && previous.left==e.left) throw std::invalid_argument("duplicate junction face");
            double norm=0;
            for(float v:e.direction) { if(!std::isfinite(v)) throw std::invalid_argument("invalid junction direction"); norm+=double(v)*v; }
            if(norm>0) for(float& v:e.direction) v/=std::sqrt(norm);
            else if(turning_loss>0) throw std::invalid_argument("direction required for turning loss");
            node.ends.push_back(e);
        }
        for(auto e:ends) claim(e.pipe,e.left);
        node.states.resize(ends.size()); node.faces.resize(ends.size());
        nodes_.push_back(std::move(node));
    }
    float open_end_mass_flow(std::size_t i) const { return opens_.at(i).last_mass_flow; }
    float open_end_volume_flow(std::size_t i) const { return opens_.at(i).last_volume_flow; }
    float open_end_pressure(std::size_t i) const { return opens_.at(i).last_face_pressure; }
    std::uint64_t limited_fluxes{}, nonlinear_iterations{}, junction_fallbacks{}, junction_iterations{}, junction_limited_fluxes{};
    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.fuel ? *b.fuel : 0));
            b.heat = b.energy ? .25f * (*b.energy) : std::numeric_limits<float>::max();
        }
        // Supply matching direct-pipe neighbors before reconstruction, preserving
        // the same stencil across artificial pipe splits and all orientations.
        for(const auto& o:connections_) {
            auto& a=pipes_[o.a]; auto& b=pipes_[o.b];
            const int ia=end(a,o.left_a), ib=end(b,o.left_b);
            if(a.solver==PipeSolver::Muscl && b.solver==PipeSolver::Muscl &&
               a.cell_area(ia)==b.cell_area(ib) && a.cell_length(ia)==b.cell_length(ib)) {
                auto sa=a.state(ia), sb=b.state(ib);
                const float sign=o.left_a==o.left_b?-1.f:1.f;
                sa.velocity*=sign; sb.velocity*=sign;
                a.set_acoustic_neighbor(o.left_a,sb); b.set_acoustic_neighbor(o.left_b,sa);
            }
        }
        for(auto& pipe:pipes_) pipe.prepare_acoustic_reconstruction(dt);
        for (auto &o : orifices_)
            resolve_orifice(o, dt);
        for (auto &o : opens_)
            resolve_open(o, dt);
        for (auto &o : connections_)
            resolve_connection(o);
        for(auto& node:nodes_) resolve_junction(node,dt);
    }

private:
    struct JunctionNode {
        std::vector<JunctionFace> ends;
        std::vector<GasState> states;
        struct Face { double pressure{},rho{},velocity{},mass{},enthalpy{},air{},fuel{},loss{}; };
        std::vector<Face> faces;
        JunctionMethod method{JunctionMethod::Characteristic};
        double turning_loss{};
    };
    std::vector<JunctionNode> nodes_;
    // Zero-storage Rusanov star state. Integrated mass, total energy and
    // species fluxes sum to zero algebraically across every connected face.
    // The node has no conserved state, timestep or storage volume of its own.
    void resolve_junction(JunctionNode& node,float dt) {
        if(node.method==JunctionMethod::Characteristic && characteristic_junction(node,dt)) return;
        if(node.method==JunctionMethod::Characteristic) ++junction_fallbacks;
        legacy_junction(node);
    }
    void legacy_junction(JunctionNode& node) {
        double speed=0;
        for(std::size_t k=0;k<node.ends.size();++k) {
            auto e=node.ends[k];auto& p=pipes_[e.pipe];auto s=p.state(end(p,e.left));
            s.velocity*=e.left?-1:1;node.states[k]=s;
            speed=std::max(speed,2*(std::sqrt(gamma*s.pressure/s.rho)+std::abs(double(s.velocity))));
        }
        double area_sum=0,rho=0,energy=0,air=0,fuel=0;
        for(std::size_t k=0;k<node.ends.size();++k) {
            auto e=node.ends[k];auto& p=pipes_[e.pipe];const auto s=node.states[k];
            const double a=p.cell_area(end(p,e.left)), u=s.velocity;
            const double en=s.pressure/pipe_gm1+.5*s.rho*u*u;
            area_sum+=a;rho+=a*s.rho*(1+u/speed);
            energy+=a*(en+(en+s.pressure)*u/speed);
            air+=a*s.rho*s.air*(1+u/speed);fuel+=a*s.rho*s.fuel*(1+u/speed);
        }
        rho/=area_sum;energy/=area_sum;air/=area_sum;fuel/=area_sum;
        const double pressure=energy*pipe_gm1;
        for(std::size_t k=0;k<node.ends.size();++k) {
            auto e=node.ends[k];auto& p=pipes_[e.pipe];const auto s=node.states[k];
            const double a=p.cell_area(end(p,e.left)), u=s.velocity;
            const double en=s.pressure/pipe_gm1+.5*s.rho*u*u;
            const double orientation=e.left?-1:1;
            const double m=.5*a*(s.rho*u-speed*(rho-s.rho));
            const double heat=.5*a*((en+s.pressure)*u-speed*(energy-en));
            const double y=.5*a*(s.rho*s.air*u-speed*(air-s.rho*s.air));
            const double f=.5*a*(s.rho*s.fuel*u-speed*(fuel-s.rho*s.fuel));
            const double momentum=.5*a*(s.rho*u*u+s.pressure+pressure+speed*s.rho*u);
            p.set_boundary_flux(e.left,{float(orientation*m),float(momentum),float(orientation*heat),float(orientation*y),float(orientation*f)});
        }
    }
    // Nonlinear outgoing acoustic-wave curves (shock or rarefaction), with
    // common node pressure, conservative stagnation-enthalpy/species mixing,
    // and outgoing contact density set from the mixed enthalpy. This is a
    // compact subsonic engineering closure, not a resolved 3D junction.
    bool characteristic_junction(JunctionNode& node,float dt) {
        double low=std::numeric_limits<double>::max(), high=0;
        for(std::size_t k=0;k<node.ends.size();++k) {
            const auto e=node.ends[k]; const auto& pipe=pipes_[e.pipe];
            auto state=pipe.state(end(pipe,e.left));
            state.velocity*=e.left?-1:1; // positive INTO the node
            const double c=std::sqrt(gamma*state.pressure/state.rho);
            if(std::abs(state.velocity)>=.95*c) return false;
            node.states[k]=state;
            low=std::min(low,double(state.pressure)); high=std::max(high,double(state.pressure));
            node.faces[k].loss=0;
        }
        const auto wave_curve=[](const GasState& q,double p) {
            if(p>q.pressure) {
                const double a=2/((gamma+1)*q.rho), b=(gamma-1)/(gamma+1)*q.pressure;
                return (p-q.pressure)*std::sqrt(a/(p+b));
            }
            const double c=std::sqrt(gamma*q.pressure/q.rho);
            return 2*c/(gamma-1)*(std::pow(p/q.pressure,(gamma-1)/(2*gamma))-1);
        };
        const auto density=[](const GasState& q,double p) {
            const double ratio=p/q.pressure, b=(gamma-1)/(gamma+1);
            return p>q.pressure ? q.rho*(ratio+b)/(b*ratio+1) : q.rho*std::pow(ratio,1/gamma);
        };
        // Each pressure evaluation writes scratch face states only, never pipes.
        const auto residual=[&](double pressure) {
            double incoming=0, heat=0, air=0, fuel=0;
            for(std::size_t k=0;k<node.ends.size();++k) {
                const auto q=node.states[k]; auto& f=node.faces[k];
                double p=pressure;
                for(int iteration=0;iteration<(f.loss>0?8:1);++iteration) {
                    const double u=q.velocity-wave_curve(q,p), rho=density(q,p);
                    const double target=std::clamp(pressure+std::copysign(.5*f.loss*rho*u*u,u),pressure*.05,pressure*20);
                    p=.5*(p+target);
                }
                f.pressure=p; f.rho=density(q,p); f.velocity=q.velocity-wave_curve(q,p);
                const auto e=node.ends[k]; const auto& pipe=pipes_[e.pipe];
                f.mass=pipe.cell_area(end(pipe,e.left))*f.rho*f.velocity;
                f.enthalpy=gamma/(gamma-1)*p/f.rho+.5*f.velocity*f.velocity;
                f.air=q.air; f.fuel=q.fuel;
                if(f.mass>0) { incoming+=f.mass; heat+=f.mass*f.enthalpy; air+=f.mass*q.air; fuel+=f.mass*q.fuel; }
            }
            double balance=0;
            for(std::size_t k=0;k<node.ends.size();++k) {
                auto& f=node.faces[k];
                if(f.mass<0 && incoming>1e-20) {
                    const double h=heat/incoming, thermal=h-.5*f.velocity*f.velocity;
                    if(!(thermal>0)) return std::numeric_limits<double>::quiet_NaN();
                    f.rho=gamma/(gamma-1)*f.pressure/thermal;
                    const auto e=node.ends[k]; const auto& pipe=pipes_[e.pipe];
                    f.mass=pipe.cell_area(end(pipe,e.left))*f.rho*f.velocity;
                    f.enthalpy=h; f.air=air/incoming; f.fuel=fuel/incoming;
                }
                balance+=f.mass;
            }
            return balance;
        };
        // Iterate the donor-weighted direction coefficients separately from the
        // scalar pressure root; K=0 remains the ideal acoustic scattering case.
        double node_pressure=.5*(low+high);
        for(int geometry=0;geometry<(node.turning_loss>0?3:1);++geometry) {
            double lo=low*.01, hi=high*100;
            double fl=residual(lo), fh=residual(hi);
            if(!std::isfinite(fl) || !std::isfinite(fh) || fl<0 || fh>0) return false;
            node_pressure=std::clamp(node_pressure,lo,hi);
            bool converged=false;
            for(int iteration=0;iteration<40;++iteration) {
                const double value=residual(node_pressure); ++junction_iterations;
                if(!std::isfinite(value)) return false;
                double magnitude=0; for(const auto& face:node.faces) magnitude+=std::abs(face.mass);
                if(std::abs(value)<1e-14+1e-8*magnitude) { converged=true; break; }
                if(value>0) lo=node_pressure; else hi=node_pressure;
                const double delta=std::max(1e-3,node_pressure*1e-5);
                const double derivative=(residual(node_pressure+delta)-value)/delta;
                const double next=node_pressure-value/derivative;
                node_pressure=std::isfinite(next) && next>lo && next<hi ? next : .5*(lo+hi);
            }
            if(!converged) return false;
            residual(node_pressure);
            if(geometry+1<(node.turning_loss>0?3:1)) {
                for(std::size_t k=0;k<node.ends.size();++k) {
                    auto& f=node.faces[k]; double weight=0, loss=0;
                    if(f.mass<0) for(std::size_t j=0;j<node.ends.size();++j) {
                        const double m=std::max(0.,node.faces[j].mass); double dot=0;
                        for(int axis=0;axis<3;++axis) dot+=node.ends[k].direction[axis]*node.ends[j].direction[axis];
                        weight+=m; loss+=m*.5*(1+std::clamp(dot,-1.,1.));
                    }
                    f.loss=weight>0?node.turning_loss*loss/weight:0;
                }
            }
        }
        double incoming=0, outgoing=0, heat=0, air=0, fuel=0, scale=1;
        for(std::size_t k=0;k<node.ends.size();++k) {
            const auto& f=node.faces[k];
            if(!std::isfinite(f.mass) || !std::isfinite(f.enthalpy)) return false;
            if(std::abs(f.velocity)>=.98*std::sqrt(gamma*f.pressure/f.rho)) return false;
            if(f.mass>0) { incoming+=f.mass;heat+=f.mass*f.enthalpy;air+=f.mass*f.air;fuel+=f.mass*f.fuel; }
            else outgoing-=f.mass;
        }
        if(incoming>1e-16 && outgoing<=1e-16) {
            if(incoming>1e-11) return false;
            // Residual tolerance at rest must not create a one-sided mass source.
            for(auto& face:node.faces) face.mass=0;
            incoming=heat=air=fuel=0;
        }
        // Remove root roundoff by scaling all outflows together. A shared donor
        // budget scales advective fluxes together, preserving node conservation.
        for(std::size_t k=0;k<node.ends.size();++k) {
            auto& f=node.faces[k];
            if(f.mass<0) { f.mass*=outgoing>0?incoming/outgoing:0; f.enthalpy=incoming>0?heat/incoming:0; f.air=incoming>0?air/incoming:0; f.fuel=incoming>0?fuel/incoming:0; }
            if(f.mass>0) {
                const auto e=node.ends[k]; const auto& pipe=pipes_[e.pipe];const int i=end(pipe,e.left);
                scale=std::min(scale,.25*pipe.cell_mass(i)/(f.mass*dt));
                scale=std::min(scale,.25*double(pipe.cells[i].energy)*pipe.cell_area(i)*pipe.cell_length(i)/(f.mass*f.enthalpy*dt));
            }
        }
        if(scale<1) { ++limited_fluxes; ++junction_limited_fluxes; }
        for(std::size_t k=0;k<node.ends.size();++k) {
            const auto e=node.ends[k]; auto& pipe=pipes_[e.pipe]; const auto& f=node.faces[k];
            const double orientation=e.left?-1:1, m=f.mass*scale;
            const double momentum=m*f.velocity+f.pressure*pipe.cell_area(end(pipe,e.left));
            pipe.set_boundary_flux(e.left,{float(orientation*m),float(momentum),float(orientation*m*f.enthalpy),float(orientation*m*f.air),float(orientation*m*f.fuel)});
        }
        return true;
    }
    struct Budget {
        float *air{}, *exhaust{}, *fuel{}, *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{}, last_volume_flow{};
        bool radiating{};
        float temperature{300};
        OpenEndRadiation reflection{};
    };
    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 fuel = from ? v.fuel_fraction : s.fuel;
        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, signed_mass * fuel});
        ghost(p, o.left, {density, direction * velocity, face_pressure, species, fuel});
        v.mass_flow = -mdot;
        v.enthalpy = h;
        if (!from)
            v.air_fraction = s.air;
        if (!from)
            v.fuel_fraction = s.fuel;
    }
    void resolve_open(Open &o, float dt) {
        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");
        if (o.radiating) {
            const float direction = o.left ? -1.f : 1.f;
            const double outward = direction*s.velocity;
            GasState face = s;
            if (outward < c) {
                const double reference_c = c*std::pow(double(o.pressure)/s.pressure, (gamma-1)/(2*gamma));
                const double incident = outward + 2*(c-reference_c)/(gamma-1);
                const double radius = std::sqrt(o.area/pi);
                const double ambient_c = std::sqrt(gamma*gas_constant*o.temperature);
                const double reflected = o.reflection.process(incident, dt*ambient_c/radius);
                const double sound = reference_c + (gamma-1)*.25*(incident+reflected);
                if (!(sound > 0)) throw std::runtime_error("radiation characteristic outside subsonic model");
                face.pressure = o.pressure*std::pow(sound/reference_c, 2*gamma/(gamma-1));
                face.velocity = direction*.5*(incident-reflected);
                face.rho = s.rho*std::pow(face.pressure/s.pressure, 1/gamma);
                if (direction*face.velocity < 0) {
                    // Fresh ambient gas enters with ambient stagnation enthalpy.
                    const float temperature = o.temperature - .5f*face.velocity*face.velocity /
                                              (gamma/(gamma-1)*gas_constant);
                    if (!(temperature > 0)) throw std::runtime_error("ambient inflow exceeds stagnation enthalpy");
                    face.rho = face.pressure/(gas_constant*temperature);
                    face.air = 1; face.fuel = 0;
                }
            }
            ghost(p, o.left, face);
            const float mass_flux = face.rho*face.velocity*o.area;
            const float enthalpy = gamma/(gamma-1)*face.pressure/face.rho + .5f*face.velocity*face.velocity;
            p.set_boundary_flux(o.left, {mass_flux, mass_flux*face.velocity+face.pressure*o.area,
                                        mass_flux*enthalpy, mass_flux*face.air, mass_flux*face.fuel});
            o.last_volume_flow = direction*face.velocity*o.area;
            o.last_mass_flow = direction*mass_flux;
            o.last_face_pressure = face.pressure;
            return;
        }
        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;
        const bool ambient_inflow = o.left ? u >= 0 : u <= 0;
        GasState g = out ? s : GasState{rho, u, o.pressure, ambient_inflow ? 1.f : s.air,
                                        ambient_inflow ? 0.f : s.fuel};
        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_volume_flow = o.last_mass_flow / (ambient_inflow ? rho : s.rho);
        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);
        if(a.solver==PipeSolver::Muscl && b.solver==PipeSolver::Muscl) {
            sa=a.acoustic_edge(o.left_a); sb=b.acoustic_edge(o.left_b);
        }
        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::numerical_flux(sa, sb, a.solver == b.solver ? a.solver : PipeSolver::Rusanov), 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, da * f.fuel});
        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, db * f.fuel});
    }
};

} // namespace fluidsim::port