#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)) {}

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

    Port& create_port() {
        ports_.emplace_back();
        auto& p = ports_.back();
        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 air_fraction() const { return air_mass_ / std::max(mass(), 1e-12f); }
    void set_pressure(float value, float temp = -1.0f) {
        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);
        sync_ports();
    }
    void update(float dt) {
        if (dt <= 0.0f) return;
        float d_air = 0.0f, d_exhaust = 0.0f, d_energy = 0.0f;
        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 += p.mass_flow * flow_fraction;
            d_exhaust += p.mass_flow * (1.0f - flow_fraction);
            d_energy += p.mass_flow * p.enthalpy;
        }
        air_mass_ += d_air * dt;
        exhaust_mass_ += d_exhaust * dt;
        energy_ += (d_energy - pressure() * volume_change_rate) * dt;
        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);
        }
        if (mass() < 1e-9f || energy_ < 0.0f) {
            air_mass_ = 1e-9f;
            exhaust_mass_ = 0.0f;
            energy_ = ambient_pressure * volume_ / (gamma - 1.0f);
        }
        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_;
    std::deque<Port> ports_;
    void sync_port(Port& p) const {
        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:
    float area{}, dx{};
    float damping_multiplier{1.0f}, energy_relaxation_rate{0.5f}, ambient_pressure{101325.0f};
    std::vector<float> cell_areas, cell_lengths;
    std::vector<GasCell> cells;
    GasState left_ghost{}, right_ghost{};
    Pipe1D(float area_m2, float length_m, int cell_count)
        : area(area_m2), dx(length_m / static_cast<float>(cell_count)),
          cell_areas(static_cast<std::size_t>(cell_count), area_m2),
          cell_lengths(static_cast<std::size_t>(cell_count), dx),
          cells(static_cast<std::size_t>(cell_count)), next_cells_(static_cast<std::size_t>(cell_count)) {}
    Pipe1D(std::vector<float> areas_m2, std::vector<float> lengths_m)
        : area(areas_m2.empty() ? 0.0f : areas_m2.front()),
          dx(lengths_m.empty() ? 0.0f : lengths_m.front()),
          cell_areas(std::move(areas_m2)), cell_lengths(std::move(lengths_m)),
          cells(cell_areas.size()), next_cells_(cell_areas.size()) {
        if (cell_areas.empty() || cell_lengths.size() != cell_areas.size())
            throw std::invalid_argument("pipe area and length profiles must be non-empty and have equal length");
        for (std::size_t i = 0; i < cell_areas.size(); ++i) {
            if (cell_areas[i] <= 0.0f || cell_lengths[i] <= 0.0f)
                throw std::invalid_argument("pipe cell area and length must be positive");
        }
    }
    GasState state(int idx) const {
        const auto& c = cells.at(static_cast<std::size_t>(idx));
        const float rho = c.rho;
        const float u = c.momentum / std::max(rho, 1e-12f);
        const float p = pipe_gm1 * (c.energy - 0.5f * c.momentum * u);
        return {rho, u, p, std::clamp(c.air, 0.0f, 1.0f)};
    }
    float cell_mass(int idx) const {
        const auto cell = static_cast<std::size_t>(idx);
        return cells.at(cell).rho * cell_areas.at(cell) * cell_lengths.at(cell);
    }
    float cell_area(int idx) const { return cell_areas.at(static_cast<std::size_t>(idx)); }
    float cell_length(int idx) const { return cell_lengths.at(static_cast<std::size_t>(idx)); }
    float minimum_cell_length() const {
        return *std::min_element(cell_lengths.begin(), cell_lengths.end());
    }
    int required_substeps(float dt, float target_cfl = 0.8f) const {
        float maximum_wave_speed = 0.0f;
        for (int i = 0; i < static_cast<int>(cells.size()); ++i) {
            const GasState s = state(i);
            maximum_wave_speed = std::max(maximum_wave_speed,
                std::abs(s.velocity) + std::sqrt(gamma * s.pressure / std::max(s.rho, 1e-12f)));
        }
        return std::max(1, static_cast<int>(std::ceil(dt * maximum_wave_speed /
            (target_cfl * std::max(minimum_cell_length(), 1e-12f)))));
    }
    void advance(float dt) {
        const int n = static_cast<int>(cells.size());
        if (n == 0 || dt <= 0.0f) return;
        const float ambient_energy = ambient_pressure * pipe_gm1_inverse;
        for (int i = 0; i < n; ++i) {
            const GasState l = i == 0 ? left_ghost : flux_state(i - 1), r = flux_state(i);
            const GasState a = flux_state(i), b = i == n - 1 ? right_ghost : flux_state(i + 1);
            const Flux fl = flux(l, r), fr = flux(a, b);
            const float factor = dt / cell_lengths[static_cast<std::size_t>(i)];
            float rho = cells[i].rho - factor * (fr.m - fl.m);
            float mom = cells[i].momentum - factor * (fr.p - fl.p);
            float en = cells[i].energy - factor * (fr.e - fl.e);
            const float rho_air = cells[i].rho * cells[i].air - factor * (fr.y - fl.y);
            if (!std::isfinite(rho) || !std::isfinite(mom) || !std::isfinite(en) || !std::isfinite(rho_air)) {
                next_cells_[i] = cells[i];
                continue;
            }
            rho = std::max(rho, 1e-12f);
            const float damping = std::exp(-0.1f * damping_multiplier * dt /
                                           std::max(cells[i].rho, 1e-12f));
            mom *= damping;
            if (energy_relaxation_rate > 0.0f) {
                const float relax = std::exp(-energy_relaxation_rate * dt);
                en = ambient_energy + (en - ambient_energy) * relax;
            }
            en = std::max(en, 100.0f * pipe_gm1_inverse + 0.5f * mom * mom / rho);
            next_cells_[i] = {rho, mom, en, std::clamp(rho_air / rho, 0.0f, 1.0f)};
        }
        cells.swap(next_cells_);
    }
private:
    struct Flux { float m{}, p{}, e{}, y{}; };
    std::vector<GasCell> next_cells_;
    GasState flux_state(int idx) const {
        const auto& c = cells[static_cast<std::size_t>(idx)];
        const float rho = c.rho;
        const float inverse_rho = 1.0f / std::max(rho, 1e-12f);
        const float u = c.momentum * inverse_rho;
        const float p = pipe_gm1 * (c.energy - 0.5f * c.momentum * u);
        return {rho, u, p, std::clamp(c.air, 0.0f, 1.0f)};
    }
    static Flux flux(const GasState& a, const GasState& b) {
        const float ua = a.velocity, ub = b.velocity;
        const float inv_rho_a = 1.0f / std::max(a.rho, 1e-12f);
        const float inv_rho_b = 1.0f / std::max(b.rho, 1e-12f);
        const float ea = a.pressure * pipe_gm1_inverse * inv_rho_a + 0.5f * ua * ua;
        const float eb = b.pressure * pipe_gm1_inverse * inv_rho_b + 0.5f * ub * ub;
        const float c = std::max(std::abs(ua) + std::sqrt(gamma * a.pressure * inv_rho_a),
                                 std::abs(ub) + std::sqrt(gamma * b.pressure * inv_rho_b));
        return {
            0.5f * (a.rho * ua + b.rho * ub) - 0.5f * c * (b.rho - a.rho),
            0.5f * (a.rho * ua * ua + a.pressure + b.rho * ub * ub + b.pressure) - 0.5f * c * (b.rho * ub - a.rho * ua),
            0.5f * ((a.rho * ea + a.pressure) * ua + (b.rho * eb + b.pressure) * ub) - 0.5f * c * (b.rho * eb - a.rho * ea),
            0.5f * (a.rho * ua * a.air + b.rho * ub * b.air) - 0.5f * 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) {
        Orifice desc;
        desc.port = &port;
        desc.pipe = pipe;
        desc.left = left;
        desc.area = &area;
        desc.cd = cd;
        orifices_.push_back(desc);
    }
    void add_orifice_with_inertance(Port& port, int pipe, bool left, float& area,
                                    float cd, float effective_length, float loss_coefficient = 0.0f) {
        add_orifice(port, pipe, left, area, cd);
        auto& desc = orifices_.back();
        desc.use_inertance = true;
        desc.effective_length = effective_length;
        desc.loss_coefficient = loss_coefficient;
    }
    void add_closed_end(int pipe, bool left) {
        Orifice desc;
        desc.pipe = pipe;
        desc.left = left;
        desc.closed_end = true;
        orifices_.push_back(desc);
    }
    void add_open_end(int pipe, bool left, float pressure, float area) { opens_.push_back({pipe, left, pressure, area}); }
    float open_end_mass_flow(std::size_t index) const { return opens_.at(index).last_mass_flow; }
    float open_end_pressure(std::size_t index) const { return opens_.at(index).last_face_pressure; }
    void resolve(float dt) {
        for (auto& o : orifices_) resolve_orifice(o, dt);
        for (auto& o : opens_) resolve_open_end(o, dt);
    }
private:
    struct Orifice {
        Port* port{};
        int pipe{};
        bool left{};
        float* area{};
        float cd{1.0f};
        bool closed_end{}, use_inertance{};
        float effective_length{}, loss_coefficient{}, current_mass_flow{};
    };
    struct Open { int pipe; bool left; float pressure, area; float last_mass_flow{}, last_face_pressure{}; };
    struct OrificeSolution { float mass_flow{}, rho{}, velocity{}, pressure{}; };
    std::vector<Pipe1D>& pipes_;
    std::vector<Orifice> orifices_;
    std::vector<Open> opens_;

    static OrificeSolution isentropic(float pup, float rho_up, float t_up, float pdown, float area, float cd) {
        OrificeSolution out{0.0f, rho_up, 0.0f, pup};
        if (area <= 0.0f || pup <= 0.0f || rho_up <= 0.0f || t_up <= 0.0f) return out;
        float ratio = std::clamp(pdown / pup, 1e-6f, 1.0f);
        const float critical = std::pow(2.0f / (gamma + 1.0f), gamma / (gamma - 1.0f));
        ratio = std::max(ratio, critical);
        const float mach = std::sqrt(std::max(0.0f, 2.0f / (gamma - 1.0f) *
            (std::pow(ratio, -(gamma - 1.0f) / gamma) - 1.0f)));
        const float a_up = std::sqrt(gamma * gas_constant * t_up);
        out.velocity = mach * a_up;
        out.rho = rho_up * std::pow(ratio, 1.0f / gamma);
        out.pressure = pup * ratio;
        out.mass_flow = out.rho * out.velocity * area * cd;
        return out;
    }
    static float available_mass(const Port& port) {
        return (port.owner_air_mass ? *port.owner_air_mass : 0.0f) +
               (port.owner_exhaust_mass ? *port.owner_exhaust_mass : 0.0f);
    }
    void resolve_orifice(Orifice& o, float dt) {
        Pipe1D& pipe = pipes_.at(static_cast<std::size_t>(o.pipe));
        const int idx = o.left ? 0 : static_cast<int>(pipe.cells.size()) - 1;
        const GasState interior = pipe.state(idx);
        const float area = o.area ? std::max(*o.area, 0.0f) : 0.0f;
        if (o.closed_end || !o.port || area <= 1e-12f) {
            GasState wall = interior;
            wall.velocity = -interior.velocity;
            if (o.left) pipe.left_ghost = wall; else pipe.right_ghost = wall;
            if (o.port) o.port->mass_flow = 0.0f;
            return;
        }

        const float dp = o.port->pressure - interior.pressure;
        OrificeSolution solution{0.0f, interior.rho, 0.0f, interior.pressure};
        // Preserve the C# tie-break: at equal pressure it treats the volume as
        // upstream, retaining the volume-side density in the ghost state. That
        // matters when equal-pressure volume and pipe states have different
        // temperatures/densities, as on the initialized intake plenum.
        if (o.port->pressure >= interior.pressure) {
            solution = isentropic(o.port->pressure, o.port->density, o.port->temperature, interior.pressure, area, o.cd);
        } else {
            solution = isentropic(interior.pressure, interior.rho,
                interior.pressure / std::max(interior.rho * gas_constant, 1e-12f), o.port->pressure, area, o.cd);
            solution.mass_flow = -solution.mass_flow;
        }
        float mdot = solution.mass_flow; // positive means port component -> pipe
        const float adjacent_mass = pipe.cell_mass(idx);
        if (o.use_inertance) {
            const float rho_up = o.current_mass_flow >= 0.0f ? o.port->density : interior.rho;
            const float inertance = rho_up * o.effective_length / area;
            const float derivative = (dp - o.loss_coefficient * o.current_mass_flow) /
                                     std::max(inertance, 1e-12f);
            mdot = o.current_mass_flow + derivative * dt;
        }
        if (mdot > 0.0f && (o.port->owner_air_mass || o.port->owner_exhaust_mass))
            mdot = std::min(mdot, available_mass(*o.port) / std::max(dt, 1e-12f));
        if (mdot < 0.0f) mdot = std::max(mdot, -adjacent_mass / std::max(dt, 1e-12f));

        if (o.use_inertance) {
            o.current_mass_flow = mdot;
            const bool leaves_port = mdot >= 0.0f;
            const float rho_up = leaves_port ? o.port->density : interior.rho;
            const float temperature_up = leaves_port ? o.port->temperature :
                interior.pressure / std::max(interior.rho * gas_constant, 1e-12f);
            const float max_velocity = 0.9f * std::sqrt(gamma * gas_constant * temperature_up);
            const float max_flow = rho_up * max_velocity * area;
            mdot = std::clamp(mdot, -max_flow, max_flow);
            o.current_mass_flow = mdot;
            if (mdot >= 0.0f) {
                solution = isentropic(o.port->pressure, o.port->density, o.port->temperature,
                                      interior.pressure, area, o.cd);
            } else {
                solution = isentropic(interior.pressure, interior.rho, temperature_up,
                                      o.port->pressure, area, o.cd);
                solution.mass_flow = -solution.mass_flow;
            }
            solution.mass_flow = mdot;
            solution.rho = rho_up;
        }

        const bool from_port = mdot >= 0.0f;
        const float rho_face = std::max(solution.rho, 1e-9f);
        const float u = std::abs(mdot) / std::max(rho_face * area, 1e-12f);
        GasState ghost{rho_face, 0.0f, std::max(solution.pressure, 1.0f), from_port ? o.port->air_fraction : interior.air};
        ghost.velocity = ((from_port == o.left) ? u : -u);
        if (o.left) pipe.left_ghost = ghost; else pipe.right_ghost = ghost;
        o.port->mass_flow = -mdot;
        if (!from_port) o.port->air_fraction = interior.air;
        o.port->enthalpy = from_port
            ? o.port->enthalpy
            : gamma / (gamma - 1.0f) * interior.pressure / std::max(interior.rho, 1e-12f);
        (void)dp;
    }
    void resolve_open_end(Open& o, float dt) {
        Pipe1D& pipe = pipes_.at(static_cast<std::size_t>(o.pipe));
        const int idx = o.left ? 0 : static_cast<int>(pipe.cells.size()) - 1;
        const GasState interior = pipe.state(idx);
        const float c = std::sqrt(gamma * interior.pressure / std::max(interior.rho, 1e-12f));
        const float gm1 = gamma - 1.0f;
        const float j_plus = interior.velocity + 2.0f * c / gm1;
        const float j_minus = interior.velocity - 2.0f * c / gm1;
        const float entropy = interior.pressure / std::pow(std::max(interior.rho, 1e-12f), gamma);
        const float rho_iso = std::pow(std::max(o.pressure / entropy, 1e-12f), 1.0f / gamma);
        const float c_iso = std::sqrt(gamma * o.pressure / std::max(rho_iso, 1e-12f));
        const float u_iso = o.left ? j_minus + 2.0f * c_iso / gm1 : j_plus - 2.0f * c_iso / gm1;
        const bool supersonic = o.left ? (interior.velocity <= -c || u_iso <= -c_iso)
                                       : (interior.velocity >= c || u_iso >= c_iso);
        GasState ghost = supersonic ? interior : GasState{rho_iso, u_iso, o.pressure,
            ((o.left ? u_iso >= 0.0f : u_iso <= 0.0f) ? 1.0f : interior.air)};

        const float cell_mass = pipe.cell_mass(idx);
        float mdot_out = ghost.rho * ghost.velocity * o.area * (o.left ? -1.0f : 1.0f);
        if (mdot_out > 0.0f) mdot_out = std::min(mdot_out, cell_mass / std::max(dt, 1e-12f));
        else mdot_out = std::max(mdot_out, -10.0f * cell_mass / std::max(dt, 1e-12f));
        const float speed = std::abs(mdot_out) / std::max(ghost.rho * o.area, 1e-12f);
        ghost.velocity = o.left ? (mdot_out >= 0.0f ? -speed : speed)
                                : (mdot_out >= 0.0f ? speed : -speed);
        if (o.left) pipe.left_ghost = ghost; else pipe.right_ghost = ghost;
        o.last_mass_flow = mdot_out;
        o.last_face_pressure = ghost.pressure;
    }
};

} // namespace fluidsim::port