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

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

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

            double fjd = fulljd(0, yr, 1);//计算当年1月0日的儒略日
            double TJD = fjd - 2415020.0;
            double T0 = TJD / 36525.0;            
            double r = 6.6460656 + 2400.051262 * T0 + 0.00002581 * T0 * T0;

            double u = r - 24.0 * (yr - 1900);
            double bigbee = 24.0 - u;
            double fjd1 = fulljd(day, yr, mo);//计算当日的儒略日
            double deltfjd = fjd1 - fjd;//当天在改年的第几天
            double GAST0 = 0.0657098 * deltfjd - bigbee + st / 3600 * 1.00273790935;//st为具体的时间以秒为单位（st=h+3600+min*60+sec）
            if (GAST0 >= 24)
                GAST0 -= 24.0;
            if (GAST0 < 0)
                GAST0 += 24.0;
            GAST0 *= 15 * DTR;
            return GAST0;