#include "FluidSim/Solver.hpp"
#include <array>
#include <cmath>
#include <iostream>
#include <stdexcept>
using namespace fluidsim::port;
namespace {
void check(bool condition,const char* message) { if(!condition) throw std::runtime_error(message); }
struct Outcome { std::array<double,3> mass{}, pressure{}; std::uint64_t fallbacks{}; };
Outcome flow(bool reverse,bool rotate,float loss,bool mixed=false,bool swap_donors=false,bool legacy=false) {
    std::vector<Pipe1D> pipes;
    for(int j=0;j<3;++j) {
        pipes.emplace_back(.001f*(mixed?j+1:1),.2f,8);
        for(auto& cell:pipes.back().cells) {
            const double rho=mixed?(j==0?1.4:1.1):ambient_density;
            const double p=mixed?(j==0?120000:101325):101325;
            const double u=(swap_donors?(j==0?0:-50):(j==0?70:0))*(reverse?-1:1);
            cell={float(rho),float(rho*u),float(p/pipe_gm1+.5*rho*u*u),j==0?.7f:1.f,j==0?.1f:0.f};
            cell.fuel_equivalent=j==0?.15f:0.f;
        }
    }
    BoundarySystem boundaries(pipes);
    std::array<std::array<float,3>,3> axes{{{-1,0,0},{1,0,0},{0,1,0}}};
    if(rotate) for(auto& axis:axes) axis={-axis[1],axis[0],axis[2]};
    boundaries.add_junction({{0,reverse,axes[0]},{1,!reverse,axes[1]},{2,!reverse,axes[2]}},legacy?BoundarySystem::JunctionMethod::LegacyRusanov:BoundarySystem::JunctionMethod::Characteristic,loss);
    boundaries.add_closed_end(0,!reverse);boundaries.add_closed_end(1,reverse);boundaries.add_closed_end(2,reverse);
    const auto totals=[&] {
        std::array<double,5> t{};
        for(auto& p:pipes) for(int i=0;i<int(p.cells.size());++i) {
            const double m=p.cell_mass(i);
            t[0]+=m;t[1]+=double(p.cells[i].energy)*p.cell_area(i)*p.cell_length(i);
            t[2]+=m*p.cells[i].air;t[3]+=m*p.cells[i].fuel;
            t[4]+=m*p.cells[i].fuel_equivalent;
        }
        return t;
    };
    const auto initial=totals();Solver solver(pipes,boundaries);solver.minimum_substeps=1;
    // End before waves reflected by the sealed far ends reach the fitting.
    for(int i=0;i<30;++i) solver.step(1e-5f);
    const auto final=totals();
    for(int j=0;j<5;++j) check(std::abs(final[j]-initial[j])<std::max(1e-10,std::abs(initial[j])*1e-5),"junction lost mass, energy, or species (including burned-fuel origin)");
    Outcome result;
    for(int j=0;j<3;++j) {
        const int endpoint=(j==0?reverse:!reverse)?0:7;
        result.pressure[j]=pipes[j].state(endpoint).pressure;
        for(int i=0;i<8;++i) result.mass[j]+=pipes[j].cell_mass(i);
        check(pipes[j].positivity_repairs==0,"junction needed positivity repair");
    }
    result.fallbacks=boundaries.junction_fallbacks;
    check(result.fallbacks==0,"subsonic junction unexpectedly fell back");
    return result;
}
}
int main() {
    try {
        const auto ideal=flow(false,false,0);
        check(std::abs(ideal.mass[1]-ideal.mass[2])<1e-10,"ideal equal-area outlets are asymmetric");
        const auto directed=flow(false,false,1);
        check(directed.mass[1]>directed.mass[2]+1e-8,"turning loss does not favor straight-through flow");
        const auto rotated=flow(false,true,1), reversed=flow(true,false,1);
        for(int j=0;j<3;++j) {
            check(std::abs(rotated.mass[j]-directed.mass[j])<1e-10,"rotating junction changes flow");
            check(std::abs(reversed.mass[j]-directed.mass[j])<1e-9,"pipe coordinate reversal changes flow");
        }
        flow(false,false,.5f,true);flow(true,true,.5f,true);
        flow(false,false,.5f,true,true);flow(true,true,.5f,true,true);
        flow(false,false,0,false,false,true);flow(true,true,0,false,false,true);
        // Supersonic endpoint data must explicitly select the conservative fallback.
        std::vector<Pipe1D> pipes;for(int j=0;j<2;++j) pipes.emplace_back(.001f,.1f,4);
        for(auto& cell:pipes[0].cells) { const float u=500;cell.momentum=cell.rho*u;cell.energy+=.5f*cell.rho*u*u; }
        BoundarySystem boundaries(pipes);
        boundaries.add_junction({{0,false},{1,true}});
        boundaries.add_closed_end(0,true);boundaries.add_closed_end(1,false);
        Solver solver(pipes,boundaries);solver.minimum_substeps=1;solver.step(1e-6f);
        check(boundaries.junction_fallbacks==1,"supersonic fallback was not reported");
        bool rejected=false;
        try { BoundarySystem invalid(pipes);invalid.add_junction({{0,false},{1,true}},BoundarySystem::JunctionMethod::Characteristic,1); }
        catch(const std::invalid_argument&) {rejected=true;}
        check(rejected,"directionless turning-loss node was accepted");
        std::cout<<"Junction conservation, species, direction, rotation, reversal, and fallback tests passed\n";
    } catch(const std::exception& e) { std::cerr<<e.what()<<'\n';return 1; }
}