#include "FluidSim/Atmosphere.hpp"
#include "FluidSim/Crankcase.hpp"
#include "FluidSim/Fluid.hpp"
#include "FluidSim/Junction.hpp"
#include "FluidSim/Solver.hpp"
#include <cmath>
#include <complex>
#include <functional>
#include <iostream>
#include <string>
using namespace fluidsim::port;
namespace {
PipeSolver selected_solver=PipeSolver::Rusanov;
void check(bool b, const std::string &s) {
    if (!b)
        throw std::runtime_error(s);
}
void near(double a, double b, double tol, const char *s) {
    check(std::isfinite(a) && std::abs(a - b) <= tol, s);
}
void set(Pipe1D &p, int i, float pressure, float rho, float u = 0, float air = 1, float fuel = 0) {
    p.cells[i] = {rho, rho * u, pressure / pipe_gm1 + .5f * rho * u * u, air, fuel};
}
void lossless(Pipe1D &p) {
    p.solver=selected_solver;
    p.damping_multiplier = 0;
    p.energy_relaxation_rate = 0;
}
double mass(const Pipe1D &p) {
    double v = 0;
    for (int i = 0; i < int(p.cells.size()); i++)
        v += p.cell_mass(i);
    return v;
}
double energy(const Pipe1D &p) {
    double v = 0;
    for (int i = 0; i < int(p.cells.size()); i++)
        v += double(p.cells[i].energy) * p.cell_area(i) * p.cell_length(i);
    return v;
}
double vol_energy(const Volume0D &v) { return v.pressure() * v.volume() / (v.gamma - 1); }
} // namespace
int main(int argc,char** argv) {
    if(argc>1) selected_solver=parse_pipe_solver(argv[1]);
    int failed = 0, total = 0;
    auto test = [&](const char *n, std::function<void()> fn) {
        ++total;
        try {
            fn();
            std::cout << "PASS " << n << '\n';
        } catch (const std::exception &e) {
            ++failed;
            std::cout << "FAIL " << n << ": " << e.what() << '\n';
        }
    };
    test("solver flux consistency and reflection symmetry", [] {
        for(auto solver : {PipeSolver::Rusanov,PipeSolver::Hlle,PipeSolver::Hllc}) {
            GasState a{1.2f,30,101325,.8f,.05f}, b{.5f,-10,60000,.1f,.3f};
            const auto equal=Pipe1D::numerical_flux(a,a,solver);
            near(equal.m,a.rho*a.velocity,1e-5,"uniform mass flux");
            near(equal.p,a.rho*a.velocity*a.velocity+a.pressure,.02,"uniform momentum flux");
            const auto forward=Pipe1D::numerical_flux(a,b,solver);
            a.velocity=-a.velocity; b.velocity=-b.velocity;
            const auto reverse=Pipe1D::numerical_flux(b,a,solver);
            near(forward.m,-reverse.m,.0001,"orientation mass symmetry");
            near(forward.p,reverse.p,.1,"orientation momentum symmetry");
            near(forward.e,-reverse.e,1,"orientation energy symmetry");
            near(forward.y,-reverse.y,.0001,"orientation air symmetry");
            near(forward.fuel,-reverse.fuel,.0001,"orientation fuel symmetry");
        }
    });
    test("HLLC preserves stationary contact and falls back at expansion vacuum", [] {
        GasState a{1,0,1,.8f,.1f}, b{.125f,0,1,.1f,.5f};
        auto f=Pipe1D::numerical_flux(a,b,PipeSolver::Hllc);
        near(f.m,0,1e-7,"stationary contact mass diffusion");
        near(f.e,0,1e-7,"stationary contact energy diffusion");
        near(f.y,0,1e-7,"stationary contact tracer diffusion");
        a.velocity=-10; b.velocity=10;
        bool fallback=false;
        f=Pipe1D::numerical_flux(a,b,PipeSolver::Hllc,std::sqrt(1.4f),std::sqrt(1.4f/.125f),&fallback);
        check(fallback && std::isfinite(f.e),"strong expansion did not use safe HLLE flux");
    });
    test("unflanged radiation reflection magnitude and phase", [] {
        for (double ka : {.1, .5, 1., 2.}) {
            OpenEndRadiation reflection;
            std::complex<double> input{}, output{};
            double t = 0;
            const double dt = 2*3.141592653589793/(ka*4096);
            for (int i=0; i<4096*32; ++i) {
                t += dt;
                const double incident = std::cos(ka*t);
                const double reflected = reflection.process(incident, dt);
                if (i >= 4096*16) {
                    auto phase = std::exp(std::complex<double>(0, -ka*t));
                    input += incident*phase;
                    output += reflected*phase;
                }
            }
            constexpr double n=.167, d=n+2*.6133, e=(d*d-n*n-1)/2;
            const std::complex<double> z(0,ka);
            const auto expected = -(1.+n*z)/(1.+d*z+e*z*z);
            near(std::abs(output/input-expected), 0, .001, "radiation amplitude/phase mismatch");
            check(std::abs(output/input) < 1, "radiation boundary creates acoustic energy");
        }
    });
    test("N-port junction conserves sealed network mass and energy", [] {
        std::vector<Pipe1D> pipes;
        for(int i=0;i<3;++i) { pipes.emplace_back(.001f*(i+1),.3f,10);lossless(pipes.back()); }
        for(int i=0;i<10;++i) set(pipes[0],i,120000,1.4f);
        BoundarySystem boundaries(pipes);
        Junction junction(.0003f);
        const std::array<JunctionEnd,3> ends{{{0,false},{1,true},{2,true}}};
        junction.connect(boundaries,pipes,ends);
        boundaries.add_closed_end(0,true);boundaries.add_closed_end(1,false);boundaries.add_closed_end(2,false);
        Solver solver(pipes,boundaries);solver.add_component(junction);
        auto totals=[&] { std::array<double,2> t{junction.mass(),junction.internal_energy()};
            for(auto& p:pipes) { t[0]+=mass(p);t[1]+=energy(p); } return t; };
        const auto before=totals();
        for(int i=0;i<1000;++i) solver.step(1.f/44100);
        const auto after=totals();
        if(std::abs(after[0]-before[0])>before[0]*2e-6) std::cerr << "mass relative drift=" << (after[0]-before[0])/before[0] << "\n";
        near(after[0],before[0],before[0]*2e-6,"junction mass drift");
        near(after[1],before[1],before[1]*2e-6,"junction energy drift");
        check(pipes[1].state(0).pressure>101325,"junction did not transmit pressure to branch");
        for(auto& p:pipes) check(p.positivity_repairs==0,"junction required positivity repairs");
    });
    test("Zero-storage N-port junction conserves sealed network mass and energy", [] {
        std::vector<Pipe1D> pipes;
        for(int i=0;i<3;++i) { pipes.emplace_back(.001f*(i+1),.3f,10);lossless(pipes.back()); }
        for(int i=0;i<10;++i) set(pipes[0],i,120000,1.4f);
        BoundarySystem boundaries(pipes);
        boundaries.add_junction({{0,false},{1,true},{2,true}});
        boundaries.add_closed_end(0,true);boundaries.add_closed_end(1,false);boundaries.add_closed_end(2,false);
        Solver solver(pipes,boundaries);
        auto totals=[&] { std::array<double,2> t{0,0};
            for(auto& p:pipes) { t[0]+=mass(p);t[1]+=energy(p); } return t; };
        const auto before=totals();
        for(int i=0;i<1000;++i) solver.step(1.f/44100);
        const auto after=totals();
        // Float cell updates accumulate roundoff over 4000 substeps (50 ppm bound).
        near(after[0],before[0],before[0]*5e-5,"junction mass drift");
        near(after[1],before[1],before[1]*5e-5,"junction energy drift");
        check(pipes[1].state(0).pressure>101325,"junction did not transmit pressure to branch");
        for(auto& p:pipes) check(p.positivity_repairs==0,"junction required positivity repairs");
    });
    test("radiating acoustic volume source is independent of donor density", [] {
        for (bool left : {false,true}) {
            std::vector<Pipe1D> hot, cold;
            hot.emplace_back(.001f,.1f,1); cold.emplace_back(.001f,.1f,1);
            BoundarySystem a(hot), b(cold);
            a.add_radiating_open_end(0,left,101325,.001f);
            b.add_radiating_open_end(0,left,101325,.001f);
            for(int i=0;i<4096;++i) {
                const float velocity=2*std::sin(2*pi*150*i/44100);
                set(hot[0],0,101325,.5f,velocity);
                set(cold[0],0,101325,ambient_density,velocity);
                a.resolve(1.f/44100); b.resolve(1.f/44100);
                near(a.open_end_volume_flow(0),b.open_end_volume_flow(0),2e-7,
                     "temperature change contaminated acoustic volume source");
            }
        }
    });
    test("radiating mouth conserves reported mass at both orientations", [] {
        for (bool left : {false,true}) {
            std::vector<Pipe1D> p;
            p.emplace_back(.001f,.1f,10);
            lossless(p[0]);
            for (int i=0;i<10;++i) set(p[0],i,105000,1.2f);
            BoundarySystem b(p);
            b.add_radiating_open_end(0,left,101325,.001f);
            b.add_closed_end(0,!left);
            const double initial = mass(p[0]);
            b.resolve(1e-6f);
            const double exchange = -b.open_end_mass_flow(0)*1e-6;
            p[0].advance(1e-6f);
            near(mass(p[0])-initial,exchange,1e-11,"radiating mass flux mismatch");
        }
    });
    test("radiating mouth reflects rarefaction at either pipe orientation", [] {
        for (bool left : {false,true}) {
            std::vector<Pipe1D> p;
            p.emplace_back(.001f,1,400);
            lossless(p[0]);
            const double c=std::sqrt(1.4*101325/1.2);
            for(int i=0;i<400;++i) {
                double x=(i+.5)/400.;
                if(left) x=1-x;
                float dp=100*std::exp(-std::pow((x-.3)/.04,2));
                set(p[0],i,101325+dp,1.2f+dp/(c*c),(left?-1:1)*dp/(1.2*c));
            }
            BoundarySystem b(p);
            b.add_closed_end(0,!left);
            b.add_radiating_open_end(0,left,101325,.001f);
            Solver solver(p,b); solver.minimum_substeps=1;
            for(int i=0;i<int(1.4/c*44100);++i) solver.step(1.f/44100);
            check(p[0].state(left?280:119).pressure < 101305,"open mouth did not return an inverted pressure pulse");
            check(p[0].positivity_repairs==0,"radiation boundary required positivity repair");
        }
    });
    test("radiating inflow supplies fresh ambient stagnation enthalpy", [] {
        std::vector<Pipe1D> p;
        p.emplace_back(.001f,.1f,10);
        for(int i=0;i<10;++i) set(p[0],i,99000,1.2f,0,.1f,.2f);
        BoundarySystem b(p);
        b.add_radiating_open_end(0,true,101325,.001f,300);
        b.add_closed_end(0,false);
        b.resolve(1e-6f);
        const auto face=p[0].left_ghost;
        check(face.velocity>0 && face.air==1 && face.fuel==0,"backflow supplied exhaust or fuel");
        near(fluidsim::port::gamma/(fluidsim::port::gamma-1)*face.pressure/face.rho+.5*face.velocity*face.velocity,
             fluidsim::port::gamma/(fluidsim::port::gamma-1)*gas_constant*300, .1,"ambient enthalpy mismatch");
    });
    test("tracer underflow preserves bulk gas and normal concentrations", [] {
        for (float fraction : {1e-40f, 1e-30f, .5f}) {
            Pipe1D p(.001f, 1, 8);
            for (int i=0; i<8; ++i) set(p, i, 101325, 1.2f, 0, fraction);
            p.left_ghost = p.state(0);
            p.right_ghost = p.state(7);
            const double before_mass = mass(p), before_energy = energy(p);
            p.advance(1e-6f);
            near(mass(p), before_mass, 0, "tracer truncation changed total mass");
            near(energy(p), before_energy, 0, "tracer truncation changed energy");
            for (auto cell : p.cells) {
                check(std::fpclassify(cell.air) != FP_SUBNORMAL, "subnormal tracer survived");
                near(cell.air, fraction < std::numeric_limits<float>::min() ? 0 : fraction,
                     fraction * 1e-6, "normal tracer concentration changed");
            }
        }
    });
    test("invalid uniform geometry", [] {
        bool caught = false;
        try {
            Pipe1D p(0, .1f, 4);
        } catch (const std::invalid_argument &) {
            caught = true;
        }
        check(caught, "zero area accepted");
    });
    test("closed pipe conservation", [] {
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, 1, 50);
        lossless(p[0]);
        for (int i = 0; i < 50; i++)
            set(p[0], i, i < 25 ? 120000 : 100000, 1.2f);
        BoundarySystem b(p);
        b.add_closed_end(0, true);
        b.add_closed_end(0, false);
        Solver s(p, b);
        auto m = mass(p[0]), e = energy(p[0]);
        for (int i = 0; i < 2000; i++)
            s.step(1.f / 44100);
        near(mass(p[0]), m, m * 2e-5, "closed mass drift");
        near(energy(p[0]), e, e * 2e-5, "closed energy drift");
    });
    test("profiled pipe conservation", [] {
        std::vector<float> a(40), l(40, .025f);
        for (int i = 0; i < 40; i++)
            a[i] = .001f * (1 + .5f * std::sin(i * .1f));
        std::vector<Pipe1D> p;
        p.emplace_back(a, l);
        lossless(p[0]);
        for (int i = 0; i < 40; i++)
            set(p[0], i, i < 20 ? 120000 : 100000, 1.2f);
        BoundarySystem b(p);
        b.add_closed_end(0, true);
        b.add_closed_end(0, false);
        Solver s(p, b);
        auto m = mass(p[0]), e = energy(p[0]);
        for (int i = 0; i < 1000; i++)
            s.step(1.f / 44100);
        near(mass(p[0]), m, m * 3e-5, "profile mass drift");
        near(energy(p[0]), e, e * 3e-5, "profile energy drift");
    });
    test("volume pipe sealed conservation", [] {
        Volume0D v(.001f, 180000, 400);
        v.energy_relaxation_rate = 0;
        auto &port = v.create_port();
        std::vector<Pipe1D> p;
        p.emplace_back(.0003f, .5f, 25);
        lossless(p[0]);
        BoundarySystem b(p);
        float a = .0001f;
        b.add_orifice(port, 0, true, a);
        b.add_closed_end(0, false);
        Solver s(p, b);
        s.add_component(v);
        auto m = mass(p[0]) + v.mass(), e = energy(p[0]) + vol_energy(v);
        for (int i = 0; i < 2000; i++)
            s.step(1.f / 44100);
        near(mass(p[0]) + v.mass(), m, m * 3e-5, "coupling creates/deletes mass");
        near(energy(p[0]) + vol_energy(v), e, e * 3e-5, "coupling creates/deletes energy");
    });
    test("equal pressure hot contact", [] {
        Volume0D v(.001f, 101325, 600);
        v.energy_relaxation_rate = 0;
        auto &port = v.create_port();
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, .2f, 20);
        lossless(p[0]);
        BoundarySystem b(p);
        float a = .001f;
        b.add_orifice(port, 0, true, a);
        b.add_closed_end(0, false);
        Solver s(p, b);
        s.add_component(v);
        double m = mass(p[0]) + v.mass();
        for (int i = 0; i < 500; i++)
            s.step(1.f / 44100);
        near(mass(p[0]) + v.mass(), m, m * 2e-5, "contact produces mass");
        near(v.pressure(), 101325, 5, "stationary contact pressure changes");
    });
    test("NASA choked nozzle mass flow", [] {
        Atmosphere v;
        v.pressure = 200000;
        v.update_state(0);
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, .2f, 10);
        for (int i = 0; i < 10; i++)
            set(p[0], i, 10000, 10000 / (gas_constant * 300));
        BoundarySystem b(p);
        float a = .0001f;
        b.add_orifice(v.port, 0, true, a, .8f);
        b.resolve(1e-8f);
        double expected = .8 * a * 200000 / std::sqrt(300.) *
                          std::sqrt(fluidsim::port::gamma / gas_constant) *
                          std::pow(2 / (fluidsim::port::gamma + 1),
                                   (fluidsim::port::gamma + 1) / (2 * (fluidsim::port::gamma - 1)));
        near(-v.port.mass_flow, expected, expected * 2e-5,
             "choked mdot differs from isentropic sonic-throat formula");
    });
    test("inertance units", [] {
        Atmosphere v;
        v.pressure = 102325;
        v.update_state(0);
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, .1f, 10);
        BoundarySystem b(p);
        float a = .0001f;
        b.add_orifice_with_inertance(v.port, 0, true, a, 1, .05f);
        b.resolve(1e-6f);
        double expected = a / .05 * 1000 * 1e-6;
        near(-v.port.mass_flow, expected, expected * 1e-4,
             "mass-flow inertance must satisfy d(mdot)/dt=A*dp/L");
    });
    test("moving volume adiabatic law", [] {
        Volume0D v(.001f);
        v.energy_relaxation_rate = 0;
        v.volume_change_rate = .0001f;
        for (int i = 0; i < 1000; i++)
            v.update(1e-3f);
        near(v.volume(), .0011, 1e-8, "volume work did not change geometric volume");
        near(v.pressure(), 101325 * std::pow(.001 / .0011, 1.4), 3, "p V^gamma not conserved");
    });
    test("profile at rest is well balanced", [] {
        std::vector<float> a(80), l(80, .01f);
        for (int i = 0; i < 80; i++)
            a[i] = .001f * (1 + .7f * std::sin(.1f * i));
        std::vector<Pipe1D> p;
        p.emplace_back(a, l);
        lossless(p[0]);
        BoundarySystem b(p);
        b.add_closed_end(0, true);
        b.add_closed_end(0, false);
        Solver s(p, b);
        for (int j = 0; j < 1000; j++)
            s.step(1.f / 44100);
        for (int i = 0; i < 80; i++) {
            near(p[0].state(i).velocity, 0, 1e-3, "stationary nozzle accelerates");
            near(p[0].state(i).pressure, 101325, 2, "stationary nozzle pressure drifts");
        }
    });
    test("pipe connection all four orientations", [] {
        for (bool la : {false, true})
            for (bool lb : {false, true}) {
                std::vector<Pipe1D> p;
                p.emplace_back(.001f, .3f, 20);
                p.emplace_back(.001f, .3f, 20);
                for (auto &q : p)
                    lossless(q);
                for (int i = 0; i < 20; i++)
                    set(p[0], i, 130000, 1.4f, 0, .2f);
                BoundarySystem b(p);
                b.add_connection(0, la, 1, lb);
                b.add_closed_end(0, !la);
                b.add_closed_end(1, !lb);
                Solver s(p, b);
                double m = mass(p[0]) + mass(p[1]), e = energy(p[0]) + energy(p[1]);
                for (int i = 0; i < 2000; i++)
                    s.step(1.f / 44100);
                near(mass(p[0]) + mass(p[1]), m, m * 4e-5, "junction mass drift");
                near(energy(p[0]) + energy(p[1]), e, e * 4e-5, "junction energy drift");
                for (auto &q : p)
                    check(q.positivity_repairs == 0, "junction positivity repair");
            }
    });
    test("split pipe matches unsplit pipe", [] {
        std::vector<Pipe1D> a, b;
        a.emplace_back(.001f, 1, 80);
        b.emplace_back(.001f, .5f, 40);
        b.emplace_back(.001f, .5f, 40);
        lossless(a[0]);
        lossless(b[0]);
        lossless(b[1]);
        for (int i = 0; i < 80; i++) {
            float pr = 101325 + 100 * std::exp(-std::pow((i - 30) / 5.f, 2));
            set(a[0], i, pr, 1.2f);
            set(b[i / 40], i % 40, pr, 1.2f);
        }
        BoundarySystem ba(a), bb(b);
        ba.add_closed_end(0, true);
        ba.add_closed_end(0, false);
        bb.add_closed_end(0, true);
        bb.add_closed_end(1, false);
        bb.add_connection(0, false, 1, true);
        Solver sa(a, ba), sb(b, bb);
        for (int j = 0; j < 300; j++) {
            sa.step(1.f / 44100);
            sb.step(1.f / 44100);
        }
        for (int i = 0; i < 80; i++)
            near(a[0].state(i).pressure, b[i / 40].state(i % 40).pressure, .1,
                 "split junction changes solution");
    });
    test("three volumes shared donor conservation", [] {
        Volume0D a(.000001f, 1e6f, 500), b(.001f), c(.001f);
        a.energy_relaxation_rate = b.energy_relaxation_rate = c.energy_relaxation_rate = 0;
        auto &ap = a.create_port();
        auto &ap2 = a.create_port();
        auto &bp = b.create_port();
        auto &cp = c.create_port();
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, .1f, 10);
        p.emplace_back(.001f, .1f, 10);
        for (auto &q : p)
            lossless(q);
        BoundarySystem bs(p);
        float ar = .001f;
        bs.add_orifice(ap, 0, true, ar);
        bs.add_orifice(ap2, 1, true, ar);
        bs.add_orifice(bp, 0, false, ar);
        bs.add_orifice(cp, 1, false, ar);
        Solver s(p, bs);
        s.add_component(a);
        s.add_component(b);
        s.add_component(c);
        double m = mass(p[0]) + mass(p[1]) + a.mass() + b.mass() + c.mass();
        double e = energy(p[0]) + energy(p[1]) + vol_energy(a) + vol_energy(b) + vol_energy(c);
        for (int j = 0; j < 3000; j++)
            s.step(1.f / 44100);
        near(mass(p[0]) + mass(p[1]) + a.mass() + b.mass() + c.mass(), m, m * 5e-5,
             "shared donor mass drift");
        near(energy(p[0]) + energy(p[1]) + vol_energy(a) + vol_energy(b) + vol_energy(c), e, e * 5e-5,
             "shared donor energy drift");
        check(a.pressure() > 0, "donor vacuum");
    });
    test("open boundary reports actual conserved flux both ends", [] {
        for (bool left : {false, true}) {
            std::vector<Pipe1D> p;
            p.emplace_back(.001f, .1f, 10);
            lossless(p[0]);
            for (int i = 0; i < 10; i++)
                set(p[0], i, 105000, 1.2f);
            BoundarySystem b(p);
            b.add_open_end(0, left, 101325, .001f);
            b.add_closed_end(0, !left);
            double m = mass(p[0]);
            b.resolve(1e-6f);
            double dm = -b.open_end_mass_flow(0) * 1e-6;
            p[0].advance(1e-6f);
            near(mass(p[0]) - m, dm, 1e-11,
                 "reported radiation flow differs from finite-volume mass exchange");
        }
    });
    test("duplicate and invalid boundary rejection", [] {
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, 1, 10);
        BoundarySystem b(p);
        b.add_closed_end(0, true);
        bool caught = false;
        try {
            b.add_open_end(0, true, 101325, .001f);
        } catch (const std::invalid_argument &) {
            caught = true;
        }
        check(caught, "duplicate face accepted");
        caught = false;
        try {
            b.add_closed_end(9, true);
        } catch (const std::invalid_argument &) {
            caught = true;
        }
        check(caught, "invalid pipe accepted");
    });
    test("supersonic outlet extrapolates", [] {
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, 1, 10);
        for (int i = 0; i < 10; i++)
            set(p[0], i, 101325, 1.2f, 700);
        BoundarySystem b(p);
        b.add_open_end(0, false, 50000, .001f);
        b.resolve(1e-6f);
        near(b.open_end_mass_flow(0), .001 * 1.2 * 700, 1e-5,
             "supersonic outlet responds to downstream pressure");
    });
    test("underdefined supersonic inlet rejected", [] {
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, 1, 10);
        set(p[0], 0, 101325, 1.2f, 700);
        BoundarySystem b(p);
        b.add_open_end(0, true, 101325, .001f);
        bool caught = false;
        try {
            b.resolve(1e-6f);
        } catch (const std::runtime_error &) {
            caught = true;
        }
        check(caught, "supersonic inlet silently accepted pressure-only data");
    });
    test("closed and open acoustic reflection signs", [] {
        for (bool closed : {false, true}) {
            std::vector<Pipe1D> p;
            p.emplace_back(.001f, 1, 400);
            lossless(p[0]);
            const double c = std::sqrt(1.4 * 101325 / 1.2);
            for (int i = 0; i < 400; i++) {
                double x = (i + .5) / 400.;
                float dp = 100 * std::exp(-std::pow((x - .3) / .04, 2));
                set(p[0], i, 101325 + dp, 1.2f + dp / (c * c), dp / (1.2 * c));
            }
            BoundarySystem b(p);
            b.add_closed_end(0, true);
            if (closed)
                b.add_closed_end(0, false);
            else
                b.add_open_end(0, false, 101325, .001f);
            Solver s(p, b);
            s.minimum_substeps = 1;
            const int steps = int((1.4 / c) * 44100);
            for (int j = 0; j < steps; j++)
                s.step(1.f / 44100);
            float dp = p[0].state(119).pressure - 101325;
            check(closed ? dp > 30 : dp < -30, "wrong reflected pressure sign/amplitude");
            check(p[0].positivity_repairs == 0, "acoustic positivity repair");
        }
    });
    test("Sod shock tube convergence", [] {
        double previous = 1e9;
        for (int n : {100, 200, 400}) {
            Pipe1D p(1, 1, n);
            lossless(p);
            for (int i = 0; i < n; i++)
                set(p, i, i < n / 2 ? 1 : .1f, i < n / 2 ? 1 : .125f);
            double t = 0;
            while (t < .2) {
                p.left_ghost = p.state(0);
                p.right_ghost = p.state(n - 1);
                float dt = std::min(.2 - t, double(p.stable_timestep(.7f)));
                p.advance(dt);
                t += dt;
            }
            double error = 0;
            for (int i = 0; i < n; i++) {
                double xi = ((i + .5) / n - .5) / .2;
                double rho;
                if (xi < -std::sqrt(1.4))
                    rho = 1;
                else if (xi < -.0702728)
                    rho = std::pow((2 / (1.4 + 1)) + (1.4 - 1) / (1.4 + 1) * (-xi / std::sqrt(1.4)),
                                   2 / (1.4 - 1));
                else if (xi < .9274526)
                    rho = .4263194;
                else if (xi < 1.7521557)
                    rho = .2655737;
                else
                    rho = .125;
                error += std::abs(p.state(i).rho - rho) / n;
            }
            std::cout << "  Sod n=" << n << " density L1=" << error << '\n';
            check(error < previous * .9, "shock density error does not decrease on refinement");
            previous = error;
            check(p.positivity_repairs == 0, "Sod positivity repair");
        }
        check(previous < .025, "Sod density error excessive");
    });
    test("strong pulse and near-vacuum rarefaction", [] {
        for (bool vacuum : {false, true}) {
            Pipe1D p(1, 1, 200);
            lossless(p);
            for (int i = 0; i < 200; i++)
                set(p, i, vacuum ? .4f : (i < 100 ? 1000 : .01f), 1, vacuum ? (i < 100 ? -2 : 2) : 0);
            double t = 0;
            while (t < .005) {
                p.left_ghost = p.state(0);
                p.right_ghost = p.state(199);
                float dt = std::min(.005 - t, double(p.stable_timestep(.5f)));
                p.advance(dt);
                t += dt;
            }
            for (int i = 0; i < 200; i++)
                check(p.state(i).rho > 0 && p.state(i).pressure > 0, "negative shock/vacuum state");
            check(p.positivity_repairs == 0, "stress case required conservation-breaking positivity repair");
        }
    });

    test("small-signal Helmholtz resonance", [] {
        std::vector<Pipe1D> p;
        float area = .000314159f;
        p.emplace_back(area, .05f, 40);
        Volume0D v(.001f, 101335, 300);
        auto &port = v.create_port();
        BoundarySystem b(p);
        b.add_orifice(port, 0, true, area, 1);
        b.add_open_end(0, false, 101325, area);
        Solver solver(p, b);
        solver.minimum_substeps = 1;
        solver.add_component(v);
        double prev = v.pressure() - 101325, first = 0, last = 0;
        int count = 0;
        for (int i = 0; i < 4410; i++) {
            solver.step(1.f / 44100);
            double now = v.pressure() - 101325;
            check(std::abs(now) < 15, "small-signal cavity grows instead of oscillating/decaying");
            if (prev > 0 && now <= 0 && i > 441) {
                double t = (i + prev / (prev - now)) / 44100;
                if (!count)
                    first = t;
                last = t;
                count++;
            }
            prev = now;
        }
        double frequency = (count - 1) / (last - first);
        double expected =
            std::sqrt(1.4 * 287 * 300) / (2 * 3.141592653589793) * std::sqrt(area / (.001 * .05));
        std::cout << "  Helmholtz frequency=" << frequency << " analytic lumped=" << expected << '\n';
        check(count > 5, "not enough cavity oscillations");
        near(frequency, expected, expected * .05, "Helmholtz mode frequency error exceeds 5%");
        check(p[0].positivity_repairs == 0, "cavity requires positivity repairs");
    });
    test("invalid volume and timestep rejection", [] {
        bool caught = false;
        try {
            Volume0D v(-1);
        } catch (const std::invalid_argument &) {
            caught = true;
        }
        check(caught, "negative volume accepted");
        Pipe1D p(.001f, .1f, 10);
        caught = false;
        try {
            p.required_substeps(std::nanf(""));
        } catch (const std::invalid_argument &) {
            caught = true;
        }
        check(caught, "NaN timestep accepted");
    });
    test("fresh-air species conservation across connection", [] {
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, .5f, 40);
        p.emplace_back(.001f, .5f, 40);
        for (int k = 0; k < 2; k++) {
            lossless(p[k]);
            for (int i = 0; i < 40; i++)
                set(p[k], i, k ? 100000 : 120000, 1.2f, 0, k ? .9f : .1f,
                    k ? .02f : .05f);
        }
        auto air_mass = [&] {
            double sum = 0;
            for (auto &q : p)
                for (int i = 0; i < 40; i++)
                    sum += q.cell_mass(i) * q.cells[i].air;
            return sum;
        };
        auto fuel_mass = [&] {
            double sum = 0;
            for (auto &q : p)
                for (int i = 0; i < 40; i++)
                    sum += q.cell_mass(i) * q.cells[i].fuel;
            return sum;
        };
        double before = air_mass();
        double fuel_before = fuel_mass();
        BoundarySystem b(p);
        b.add_closed_end(0, true);
        b.add_closed_end(1, false);
        b.add_connection(0, false, 1, true);
        Solver solver(p, b);
        for (int i = 0; i < 3000; i++)
            solver.step(1.f / 44100);
        near(air_mass(), before, before * 4e-5, "species created at shared face");
        near(fuel_mass(), fuel_before, fuel_before * 4e-5, "fuel species created at shared face");
    });
    test("unsupported unequal junction rejected", [] {
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, .5f, 20);
        p.emplace_back(.002f, .5f, 20);
        BoundarySystem b(p);
        bool caught = false;
        try {
            b.add_connection(0, false, 1, true);
        } catch (const std::invalid_argument &) {
            caught = true;
        }
        check(caught, "unvalidated sudden area junction silently accepted");
    });
    test("unconnected endpoint rejected by solver", [] {
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, .5f, 20);
        BoundarySystem b(p);
        b.add_closed_end(0, true);
        Solver solver(p, b);
        bool caught = false;
        try {
            solver.step(1.f / 44100);
        } catch (const std::runtime_error &) {
            caught = true;
        }
        check(caught, "unconnected face silently used default ghost");
    });
    test("closed duct acoustic eigenfrequency", [] {
        for (int n : {50, 100, 200}) {
            std::vector<Pipe1D> p;
            p.emplace_back(.001f, 1, n);
            double c = std::sqrt(1.4 * 101325 / 1.2);
            for (int i = 0; i < n; i++) {
                double dp = 100 * std::cos(3.141592653589793 * (i + .5) / n);
                set(p[0], i, 101325 + dp, 1.2f + dp / (c * c));
            }
            BoundarySystem b(p);
            b.add_closed_end(0, true);
            b.add_closed_end(0, false);
            Solver solver(p, b);
            solver.minimum_substeps = 1;
            double prev = 1, first = 0, last = 0;
            int count = 0;
            for (int j = 0; j < 1764; j++) {
                solver.step(1.f / 44100);
                double projection = 0;
                for (int i = 0; i < n; i++)
                    projection +=
                        (p[0].state(i).pressure - 101325) * std::cos(3.141592653589793 * (i + .5) / n);
                if (prev > 0 && projection <= 0) {
                    double t = (j + prev / (prev - projection)) / 44100;
                    if (!count)
                        first = t;
                    last = t;
                    count++;
                }
                prev = projection;
            }
            double f = (count - 1) / (last - first);
            std::cout << "  closed mode n=" << n << " f=" << f << " analytic=" << c / 2 << '\n';
            check(count > 4, "insufficient acoustic cycles");
            near(f, c / 2, c / 2 * .01, "closed duct frequency exceeds 1% error");
        }
    });
    test("inertance impulse and dissipative resistance", [] {
        Atmosphere v;
        v.pressure = 102325;
        v.update_state(0);
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, .1f, 10);
        BoundarySystem b(p);
        float area = .0001f;
        b.add_orifice_with_inertance(v.port, 0, true, area, 1, .05f, 1);
        for (int j = 0; j < 10000; j++)
            b.resolve(.0001f);
        double equilibrium = std::sqrt(2 * v.density() * area * area * 1000);
        near(-v.port.mass_flow, equilibrium, equilibrium * .001,
             "inertance resistance steady flow has wrong units or sign");
    });
    test("incompatible gas coupling rejected", [] {
        Volume0D v(.001f, 101325, 300, 300, 1.3f);
        auto &port = v.create_port();
        std::vector<Pipe1D> p;
        p.emplace_back(.001f, .1f, 10);
        BoundarySystem b(p);
        float area = .0001f;
        bool caught = false;
        try {
            b.add_orifice(port, 0, true, area);
        } catch (const std::invalid_argument &) {
            caught = true;
        }
        check(caught, "different gamma gases silently coupled");
    });
    test("thermal relaxation matches analytic energy decay", [] {
        Volume0D v(.001f, 101325, 600);
        v.energy_relaxation_rate = 2;
        double initial = vol_energy(v), target = v.mass() * 287 * 300 / (1.4 - 1);
        for (int i = 0; i < 1000; i++)
            v.update(.001f);
        near(vol_energy(v), target + (initial - target) * std::exp(-2), initial * 3e-5,
             "thermal relaxation differs from exact exponential");
    });

    test("sealed crankcase compression conserves mass above five bar", [] {
        Crankshaft crank;
        crank.angle = 0;
        Crankcase v(crank, .02725f, .11f, .054f, 1e-5f);
        double m = v.mass(), invariant = v.pressure() * std::pow(v.volume(), 1.4);
        for (int i = 1; i <= 1000; i++) {
            crank.angle = pi * i / 1000;
            v.pre_step(1e-5f);
            v.update_state(1e-5f);
        }
        near(v.mass(), m, m * 1e-6, "crankcase pressure cap deletes mass");
        near(v.pressure() * std::pow(v.volume(), 1.4), invariant, invariant * 1e-4,
             "crankcase compression violates adiabatic invariant");
        check(v.pressure() > 5e5, "compression did not exercise old hard pressure cap");
    });
    std::cout << total - failed << "/" << total << " passed\n";
    return failed ? 1 : 0;
}