using System;

namespace FluidSim.Core
{
    /// <summary>
    /// Compressible flow through an orifice, modelled as an isentropic nozzle.
    /// Supports choked and unchoked flow, forward and reverse.
    /// </summary>
    public static class IsentropicOrifice
    {
        /// <summary>
        /// Compute mass flow and face primitive state for an orifice.
        /// </summary>
        /// <param name="pUp">Upstream stagnation pressure (Pa).</param>
        /// <param name="rhoUp">Upstream stagnation density (kg/m³).</param>
        /// <param name="gamma">Ratio of specific heats.</param>
        /// <param name="R">Specific gas constant (J/kg·K).</param>
        /// <param name="pDown">Downstream static pressure (Pa).</param>
        /// <param name="area">Effective orifice area (m²).</param>
        /// <param name="Cd">Discharge coefficient (default 0.62).</param>
        /// <param name="mdot">Mass flow rate (kg/s), positive from upstream to downstream.</param>
        /// <param name="rhoFace">Face density (kg/m³).</param>
        /// <param name="uFace">Face velocity (m/s).</param>
        /// <param name="pFace">Face pressure (Pa).</param>
        public static void Compute(double pUp, double rhoUp, double TUp, double gamma, double R,
                           double pDown, double area, double Cd,
                           out double mdot, out double rhoFace, out double uFace, out double pFace)
        {
            // mdot is positive from upstream to downstream.
            double pr = Math.Max(pDown / pUp, 1e-6);
            double prCrit = Math.Pow(2.0 / (gamma + 1.0), gamma / (gamma - 1.0));
            if (pr < prCrit) pr = prCrit;

            double M = Math.Sqrt((2.0 / (gamma - 1.0)) * (Math.Pow(pr, -(gamma - 1.0) / gamma) - 1.0));
            uFace = M * Math.Sqrt(gamma * R * TUp);
            rhoFace = rhoUp * Math.Pow(pr, 1.0 / gamma);
            pFace = pUp * pr;
            mdot = rhoFace * uFace * area * Cd;   // mass flow from upstream to downstream
        }
    }
}