using System;

namespace FluidSim.Components
{
    public static class NozzleFlow
    {
        /// <summary>
        /// Computes the nozzle‑exit primitive state and mass flow rate from a
        /// volume to a pipe, using isentropic relations. Follows ensim4's flow() logic.
        /// </summary>
        public static void Compute(double Pt_high, double Tt_high,
                                   double P_low, double gamma, double R, double area,
                                   out double rhoExit, out double uExit,
                                   out double pExit, out double mdot)
        {
            double gm1 = gamma - 1.0;
            double Pt_over_Ps = Pt_high / P_low;

            // Mach number (subsonic, clamped to 1)
            double M = Math.Sqrt(Math.Max(0.0,
                (2.0 / gm1) * (Math.Pow(Pt_over_Ps, gm1 / gamma) - 1.0)));
            if (M > 1.0) M = 1.0;

            double T_star = Tt_high / (1.0 + 0.5 * gm1 * M * M);
            double a_star = Math.Sqrt(gamma * R * T_star);
            double u_star = M * a_star;
            pExit = Pt_high * Math.Pow(1.0 + 0.5 * gm1 * M * M, -gamma / gm1);
            rhoExit = pExit / (R * T_star);
            uExit = u_star;   // positive away from high‑pressure side
            mdot = rhoExit * uExit * area;
        }

        /// <summary>
        /// Ambient cell for non‑reflecting open end (ensim4 calc_ambient_cell).
        /// </summary>
        public static void ComputeAmbientCell(double rhoInt, double uInt, double pInt,
                                               double pAmbient, double gamma,
                                               out double rhoAmb, out double uAmb,
                                               out double pAmb)
        {
            pAmb = pAmbient;
            uAmb = uInt;
            rhoAmb = rhoInt * Math.Pow(pAmb / pInt, 1.0 / gamma);
        }
    }
}