#include "FluidSim/Carburetor.hpp"

#include <algorithm>
#include <cmath>

namespace fluidsim::port {
namespace {
constexpr float pi = 3.14159265358979323846f;
float circle_area(float diameter) { return .25f * pi * diameter * diameter; }
float smoothstep(float value) {
    value = std::clamp(value, 0.0f, 1.0f);
    return value * value * (3.0f - 2.0f * value);
}
}

Carburetor::Carburetor(float throat_diameter)
    : bowl_mass_(bowl_capacity), throat_area_(circle_area(std::max(throat_diameter, 1e-4f))) {}

void Carburetor::step(float dt, float pedal, float slide, float rpm, float air_mass_flow,
                      float air_density, float upstream_pressure, float downstream_pressure) {
    fault_ = !std::isfinite(dt) || dt <= 0.0f || !std::isfinite(pedal) ||
             !std::isfinite(slide) || !std::isfinite(rpm) || rpm < 0.0f ||
             !std::isfinite(air_mass_flow) || !std::isfinite(air_density) || air_density <= 0.0f ||
             !std::isfinite(upstream_pressure) ||
             !std::isfinite(downstream_pressure) || upstream_pressure <= 0.0f ||
             downstream_pressure <= 0.0f;
    if (fault_) {
        delivered_fuel_flow_ = 0.0f;
        overrun_ = false;
        return;
    }

    pedal = std::clamp(pedal, 0.0f, 1.0f);
    slide = std::clamp(slide, 0.0f, 1.0f);
    air_mass_flow = std::max(air_mass_flow, 0.0f);
    if (!initialized_) {
        previous_pedal_ = pedal;
        initialized_ = true;
    }
    const bool throttle_lift = previous_pedal_ > 0.08f && pedal <= 0.005f && rpm > overrun_rpm;
    if (throttle_lift)
        lift_enrichment_ = 1.0f;
    previous_pedal_ = pedal;
    lift_enrichment_ *= std::exp(-dt / 0.32f);
    if (pedal > 0.005f || rpm < overrun_exit_rpm)
        overrun_ = false;
    else if (throttle_lift)
        overrun_ = true;

    const float air_alpha = -std::expm1(-dt / 0.012f);
    filtered_air_flow_ += (air_mass_flow - filtered_air_flow_) * air_alpha;

    // The open slide controls the effective venturi area. Its local dynamic
    // pressure and the measured static-pressure loss provide the jet signal.
    const float venturi_area = throat_area_ * 0.55f * (0.10f + 0.90f * slide);
    const float velocity = filtered_air_flow_ / std::max(air_density * venturi_area, 1e-9f);
    const float dynamic_pressure = 0.5f * air_density * velocity * velocity;
    const float static_drop = std::max(upstream_pressure - downstream_pressure, 0.0f);
    venturi_depression_ = std::clamp(dynamic_pressure + 0.18f * static_drop, 0.0f, 60000.0f);

    const float pilot_area = circle_area(pilot_jet_diameter);
    const float main_area = circle_area(main_jet_diameter);
    const float needle_jet_area = circle_area(needle_jet_diameter);
    const float needle_lift = smoothstep((slide - 0.04f) / 0.86f);
    const float needle_area = needle_jet_area * needle_lift;
    const float main_path_area = needle_area > 1e-12f ?
        1.0f / std::sqrt(1.0f / (main_area * main_area) + 1.0f / (needle_area * needle_area)) : 0.0f;
    const float pilot_depression = bowl_head_pressure + (0.62f * static_drop +
                                   0.25f * dynamic_pressure) * (1.0f + 0.25f * lift_enrichment_);
    const float main_depression = bowl_head_pressure + venturi_depression_;
    const float pilot_flow = pilot_discharge_coefficient * pilot_area *
        std::sqrt(2.0f * fuel_density * std::max(pilot_depression, 0.0f));
    const float main_flow = main_discharge_coefficient * main_path_area *
        std::sqrt(2.0f * fuel_density * std::max(main_depression, 0.0f));
    float demanded_fuel_flow = pilot_flow + main_flow;

    const float deposit_fraction = std::clamp(film_fraction * (0.55f + 0.45f * slide), 0.0f, 0.8f);
    const float deposited = demanded_fuel_flow * deposit_fraction;
    const float evaporation_time = film_time_constant * (overrun_ ? 0.58f : 1.0f);
    const float evaporated = std::min(film_mass_ / dt,
                                      film_mass_ / std::max(evaporation_time, 1e-4f));
    float jet_flow = demanded_fuel_flow;
    const float available_flow = bowl_mass_ / dt;
    if (jet_flow > available_flow) {
        const float scale = available_flow / std::max(jet_flow, 1e-12f);
        jet_flow *= scale;
    }
    bowl_mass_ -= jet_flow * dt;
    bowl_mass_ += std::min(bowl_refill_rate * dt, bowl_capacity - bowl_mass_);
    bowl_mass_ = std::clamp(bowl_mass_, 0.0f, bowl_capacity);

    // Deposition cannot exceed the fuel that actually left the bowl.
    const float actual_deposit = std::min(deposited, jet_flow);
    film_mass_ = std::max(film_mass_ + (actual_deposit - evaporated) * dt, 0.0f);
    delivered_fuel_flow_ = std::max(jet_flow - actual_deposit + evaporated, 0.0f);
    const float fuel_alpha = -std::expm1(-dt / 0.025f);
    filtered_fuel_flow_ += (delivered_fuel_flow_ - filtered_fuel_flow_) * fuel_alpha;
    if (filtered_air_flow_ > 1e-7f && filtered_fuel_flow_ > 1e-9f)
        indicated_afr_ = std::clamp(filtered_air_flow_ / filtered_fuel_flow_, 4.0f, 40.0f);
}

} // namespace fluidsim::port