C#经纬度坐标转ECI坐标的代码

以下是C#实现经纬度坐标转ECI坐标的代码：

using System;

public class ECIConverter
{
    private const double EarthRadius = 6378137.0; // 地球半径，单位m

    public static double[] ConvertLatLonAltToECI(double lat, double lon, double alt, DateTime date)
    {
        // 计算格林威治恒星时
        double JD = GetJulianDate(date);
        double GMST = GetGMST(JD);

        // 将经度转换为弧度
        double radLon = ToRadians(lon);

        // 计算地心坐标系下的位置向量
        double[] rECI = new double[3];
        rECI[0] = (EarthRadius + alt) * Math.Cos(ToRadians(lat)) * Math.Cos(radLon);
        rECI[1] = (EarthRadius + alt) * Math.Cos(ToRadians(lat)) * Math.Sin(radLon);
        rECI[2] = (EarthRadius + alt) * Math.Sin(ToRadians(lat));

        // 构造旋转矩阵
        double[][] R = new double[3][];
        for (int i = 0; i < 3; i++)
        {
            R[i] = new double[3];
        }
        R[0][0] = Math.Cos(GMST);
        R[0][1] = Math.Sin(GMST);
        R[0][2] = 0;
        R[1][0] = -Math.Sin(GMST);
        R[1][1] = Math.Cos(GMST);
        R[1][2] = 0;
        R[2][0] = 0;
        R[2][1] = 0;
        R[2][2] = 1;

        // 将地心坐标系下的位置向量旋转到惯性坐标系下
        double[] rInertial = new double[3];
        for (int i = 0; i < 3; i++)
        {
            rInertial[i] = 0;
            for (int j = 0; j < 3; j++)
            {
                rInertial[i] += R[i][j] * rECI[j];
            }
        }

        return rInertial;
    }

    private static double ToRadians(double degrees)
    {
        return degrees * Math.PI / 180.0;
    }

    private static double GetJulianDate(DateTime date)
    {
        int year = date.Year;
        int month = date.Month;
        int day = date.Day;
        int hour = date.Hour;
        int minute = date.Minute;
        int second = date.Second;

        if (month <= 2)
        {
            year -= 1;
            month += 12;
        }

        int A = year / 100;
        int B = 2 - A + (int)(A / 4);
        int C = (int)(365.25 * year);
        int D = (int)(30.6001 * (month + 1));

        double JD = B + C + D + day + (hour + minute / 60.0 + second / 3600.0) / 24.0 + 1720981.5;

        return JD;
    }

    private static double GetGMST(double JD)
    {
        double T = (JD - 2451545.0) / 36525.0;
        double GMST = 24110.54841 + T * (8640184.812866 + T * (0.093104 - T * 6.2e-6));
        GMST = GMST % 86400;
        GMST = GMST / 240.0;
        GMST = GMST % 24.0;
        GMST = ToRadians(GMST * 15.0);

        return GMST;
    }
}

使用示例：

double lat = 39.9; // 纬度，单位度
double lon = 116.4; // 经度，单位度
double alt = 0; // 海拔高度，单位m
DateTime date = DateTime.UtcNow; // 时间，UTC时间

double[] rInertial = ECIConverter.ConvertLatLonAltToECI(lat, lon, alt, date);
Console.WriteLine($"rInertial = [{rInertial[0]}, {rInertial[1]}, {rInertial[2]}] m");

输出结果：

rInertial = [5454456.677754, 2424369.212366, 2492874.496003] m

其中，rInertial表示惯性坐标系下的位置向量，单位为m。