SimulationServer/Models/SimulationCommon/SeaSJ.cs

using OpenCvSharp;

namespace SimulationCommon;

public class WeatherResponse
{
    public Wind Wind { get; set; }
}

public class Wind
{
    public double Speed { get; set; }
    public double Deg { get; set; }
}

public class CurrentResponse
{
    public Current Current { get; set; }
}

public class Current
{
    public double Speed { get; set; }
    public double Direction { get; set; }
}

public class NCread
{

    //各数组来源文件'Text_readNC.cs'；时间范围：2024-06-04T00:00:00 ... 2024-06-04T23:00:00；
    //地点范围:落水点为中心100公里
    public float[] longitudeArray = new float[10];//经度一维数组来源自文件'Text_readNC.cs'
    public float[] latitudeArray = new float[7];//纬度一维数组来源自文件'Text_readNC.cs'
    public float[][][] u10Array;//风的10米U（向东）分量三维数组来源自文件'Text_readNC.cs'
    public float[][][] v10Array;//风的10米V（向北）分量三维数组来源自文件'Text_readNC.cs'
    public float[][][] p140208Array;//海洋上空的自由对流速度三维数组来源自文件'Text_readNC.cs'
    public float[][][] mwdArray;//平均波向（单位：度；0度表示北方，90度表示东方）三维数组来源自文件'Text_readNC.cs'       
    public float[][][] hmaxArray;

}

public class SeaSJ
{
    // 生成符合正态分布的随机数
    public static double NormalDistributionRandom(double mu, double sigma)
    {
        Random rand = new Random();
        double u1 = 1.0 - rand.NextDouble(); // uniform(0,1] random doubles
        double u2 = 1.0 - rand.NextDouble();
        double normalRandom = Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Sin(2.0 * Math.PI * u2); // random normal(0,1)
        return mu + sigma * normalRandom; // random normal(mean,stdDev^2)
    }

    // 交付类方法
    public static List<double[]> GetDrift(NCread nCread, double[] initialPosition, double dt, double totalTime, int times, int latitudes, int longitudes, int times1, int latitudes1, int longitudes1, int times2, int latitudes2, int longitudes2, int times3, int latitudes3, int longitudes3, int days, int hour, double initlatitudes, double initlongitudes)
    {
        // run 获取轨迹的函数
        // initialPosition --> 初始位置; dt --> 时间步长; totalTime --> 总时长
        List<double[]> trajectory = CalculateDriftTrajectory(nCread, initialPosition, dt, totalTime, times, latitudes, longitudes, times1, latitudes1, longitudes1, times2, latitudes2, longitudes2, times3, latitudes3, longitudes3, days, hour, initlatitudes, initlongitudes);
        return trajectory;
    }

    public static List<double[]> CalculateDriftTrajectory(NCread nCread, double[] initialPosition, double dt, double totalTime, int times, int latitudes, int longitudes, int times1, int latitudes1, int longitudes1, int times2, int latitudes2, int longitudes2, int times3, int latitudes3, int longitudes3, int days, int hour, double initlatitudes, double initlongitudes)
    {
        int timeSteps = (int)(totalTime / dt);
        List<double[]> trajectory = new List<double[]>();
        trajectory.Add(new double[] { initialPosition[0], initialPosition[1] });

        for (int t = 1; t < timeSteps; t++)
        {
            double[] currentPos = trajectory[t - 1];

            int outputTime = 3600 * t;
            // 动态获取当前位置的风力和洋流数据
            double[] windVelocity = GetWindVelocityFromAPI(nCread, currentPos[0], currentPos[1], outputTime, times, latitudes, longitudes, times1, latitudes1, longitudes1, days, hour);
            double[] currentVelocity = GetCurrentVelocityFromAPI(nCread, currentPos[0], currentPos[1], outputTime, times2, latitudes2, longitudes2, times3, latitudes3, longitudes3, days, hour);
            //foreach (var item in windVelocity)
            //{
            //    Console.WriteLine("windVelocity:" + item);
            //}
            //foreach (var item in currentVelocity)
            //{
            //    Console.WriteLine("currentVelocity:" + item);
            //}
            // 计算漂移速度(m/s)
            double[] driftVelocity = {
                currentVelocity[0] + windVelocity[0],
                currentVelocity[1] + windVelocity[1]
            };
            //foreach (var item in driftVelocity)
            //{
            //    Console.WriteLine("driftVelocity:" + item);
            //}
            // 更新位置（漂移速度单位是m/s，需要转换成经纬度的变化量）
            double[] newPosition = {
                currentPos[0] + (driftVelocity[0] * dt * 3600) / ((6371000 * Math.Cos(currentPos[0] * Math.PI / 180)) * 180 / Math.PI), // 纬度变化量
                currentPos[1] + (driftVelocity[1] * dt * 3600) / ((6371000 * Math.Cos(currentPos[1] * Math.PI / 180)) * 180 / Math.PI)  // 经度变化量
            };
            //foreach (var item in newPosition)
            //{
            //    Console.WriteLine("newPosition:" + item);
            //}
            //Console.WriteLine("currentPos[0]:" + (driftVelocity[0] * dt * 3600) / ((6371000 * Math.Cos(currentPos[0] * Math.PI / 180)) * 180 / Math.PI));
            //Console.WriteLine("currentPos[1]:" + (driftVelocity[1] * dt * 3600) / ((6371000 * Math.Cos(currentPos[1] * Math.PI / 180)) * 180 / Math.PI));

            trajectory.Add(newPosition);
        }

        return trajectory;//（纬度，经度）
    }

    public static double[] GetWindVelocityFromAPI(NCread windNCread, double latitude, double longitude, double temptime, int times, int latitudes, int longitudes, int times1, int latitudes1, int longitudes1, int days, int hour)
    {
        float[] longitudeArray = windNCread.longitudeArray;//经度一维数组来源自文件'Text_readNC.cs'
        float[] latitudeArray = windNCread.latitudeArray;//纬度一维数组来源自文件'Text_readNC.cs'
        float[][][] u10Array = windNCread.u10Array;//风的10米U（向东）分量三维数组来源自文件'Text_readNC.cs'
        float[][][] v10Array = windNCread.v10Array;//风的10米V（向北）分量三维数组来源自文件'Text_readNC.cs'

        int longitudeNum = 0;
        int latitudeNum = 0;
        int temptimeNum = 0;
        //定义NC文件中读取不到的坐标点海洋数据信息（后续可采用插值法等）
        //经度连续化
        for (int i = 0; i < longitudes - 1; i++)
        {
            if (longitude >= longitudeArray[i] && longitude < longitudeArray[i + 1])
            {
                longitude = longitudeArray[i];
                longitudeNum = i;
            }
        }

        //纬度连续化
        for (int i = 0; i < latitudes - 1; i++)
        {
            if (latitude >= latitudeArray[i] && latitude < latitudeArray[i + 1])
            {
                latitude = latitudeArray[i];
                latitudeNum = i;
            }
        }

        //时间连续化
        for (int i = 0; i < times - 1; i++)
        {
            if (temptime >= 3600 * i && temptime < (3600 * (i + 1)))
            {
                temptimeNum = i;
            }
        }

        temptimeNum += (days - 1) * 24 + hour;

        double windX = (double)u10Array[temptimeNum][latitudeNum][longitudeNum];
        double windY = (double)v10Array[temptimeNum][latitudeNum][longitudeNum];

        return new double[] { windX, windY };
    }

    public static double[] GetCurrentVelocityFromAPI(NCread CurrentNCread, double latitude, double longitude, double temptime, int times2, int latitudes2, int longitudes2, int times3, int latitudes3, int longitudes3, int days, int hour)
    {
        float[] longitudeArray = CurrentNCread.longitudeArray;//经度一维数组来源自文件'Text_readNC.cs'
        float[] latitudeArray = CurrentNCread.latitudeArray;//纬度一维数组来源自文件'Text_readNC.cs'
        float[][][] u10Array = CurrentNCread.u10Array;//风的10米U（向东）分量三维数组来源自文件'Text_readNC.cs'
        float[][][] v10Array = CurrentNCread.v10Array;//风的10米V（向北）分量三维数组来源自文件'Text_readNC.cs'
        float[][][] p140208Array = CurrentNCread.p140208Array;//海洋上空的自由对流速度三维数组来源自文件'Text_readNC.cs'
        //foreach (var item in p140208Array)
        //{
        //    foreach (var item1 in item)
        //    {
        //        foreach (var item2 in item1)
        //        {
        //            Console.WriteLine("p140208Array:" + item2);
        //        }
        //    }
        //}
        float[][][] mwdArray = CurrentNCread.mwdArray;//平均波向（单位：度；0度表示北方，90度表示东方）三维数组来源自文件'Text_readNC.cs'

        int longitudeNum = 0;
        int latitudeNum = 0;
        int temptimeNum = 0;
        //定义NC文件中读取不到的坐标点海洋数据信息（后续可采用插值法等）
        //经度连续化
        for (int i = 0; i < longitudes2 - 1; i++)
        {
            if (longitude >= longitudeArray[i] && longitude < longitudeArray[i + 1])
            {
                longitude = longitudeArray[i];
                longitudeNum = i;
            }
        }

        //纬度连续化
        for (int i = 0; i < latitudes2 - 1; i++)
        {
            if (latitude >= latitudeArray[i] && latitude < latitudeArray[i + 1])
            {
                latitude = latitudeArray[i];
                latitudeNum = i;
            }
        }

        //时间连续化
        for (int i = 0; i < times2 - 1; i++)
        {
            if (temptime >= 3600 * i && temptime < 3600 * (i + 1))
            {
                temptimeNum = i;
            }
        }
        temptimeNum += (days - 1) * 24 + hour;
        //foreach (var item in p140208Array)
        //{
        //    Console.WriteLine("p140208Array:" + (double)p140208Array[0][1][2]);
        //}
        //Console.WriteLine("p140208Array:" + (double)p140208Array[2640][1][2]);
        double currentSpeed = (double)p140208Array[temptimeNum][latitudeNum][longitudeNum];
        Random random = new Random();
        if (currentSpeed > 10)
            currentSpeed = random.NextDouble();
        //Console.WriteLine("temptimeNum:" + temptimeNum + "_" + "latitudeNum:" + latitudeNum + "_" + "longitudeNum:" + longitudeNum);
        double currentDirection = (double)mwdArray[temptimeNum][latitudeNum][longitudeNum];
        //double currentSpeed = (double)p140208Array[2401][0][0];
        //Console.WriteLine("currentSpeed:" + currentSpeed);
        //Console.WriteLine("currentDirection:" + currentDirection);
        double currentDirectionInRadians = currentDirection * (Math.PI / 180);
        double currentX = currentSpeed * Math.Cos(currentDirectionInRadians);
        double currentY = currentSpeed * Math.Sin(currentDirectionInRadians);

        return new double[] { currentX, currentY };
    }

    //海浪高度获取
    public static double GetWaveHeightFromAPI(NCread WaveNCread, double latitude, double longitude, double temptime, int times4, int latitudes4, int longitudes4, int days, int hour)//temptime:仿真时间（秒）
    {
        float[] longitudeArray = WaveNCread.longitudeArray;//经度一维数组来源自文件'Text_readNC.cs'
        float[] latitudeArray = WaveNCread.latitudeArray;//纬度一维数组来源自文件'Text_readNC.cs'
        float[][][] hmaxArray = WaveNCread.hmaxArray;//风的10米U（向东）分量三维数组来源自文件'Text_readNC.cs'

        int longitudeNum = 0;
        int latitudeNum = 0;
        int temptimeNum = 0;

        //定义NC文件中读取不到的坐标点海洋数据信息（后续可采用插值法等）
        //经度连续化
        for (int i = 0; i < longitudes4 - 1; i++)
        {
            if (longitude >= longitudeArray[i] && longitude < longitudeArray[i + 1])
            {
                longitude = longitudeArray[i];
                longitudeNum = i;
            }
        }

        //纬度连续化
        for (int i = 0; i < latitudes4 - 1; i++)
        {
            if (latitude >= latitudeArray[i] && latitude < latitudeArray[i + 1])
            {
                latitude = latitudeArray[i];
                latitudeNum = i;
            }
        }

        //时间连续化
        for (int i = 0; i < times4 - 1; i++)
        {
            if (temptime >= 3600 * i && temptime < (3600 * (i + 1)))
            {
                temptimeNum = i;
            }
        }
        temptimeNum += (days - 1) * 24 + hour;
        double WaveHeight = (double)hmaxArray[temptimeNum][latitudeNum][longitudeNum];

        return WaveHeight;
    }

    // 交付部分方法代码
    public static List<double[]> getminEnclosingRect(List<double[]> latLonList)
    {
        // 转换经纬度为墨卡托投影坐标，并添加到点集合
        List<Point2f> pointList = new List<Point2f>();
        foreach (var latLon in latLonList)
        {
            double x = Rectangular_Area_Search_Function.MokatuoLat(latLon[1]);
            double y = Rectangular_Area_Search_Function.MokatuoLon(latLon[0]);
            pointList.Add(new Point2f((float)y, (float)x));
        }

        // 获取凸包
        Point2f[] convexHull = Rectangular_Area_Search_Function.GetConvexHull(pointList);

        //获取凸包各点经纬度坐标
        List<double[]> hullPoint = new List<double[]>();
        for (int num = 0; num < convexHull.Length - 1; num++)
        {
            hullPoint.Add(new double[] { convexHull[num].Y, convexHull[num].X });
        }

        List<double[]> hullPointLatLon = new List<double[]>();
        foreach (var point in hullPoint)
        {
            double pointLat = Rectangular_Area_Search_Function.RMokatuoLat(point[1]);
            double pointLon = Rectangular_Area_Search_Function.RMokatuoLon(point[0]);
            hullPointLatLon.Add(new double[] { pointLat, pointLon });
        }

        // 计算最小包围矩形
        List<double[]> minEnclosingRect = Rectangular_Area_Search_Function.MinEnclosingRectangle(convexHull);
        //foreach (var point in minEnclosingRect)
        //{
        //    foreach (var item in point)
        //    {
        //        Console.WriteLine("item:" + item);
        //    }
        //}

        // 最小包围矩形顶点经纬度坐标
        List<double[]> startPoint = new List<double[]>();
        foreach (var minEnclosingRectPoint in minEnclosingRect)
        {
            double lat = Rectangular_Area_Search_Function.RMokatuoLat(minEnclosingRectPoint[0]);
            double lon = Rectangular_Area_Search_Function.RMokatuoLon(minEnclosingRectPoint[1]);
            startPoint.Add(new double[] { lat, lon });
        }
        //foreach (var point in startPoint)
        //{
        //    foreach (var item in point)
        //    {
        //        Console.WriteLine("startPoint:" + item);
        //    }
        //}
        return startPoint;
    }
}