using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Assets.Scripts.Devices.Ant
{
    internal static class AirDensity
    {
        static double pi = 3.141592653589793;
        static double radPerDeg = 0;

        static double ae = 6378140.0;
        //static double alphae = 1 / 298.257;
        static double be = 0;
        static double phi0 = 0;
        static double R0 = 0;


        static double meanRadius = 6371004.0;  //all use SI units
        //static double semiMajorAxis = 6378140.0;
        static double gravityConst = 3.98600445e14;
        //static double J2 = 1.08263e-3;

        static AirDensity()
        {
            radPerDeg = pi / 180.0;
            be = 297.257 / 298.257 * ae;
            phi0 = 29.8336 * radPerDeg;
            R0 = ae * be / Math.Sqrt(ae * ae * Math.Sin(phi0) * Math.Sin(phi0) + be * be * Math.Cos(phi0) * Math.Cos(phi0));

        }


        /// <summary>
        /// 计算空气密度
        /// </summary>
        /// <param name="geoHeight"></param>
        /// <returns></returns>
        public static double GetAirDensity(double geoHeight)
        {
            //tested
            var H = geoHeightToPotentialH(geoHeight);
            var z = H / (1 - H / 6356766);
            var w = 0D;
            if (z < 11019.1)
            {
                w = 1.0 - H / 44330.8;
                return 1.225 * Math.Pow(w, 4.2559);
            }
            if (z < 20063.1)
            {
                w = Math.Exp((14964.7 - H) / 6341.6);
                return 1.225 * .15898 * w;
            }
            if (z < 32161.9)
            {
                w = 1.0 + (H - 24902.1) / 221552;
                return 1.225 * 3.2722e-2 * Math.Pow(w, -35.1629);
            }
            if (z < 47350.1)
            {
                w = 1.0 + (H - 39749.9) / 89410.7;
                return 1.225 * 3.2618e-3 * Math.Pow(w, -13.2011);
            }
            if (z < 51412.5)
            {
                w = Math.Exp((48625.2 - H) / 7922.3);
                return 1.225 * 9.4920e-4 * w;
            }
            if (z < 71802.0)
            {
                w = 1.0 - (H - 59439.0) / 88221.8;
                return 1.225 * 2.5280e-4 * Math.Pow(w, 11.2011);
            }
            if (z < 86000.0)
            {
                w = 1.0 - (H - 78030.3) / 100295.0;
                return 1.225 * 1.7632e-5 * Math.Pow(w, 16.0816);
            }
            if (z < 91000.0)
            {
                w = Math.Exp((87284.8 - H) / 5470.0);
                return 1.225 * 3.6411e-6 * w;
            }
            return .0;//linearInter(90000.0,13000.0,3.416e-6,8.152e-9,z);
        }


        private static double getRadiusFromLatitude(double latitude)
        {
            var s = Math.Sin(latitude);
            var c = Math.Cos(latitude);
            return (ae * be / Math.Pow((ae * ae * s * s + be * be * c * c), .5));
        }

        private static double geoHeightToPotentialH(double h)
        {
            return gravityConst * h / meanRadius / (meanRadius + h) / 9.8066;
        }
    }
}