using System;
using Orbiter.Interfaces;
using Orbiter.Wrapper;
using Vessel=Orbiter.Interfaces.Vessel;

namespace Orbiter.StockVessels
{
    public class Ramjet
    {
        private Vessel me;
        private uint nthdef=0;
        
        public Ramjet(Vessel vessel)
        {
            me = vessel;
        }
        
        private struct THDEF
        {
            //Static parameters
            
            /// <summary>
            /// thruster handle
            /// </summary>
		    public IntPtr th;
            /// <summary>
            /// fuel heating parameter [J/kg]
            /// </summary>
		    public double Qr;
            /// <summary>
            /// air intake cross section [m^2]
            /// </summary>
		    public double Ai;
            /// <summary>
            /// max. burner temperature [K]
            /// </summary>
		    public double Tb_max;
            /// <summary>
            /// max. fuel flow rate [kg/s]
            /// </summary>
		    public double dmf_max;
            
            //Dynamic parameters
            
            /// <summary>
            /// current fuel mass rate [kg/s]
            /// </summary>
		    public double dmf;
            /// <summary>
            /// current thrust [N]
            /// </summary>
		    public double F;
            /// <summary>
            /// temperatures
            /// </summary>
		    public double[] T;

            public THDEF(IntPtr th, params double[] values)
            {
                this.th = th;
                if (values.Length > 0) Qr = values[0]; else Qr = 0;
                if (values.Length > 1) Ai = values[1]; else Ai = 0;
                if (values.Length > 2) Tb_max = values[2]; else Tb_max = 0;
                if (values.Length > 3) dmf_max = values[3]; else dmf_max = 0;
                if (values.Length > 4) dmf = values[4]; else dmf = 0;
                if (values.Length > 5) F = values[5]; else F = 0;
                T=new double[3];
                if (values.Length > 6) T[0] = values[6]; else T[0] = 0;
                if (values.Length > 7) T[1] = values[7]; else T[1] = 0;
                if (values.Length > 8) T[2] = values[8]; else T[2] = 0;
            }
	    }
        
        private THDEF[] thdef;
	    
        public void AddThrusterDefinition (IntPtr th,double Qr, double Ai, double Tb_max, double dmf_max)
        {
            THDEF[] tmp=new THDEF[nthdef+1];
            thdef.CopyTo(tmp, 0);
            thdef = tmp;
            thdef[nthdef++] = new THDEF(th, Qr, Ai, Tb_max, dmf_max);
        }
        
        public void Thrust (double[] F)
        {
            IntPtr hBody = This.AtmosphereReference;
            
            if (hBody!=IntPtr.Zero)
            {
                // atmospheric parameters available
                AtmosphereConstants atm = OAPI.GetPlanetAtmConstants(hBody);
                double M, Fs, T0, Td, Tb, Tb0, Te, p0, pd, D, rho, cp, v0, ve, tr, lvl, dma, dmf, precov, dmafac;
                const double eps = 1e-4;
                const double dma_scale = 2.7e-4;

                M = This.MachNumber;                     // Mach number
                T0 = This.AtmosphereTemperature;                 // freestream temperature
                p0 = This.AtmospherePressure;                    // freestream pressure
                rho = This.AtmosphereDensity;                     // freestream density
                cp = atm.Gamma * atm.R / (atm.Gamma - 1.0);      // specific heat (pressure)
                v0 = M * Math.Sqrt(atm.Gamma * atm.R * T0);         // freestream velocity
                tr = (1.0 + 0.5 * (atm.Gamma - 1.0) * M * M);          // temperature ratio
                Td = T0 * tr;                                     // diffuser temperature
                pd = p0 * Math.Pow(Td / T0, atm.Gamma / (atm.Gamma - 1.0)); // diffuser pressure
                precov = Math.Max(0.0, 1.0 - 0.075 * Math.Pow(Math.Max(M, 1.0) - 1.0, 1.35)); // pressure recovery
                dmafac = dma_scale * precov * pd;

                for (uint i = 0; i < nthdef; i++)
                {
                    Tb0 = thdef[i].Tb_max;                        // max burner temperature
                    if (Tb0 > Td)
                    {                                // we are within operational range
                        lvl = This.GetThrusterLevel(thdef[i].th); // throttle level
                        D = (Tb0 - Td) / (thdef[i].Qr / cp - Tb0); // max fuel-to-air ratio (what if negative?)
                        D *= lvl;                                // actual fuel-to-air ratio

                        dma = dmafac * thdef[i].Ai;               // air mass flow rate [kg/s]
                        //dma  = rho * v0 * thdef[i]->Ai;            // air mass flow rate
                        dmf = D * dma;                            // fuel mass flow rate
                        if (dmf > thdef[i].dmf_max)
                        {             // max fuel rate exceeded
                            dmf = thdef[i].dmf_max;
                            D = dmf / dma;
                        }
                        Tb = (D * thdef[i].Qr / cp + Td) / (1.0 + D); // actual burner temperature
                        Te = Tb * Math.Pow(p0 / pd, (atm.Gamma - 1.0) / atm.Gamma); // exhaust temperature
                        ve = Math.Sqrt(2.0 * cp * (Tb - Te));              // exhaust velocity
                        Fs = (1.0 + D) * ve - v0;                     // specific thrust
                        thdef[i].F = F[i] = Math.Max(0.0, Fs * dma);    // thrust force
                        thdef[i].dmf = dmf;
                        thdef[i].T[1] = Tb;
                        thdef[i].T[2] = Te;

                    }
                    else
                    {                                       // overheating!

                        thdef[i].F = F[i] = 0.0;
                        thdef[i].dmf = 0.0;
                        thdef[i].T[1] = thdef[i].T[2] = Td;

                    }
                    thdef[i].T[0] = Td;
                }

            }
            else
            {
                // no atmospheric parameters
                for (uint i = 0; i < nthdef; i++)
                {
                    thdef[i].F = F[i] = 0.0;
                    thdef[i].dmf = 0.0;
                }
            }
        }

        public double TSFC(uint idx)
        {
            const double eps = 1e-5;
            return thdef[idx].dmf / (thdef[idx].F + eps);
        }
        
        public double DMF (uint idx)
        {
            return thdef[idx].dmf;
        }
        
        public double Temp (uint idx, uint which)
        {
            return thdef[idx].T[which];
        }

        public Vessel This
        {
            get { return me; }
        }
    }
}
