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

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

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

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