2021-03-30 14:23:41 +08:00

104 lines
3.1 KiB
C#

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;
}
}
}