SimulationServer/Models/SimulationCommon/EquationHelper.cs

using BHJD.DEMdll.Entity;
using BHJD.DEMdll.Public;
using System.Diagnostics;
using static BHJD.DEMdll.Public.IHttpHelper;
using Unity.Mathematics;
using Newtonsoft.Json;
using Newtonsoft.Json.Linq;
using SJ;
using Model;
using MathNet.Numerics;

namespace MuShiApp
{
    public class EquationHelper
    {
        //数据库访问类，dll有
        //IDbHelper _db = null;
        //http请求类，dll有
        IHttpHelper m_HttpHelper = null;

        public string baseUrl;

        public double x;
        public EquationHelper(string baseUrl)
        {
            this.baseUrl = baseUrl;
            m_HttpHelper = new HttpHelper();
            //this._db = db;
            //m_HttpHelper = Factory.Load();
        }


        //a、b、c、d、r、e均为修正系数
        double[] a = { 9.5128, 1.3738, 0.0875 };
        double[] b = { -49.6937, -13.2945, -2.5112 };
        double[] c = { 31.6031, 13.1456, 6.1604 };
        double[] d = { 95.0664, 97.0952, 94.9145 };

        double[][] r = new double[3][];
        double[][] e = new double[3][];
        double[] posibility = { 0, 0, 0 };

        public double GetMushiSeaProbability(double3 aircraftPoint, double visibility, double wavehigh, TargetPoint targetPoint)
        {
            r[0] = new double[] { 4.375, 1.18, 1.01, 1 };
            r[1] = new double[] { 2.85, 1.51, 1.18, 1 };
            r[2] = new double[] { 2.75, 2.1, 1.2, 1 };
            e[0] = new double[] { -10, -30, -50, -65, -80, -90 };
            e[1] = new double[] { -2, -15, -30, -50, -70, -80 };
            e[2] = new double[] { 0, 0, -10, -20, -30, -40 };
            
            double[] posibility = { 0, 0, 0 };

            var distance = GetDistance(aircraftPoint.x, targetPoint.TargetPointLongitude,aircraftPoint.y,targetPoint.TargetPointLatitude);

            if (distance >= 0 && distance < 2)
            {
                posibility = MuShiModel.GetPosibilityClear(distance, a, b, c, d, r, e, visibility, wavehigh);
                while (posibility[0] <= 0 || posibility[1] <= 0 || posibility[2] <= 0)
                {
                    distance += 0.1;
                    posibility = MuShiModel.GetPosibilityClear(distance, a, b, c, d, r, e, visibility, wavehigh);
                }
            }
            if (distance >= 2)
            {
                posibility = MuShiModel.GetPosibilityFar(distance, a, b, c, d);
                while (posibility[0] <= 0 || posibility[1] <= 0 || posibility[2] <= 0)
                {
                    distance += 0.5;
                    posibility = MuShiModel.GetPosibilityFar(distance, a, b, c, d);
                }
                if (distance > 3.5)
                {
                    posibility[0] = 0;
                }
                if (distance > 5)
                {
                    posibility[1] = 0;
                }
                if (distance > 10)
                {
                    posibility[2] = 0;
                }
            }

            return posibility[GetType(targetPoint.TargetType.Type)];
        }
        
        public double GetMushiLandProbability(double3 aircraftPoint, double visibility, TargetPoint targetPoint)
        {
            r[0] = new double[] { 4.375, 1.18, 1.01, 1 };
            r[1] = new double[] { 2.85, 1.51, 1.18, 1 };
            r[2] = new double[] { 2.75, 2.1, 1.2, 1 };
            e[0] = new double[] { 0, -30, -50, -65, -80, -90 };
            e[1] = new double[] { 0, -15, -30, -50, -70, -80 };
            e[2] = new double[] { 0, 0, -10, -20, -30, -40 };
            
            double[] posibility = { 0, 0, 0 };

            var distance = GetDistance(aircraftPoint.x, targetPoint.TargetPointLongitude,aircraftPoint.y,targetPoint.TargetPointLatitude);

            if (distance >= 0 && distance < 2)
            {
                posibility = MuShiModel.GetPosibilityClear(distance, a, b, c, d, r, e, visibility, 0);
                while (posibility[0] <= 0 || posibility[1] <= 0 || posibility[2] <= 0)
                {
                    distance += 0.1;
                    posibility = MuShiModel.GetPosibilityClear(distance, a, b, c, d, r, e, visibility, 0);
                }
            }
            if (distance >= 2)
            {
                posibility = MuShiModel.GetPosibilityFar(distance, a, b, c, d);
                while (posibility[0] <= 0 || posibility[1] <= 0 || posibility[2] <= 0)
                {
                    distance += 0.5;
                    posibility = MuShiModel.GetPosibilityFar(distance, a, b, c, d);
                }
                if (distance > 3.5)
                {
                    posibility[0] = 0;
                }
                if (distance > 5)
                {
                    posibility[1] = 0;
                }
                if (distance > 10)
                {
                    posibility[2] = 0;
                }
            }

            return posibility[GetType(targetPoint.TargetType.Type)];
        }

        // Pd0 = 0.5 / Pf0 = Math.Pow(10,-6) / Pf = Math.Pow(10,-6) / R0 = 23645 / sigma0 = 5000 / sigma = Editor雷达截面面积 // R 单位m 
        public double GetRadarPossibility(double Pd0, double Pf0, double Pf, double R0, double3 aircraftPoint , TargetPoint targetPoint, double sigma, double sigma0, double visibility)
        {
            var distance = GetDistance(aircraftPoint.x, targetPoint.TargetPointLongitude, aircraftPoint.y, targetPoint.TargetPointLatitude);
            double R = 1000 * distance;
            //修正系数
            double RadarC = 1;

            //能见度单位:km
            if (visibility >= 1)
            {
                RadarC = 1;
            }
            else if (visibility >= 0.3 && visibility < 1)
            {
                RadarC = 0.8;
            }
            else if (visibility >= 0.09 && visibility < 0.3)
            {
                RadarC = 0.6;
            }
            else
            {
                RadarC = 0.2;
            }

            double Possibility0 = 0;
            double SNR = 1;

            //SNR为信噪比；sigma为实际目标雷达截面积，单位：平方米；sigma0为理论目标雷达截面积，单位：平方米
            SNR = sigma * Math.Pow(R0, 4) / (sigma0 * Math.Pow(R, 4));

            Possibility0 = Math.Pow(Math.E, Math.Log(Pd0, Math.E) * Math.Log(Pf, Math.E) / (SNR * Math.Log(Pf0, Math.E) + Math.Log(Pd0, Math.E)));

            double Possibility = RadarC * Possibility0;

            return Possibility;
        }

        // Lt = 探测目标亮度 Editor / At = 探测目标面积 Editor / τa = 1 / Lb = 3 / A0 = 1 / D0 = 0.075 / Dstar = 3 / τo = 0.8 / Ad = 0.0073728 / Δf = 0.125 / δ = 0.5 / Pf0 = Math.Pow(10, -6);
        public double GetInfraredDetectionProbability(double Lt, double At, double τa, double Lb, double A0, double3 aircraftPoint, TargetPoint targetPoint, double D0, double Dstar, double τo, double Ad, double Δf, double δ, double Pf0)
        {
            var distance = GetDistance(aircraftPoint.x, targetPoint.TargetPointLongitude, aircraftPoint.y, targetPoint.TargetPointLatitude);
            double R = 1000 * distance;
            //确定由Pfa确定的最小信噪比
            double minSNR = Math.Sqrt(-2 * Math.Log(Pf0, Math.E));

            // 计算信噪比
            double a = (Lt * At * τa - Lb * A0) * Math.PI * Math.Pow(D0, 2) * Dstar * δ * τo;
            double b = 4 * Math.Sqrt(Ad * Δf) * Math.Pow(R, 2);
            double SNR = a / b;

            // 使用 Math.NET Numerics 库计算高斯积分
            double Pd = GaussianIntegral(SNR - minSNR);

            return Pd;
        }

        public static double GaussianIntegral(double x)
        {
            double integral = 1;

            if (x < 0)
            {
                integral = 0.5 - Math.Pow(Math.Sqrt(2 * Math.PI), -1) * Integrate.OnClosedInterval(t => Math.Exp(-Math.Pow(t, 2) / 2), 0, -x, 1e-2);
            }
            else
            {
                integral = 0.5 + Math.Pow(Math.Sqrt(2 * Math.PI), -1) * Integrate.OnClosedInterval(t => Math.Exp(-Math.Pow(t, 2) / 2), 0, x, 1e-2);
            }

            return integral;
        }

        public int GetType(string type)
        {
            switch (type)
            {
                case "遇险人员":
                case "落水人员":
                    return 0;
                case "车辆":
                case "救生筏":
                case "小于5700kg航空器":
                case "小于20m船舶":
                    return 1;
                case "大于20m船舶":
                case "大于5700kg航空器":   
                    return 2;
            }
            return -1;
        }
        

        public static double GetDistance(double lon1, double lon2, double lat1, double lat2)
        {
            double R = 6371; // 地球的半径（公里）
            double dLat = (lat2 - lat1) * Math.PI / 180.0;
            double dLon = (lon2 - lon1) * Math.PI / 180.0;
            double a = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) +
                       Math.Cos(lat1 * Math.PI / 180.0) * Math.Cos(lat2 * Math.PI / 180.0) *
                       Math.Sin(dLon / 2) * Math.Sin(dLon / 2);
            double c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
            double distance = R * c;

            return distance;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="aircraftPoint">飞机坐标</param>
        /// <param name="targetPoint">目标点坐标</param>
        /// <param name="course">航向（度数）</param>
        /// <param name="currentDate">（当前日期）</param>
        /// <param name="windspeed">风速（km/h）</param>
        /// <param name="wavehigh">浪高（m）</param>
        /// <param name="targetpye">遇险目标类型</param>
        /// <param name="typee">救援场景类型</param>
        /// <returns></returns>
        public double getProbability(double3 aircraftPoint, double3 targetPoint, double pb, double course, double windspeed, double wavehigh, string targetpye, string type)
        {
            try
            {
                double Cw = 1;

                if (windspeed > 46 && wavehigh > 1.5)
                {
                    if (targetpye == "落水人员")
                    {
                        Cw = 0.25;
                    }
                    else
                    {
                        Cw = 0.6;
                    }
                }
                else if (28 < windspeed && windspeed < 46 && 1 < wavehigh && wavehigh < 1.5)
                {
                    if (targetpye == "船")
                    {
                        Cw = 0.5;
                    }
                    else
                    {
                        Cw = 0.9;
                    }
                }
                else if (0 < windspeed && windspeed < 28 && 0 < wavehigh && wavehigh < 1)
                {
                    Cw = 1;
                }
                double Cv = 1;

                System.Random ran = new System.Random();
                double x = getx(aircraftPoint, targetPoint, course);

                double V = 15;//getVisibility(getCityName(targetPoint.x, targetPoint.y), DateTime.Now.ToString("yyyy-MM-dd HH"));


                if (V < 6)
                {
                    Cv = 0.4;
                }
                else if (6 < V && V < 9)
                {
                    Cv = 0.6;
                }
                else if (9 < V && V < 19)
                {
                    Cv = 0.8;
                }
                else if (19 < V && V < 28)
                {
                    Cv = 0.9;
                }
                else if (37 < V)
                {
                    Cv = 1;
                }


                double pt = this.nextDouble(ran, 0.800, 0.999, 8);
                //double pb = targetPoint;//getPb(targetPoint.x, targetPoint.y);
                double D = 0;
                double C = Math.Round((double)((pt - pb) / pt), 8);
                if (type == "陆地")
                {
                    D = (V / 3.912) * Math.Log10(C / 0.02);
                }
                else if (type == "海上")
                {
                    if (targetpye == "船")
                    { D = 31 * Cv * Cw; }
                    else if (targetpye == "落水人员")
                    { D = 0.2 * Cv * Cw; }
                }
                double px = 1 - Math.Exp((-1 * D * D) / (4 * Math.PI * x * x));
                return px;
            }
            catch (Exception ex)
            {
                Debug.Print("error!!!!!getProbability:" + ex.ToString());
                return -1;
            }
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="aircraftPoint">飞机坐标</param>
        /// <param name="targetPoint">目标点坐标</param>
        /// <param name="course">航向（度数）</param>
        /// <param name="currentDate">（当前日期）</param>
        /// <returns></returns>
        public double getProbability(double3 aircraftPoint, double3 targetPoint, double course)
        {
            try
            {
                x = getx(aircraftPoint, targetPoint, course);

                double V = getVisibility(getCityName(targetPoint.x, targetPoint.y), DateTime.Now.ToString("yyyy-MM-dd HH"));

                System.Random ran = new System.Random();
                double pt = this.nextDouble(ran, 0.800, 0.999, 8);

                double pb = getPb(targetPoint.x, targetPoint.y);

                double C = Math.Round((double)((pt - pb) / pt), 8);
                double D = (V / 3.912) * Math.Log10(C / 0.02);

                double px = 1 - Math.Exp((-1 * D * D) / (4 * Math.PI * x * x));
                return px;
            }
            catch (Exception ex)
            {
                Console.WriteLine("error!!!!!getProbability:" + ex.ToString());
                return -1;
            }
        }

        public double getX(double3 targetPoint)
        {
            try
            {
                double result = -1;

                double V = getVisibility(getCityName(targetPoint.x, targetPoint.y), DateTime.Now.ToString("yyyy-MM-dd HH"));

                double pt = 0.85;

                double pb = getPb(targetPoint.x, targetPoint.y);

                double C = Math.Round((double)((pt - pb) / pt), 3);
                double D = (V / 3.912) * Math.Log10(C / 0.02);

                double x = Math.Sqrt(-1 * D * D / Math.Log(1 - 0.85) / 4 / Math.PI);
                return x;
            }
            catch (Exception ex)
            {
                Debug.Print("error!!!!!getX:" + ex.ToString());
                return -1;
            }
        }

        private double nextDouble(System.Random ran, double minValue, double maxValue, int decimalPlace)
        {
            double randNum = ran.NextDouble() * (maxValue - minValue) + minValue;
            return Convert.ToDouble(randNum.ToString("f" + decimalPlace));
        }

        public double getVisibility(string city_name, string weather_date)
        {
            try
            {
                //调在线气象接口，封装到服务中了
                HttpCmd cmd = new HttpCmd
                {
                    m_RequestType = HttpRequestType.GET,
                    m_Addr = $"{baseUrl}rescue-platform-service/api/v1/tbMeteorology/getInfosByNowApi",
                    m_Args = new List<string> { "city_name", "weather_date" }
                };
                string response = m_HttpHelper.Request(cmd, new List<string> { city_name, weather_date });
                R data = JsonConvert.DeserializeObject<R>(response);
                if (data != null && data.code == 200)
                {
                    NAdtListDto rt = JsonConvert.DeserializeObject<NAdtListDto>(data.data.ToString());
                    double fl = Convert.ToDouble(rt.wtVisibility);
                    if (fl.Equals(0))
                    {
                        fl = 15;
                    }
                    return fl;
                }
                return 15;
            }
            catch (Exception ex)
            {
                Debug.Print("error!!!!!getVisibility");
                Debug.Print(ex.ToString());
                return 15;
            }
        }

        public double getVisibilityByDb(double centerLon, double centerLat, string visibility_date)
        {
            try
            {
                HttpCmd cmd = new HttpCmd
                {
                    m_RequestType = HttpRequestType.GET,
                    m_Addr = "http://10.130.100.5:7785//rescue-platform-service/api/v1/attribute/getVisibility", //10.130.100.5 127.0.0.1
                    m_Args = new List<string> { "centerLon", "centerLat", "visibility_date" }
                };
                string response = m_HttpHelper.Request(cmd, new List<string> { centerLon.ToString(), centerLat.ToString(), visibility_date });
                //Console.WriteLine(" centerLon:" + centerLon + " " + "centerLat:" + centerLat + " " + "visibility_date:" + visibility_date);
                //Console.WriteLine("response:" + response);
                R data = JsonConvert.DeserializeObject<R>(response);
                if (data != null && data.code == 200)
                {
                    string vi = data.data.ToString();
                    double fl = Convert.ToDouble(vi);
                    if (fl.Equals(0))
                    {
                        fl = 15;
                        //Console.WriteLine("1:" + fl);
                    }
                    //Console.WriteLine("2:" + fl);
                    return fl;
                }
                //Console.WriteLine("3:" + 15);
                return 15;
            }
            catch (Exception ex)
            {
                //Console.WriteLine("4:" + 15);
                Debug.Print("error!!!!!getVisibility");
                Debug.Print(ex.ToString());
                return 15;
            }
        }

        //rescue-platform-service/api/v1/dem/getCityName

        // {
        //     "msg": "success",
        //     "code": 200,
        //     "cityName": "承德市"
        // }
        public string getCityName(double lon, double lat)
        {
            string cityName = "北京";
            // try
            // {
            //     string sql = string.Format("select a.f_xzqdm cityCode, a.f_xzqmc cityName from xzq_ds a where ST_Contains(shape,st_geomfromtext('POINT({0} {1})',4326)) limit 1", lon, lat);
            //     DataTable dt = this._db.DoQueryEx(sql);
            //     if (dt != null && dt.Rows.Count > 0)
            //     {
            //         cityName = dt.Rows[0][1].ToString().Replace("市", "").Replace("自治区", "");
            //     }
            // }
            // catch (Exception ex)
            // {
            //     Debug.Print("error!!!!!getCityName");
            //     Debug.Print(ex.ToString());
            // }

            try
            {
                //调在线气象接口，封装到服务中了
                HttpCmd cmd = new HttpCmd
                {
                    m_RequestType = HttpRequestType.GET,
                    m_Addr = $"{baseUrl}rescue-platform-service/api/v1/dem/getCityName",
                    m_Args = new List<string> { "lon", "lat" }
                };
                string response = m_HttpHelper.Request(cmd, new List<string> { lon.ToString(), lat.ToString() });
                JObject jObject = JObject.Parse(response);
                cityName = jObject["cityName"].ToString().Replace("市", "").Replace("自治区", "");

            }
            catch (Exception ex)
            {
                Debug.Print("error!!!!!getVisibility");
                Debug.Print(ex.ToString());
            }
            return cityName;
        }

        //rescue-platform-service/api/v1/dem/getCategory
        // {
        //     "msg": "success",
        //     "code": 200,
        //     "category": "7 草丛"
        // }
        public double getPb(double lon, double lat)
        {
            double pb = 0.13;
            try
            {
                // string sql_zb = string.Format("select f_category from ly_vegetation where  ST_Contains(st_geomfromtext(st_astext(shape) ),st_geomfromtext('POINT({0} {1})'))", lon, lat);
                // DataTable dt_zb = _db.DoQueryEx(sql_zb);

                //调在线气象接口，封装到服务中了
                HttpCmd cmd = new HttpCmd
                {
                    m_RequestType = HttpRequestType.GET,
                    m_Addr = $"{baseUrl}rescue-platform-service/api/v1/dem/getCategory",
                    m_Args = new List<string> { "lon", "lat" }
                };
                string response = m_HttpHelper.Request(cmd, new List<string> { lon.ToString(), lat.ToString() });

                JObject jObject = JObject.Parse(response);
                string fcategory = jObject["category"].ToString();
                switch (fcategory)
                {
                    case "1 针叶林":
                        pb = 0.13;
                        break;
                    case "10 高山植被":
                        pb = 0.13;
                        break;
                    case "11 栽培植被":
                        pb = 0.13;
                        break;
                    case "2 针阔叶混交林":
                        pb = 0.13;
                        break;
                    case "3 阔叶林":
                        pb = 0.13;
                        break;
                    case "4 灌丛":
                        pb = 0.13;
                        break;
                    case "5 荒漠":
                        pb = 0.13;
                        break;
                    case "6 草原":
                        pb = 0.13;
                        break;
                    case "7 草丛":
                        pb = 0.13;
                        break;
                    case "8 草甸":
                        pb = 0.13;
                        break;
                    case "9 沼泽":
                        pb = 0.13;
                        break;
                    default:
                        pb = 0.13;
                        break;
                }
            }
            catch (Exception ex)
            {
                Debug.Print("error!!!!!getPb");
                Debug.Print(ex.ToString());
            }
            return pb;
        }

        //rescue-platform-service/api/v1/dem/getDistance

        // {
        //     "msg": "success",
        //     "code": 200,
        //     "data": {
        //         "differ": 328.31408076,
        //         "azimuth": 2.400469679065642
        //     }
        // }
        private double getx(double3 aircraftPoint, double3 targetPoint, double course)
        {
            double result = -1;
            // string sql = string.Format("select ST_Distance( ST_SetSRID(ST_MakePoint({0},{1}),4326)::geography, ST_SetSRID(ST_MakePoint({2},{3}),4326)::geography) as differ,ST_azimuth(st_geomfromtext('POINT({0} {1})'),st_geomfromtext('POINT({2} {3})')) as azimuth", aircraftPoint.x, aircraftPoint.y, targetPoint.x, targetPoint.y);
            // DataTable dt = _db.DoQueryEx(sql);
            HttpCmd cmd = new HttpCmd
            {
                m_RequestType = HttpRequestType.GET,
                m_Addr = $"{baseUrl}rescue-platform-service/api/v1/dem/getDistance",
                m_Args = new List<string> { "lon_air", "lat_air", "lon_target", "lat_target" }
            };
            string response = m_HttpHelper.Request(cmd, new List<string> { aircraftPoint.x.ToString(), aircraftPoint.y.ToString(), targetPoint.x.ToString(), targetPoint.y.ToString() });

            JObject jObject = JObject.Parse(response);


            double differ = double.Parse(jObject["data"]["differ"].ToString());
            double azimuth = double.Parse(jObject["data"]["azimuth"].ToString());

            // differ = Convert.ToDouble(dt.Rows[0][0].ToString());
            azimuth = 180 / Math.PI * azimuth;//Convert.ToDouble(dt.Rows[0][1].ToString());
            double angle = course - azimuth;
            if (angle < 0)
            {
                angle = -1 * angle;
            }
            double x = differ * Math.Sin(angle * Math.PI / 180) / 1000;
            return x;
        }
    }
}