#pragma once
#include <algorithm>
#include <cmath>
#include <stdexcept>

namespace fluidsim::port {
// Lumped cantilever-petal model; geometry/material values require calibration.
class ReedValve {
public:
    float length{0.030f}, width{0.016f}, thickness{0.00035f};
    float youngs_modulus{70e9f}, material_density{1600.0f};
    float damping_ratio{0.30f}, maximum_lift{0.004f};
    int petals{4};
    float lift() const { return lift_; }
    float step(float dt, float pressure_difference) {
        if (!std::isfinite(dt) || dt <= 0 || !std::isfinite(pressure_difference))
            throw std::invalid_argument("invalid reed timestep or pressure");
        if (!(length > 0 && width > 0 && thickness > 0 && youngs_modulus > 0 && material_density > 0 &&
              maximum_lift > 0 && petals > 0 && damping_ratio >= 0))
            throw std::invalid_argument("invalid reed geometry/material");
        const float mass = .236f * material_density * length * width * thickness * petals;
        const float stiffness = petals * youngs_modulus * width * thickness * thickness * thickness /
                                (4 * length * length * length);
        const float damping = 2 * damping_ratio * std::sqrt(mass * stiffness);
        const float force = pressure_difference * (3.0f / 8.0f) * length * width * petals;
        const int steps = std::max(1, int(std::ceil(dt * std::sqrt(stiffness / mass) / .15f)));
        const float h = dt / steps;
        for (int i = 0; i < steps; i++) {
            velocity_ += (force - stiffness * lift_ - damping * velocity_) / mass * h;
            lift_ += velocity_ * h;
            if (lift_ < 0) {
                lift_ = 0;
                velocity_ = 0;
            }
            if (lift_ > maximum_lift) {
                lift_ = maximum_lift;
                velocity_ = 0;
            }
        }
        // Front-edge plus triangular side-gap areas for all petals.
        return petals * (width + length) * lift_;
    }

private:
    float lift_{}, velocity_{};
};
} // namespace fluidsim::port