计算任意时刻格林尼治视恒星时角

近期学习卫星轨道方面的一些知识，遇到计算任意时刻格林尼治视恒星时角的问题，在网上搜了好久也没有一个完整的解决方案，后来通过，网上的一些零碎的信息，终于完成了计算格林尼治视恒星时角的程序，先整理如下。

计算格林尼治视恒星时角，首先需要计算当前时间的儒略日，计算方法如下

设Y为给定的年份，M为给定的月份，D为给定的日期。运算符INT为取所给数据的整数部分，若M大于2则Y与M不变，否则M加上12，Y减去1，换句话说，如果M为1月或二月，则被看做是前一年的13月和14月。

对于格里高利历（1582年10月15日以后），有

A=INT（Y/100）

B=2-A+INT（A/4）。

对于儒略历（即1582年10月15日以前），B=0。

所求儒略日即为

JD=INT（365.25（Y+4716））+INT（30.6001*（M+1））+D+B-1524.5

实现代码如下：

            const long IGREG = (15 + 31 * (10 + 12 * 1582));
            int jul;
            int ja, jy = iyyy, jm;

            if (jy == 0)
            {
                MessageBox.Show("ERROR in subroutine julday: there is no year zero!");
                return 0;

            }

            if (jy < 0)
                ++jy;
            if (mm > 2)
            {
                jm = mm + 1;
            }
            else
            {
                --jy;
                jm = mm + 13;
            }
            jul = (int)(Math.Floor(365.25 * (jy+4716)) + Math.Floor(30.6001 * jm) + id - 1524.5);
            if (id + 31 * (mm + 12 * iyyy) >= IGREG)//判断是否为格里高利历日IGREG
            {
                ja = (int)(0.01 * jy);
                jul += 2 - ja + (int)(0.25 * ja);//加百年闰
            }
            return jul;

然后计算任意时刻的格林尼治视恒星时角，需要首先计算格林尼治平恒星时Θm，计算公式如下

Θm=100.46061837 + 36000.770053608T + 0.000387933T*T  - T*T*T / 38710000

其中T =(JD - 2451545)/36525 ，JD为相应的儒略日，实现代码如下

double calGMST(double jd，double st)
{
            double D = jd - 2451545.0;//Compute the number of days since Jully2000
            
            double T = D / 36525.0;
            double GMST = 6.697374558 + 2400.051336*T + 0.000025862*T*T+st / 3600 *1.00273790935
            GMST = GMST % 24 * 15;
            return GMST;

}

然后计算格林尼治视恒星时角Θ，计算公式如下

Θ = Θm + ΔΨ cos(εm + Δε) 

其中εm为黄道倾角，计算公式如下：

εm=23°26′21.″448 - 46.″8150T - 0.″00059T*T  + 0.″001813T*T*T

倾角章动和经度设计一下三个三角参数

L = 280.4665 + 36000.7698T
L′ = 218.3165 + 481267.8813T
Ω = 125.04452 - 1934.136261T

计算倾角章动用如下公式

ΔΨ= -17.20 sinΩ - 1.32 sin2 L - 0.23 sin2 L′ + 0.21sin2 Ω

Δε= 9.20 cos Ω + 0.57cos2 L + 0.10 cos 2 L′ - 0.09 cos 2 Ω

这些参数均以弧秒为单位

实现代码如下

double calGAST(double day,double yr,double mo,double st)
            double fjd = fulljd(day, yr, mo);
            double TJD = fjd - 2451545.0;//2451545为2000 1月1.5的儒略日
            double T0 = TJD / 36525.0;
            double THETAm = calGMST(fjd,st);
            double EPSILONm = 23.4392911111111111 - 0.0130041666667* T0 - 1.638e-07 * T0 * T0 + 5.0361e-07 * T0 * T0 * T0;

            double L = 280.4665 + 36000.7698 * T0;
            double dL = 218.3165 + 481267.8813 * T0;
            double OMEGA = 125.04452 - 1934.136261 * T0;
            double dPSI = -17.20 * Math.Sin(OMEGA) - 1.32 * Math.Sin(2 * L) - 0.23 * Math.Sin(2 * dL) + 0.21 * Math.Sin(2 * OMEGA);
            double dEPSILON = 9.20 * Math.Cos(OMEGA) + 0.57 * Math.Cos(2 * L) + 0.10 * Math.Cos(2 * dL) - 0.09 * Math.Cos(2 * OMEGA);

            //Convert the units from arc-seconds to degrees
            dPSI /= 3600;
            dEPSILON /= 3600;

            double GAST = (THETAm + dPSI * Math.Cos(EPSILONm + dEPSILON)) % 360;
            return GAST;
}