近期学习卫星轨道方面的一些知识,遇到计算任意时刻格林尼治视恒星时角的问题,在网上搜了好久也没有一个完整的解决方案,后来通过,网上的一些零碎的信息,终于完成了计算格林尼治视恒星时角的程序,先整理如下。
计算格林尼治视恒星时角,首先需要计算当前时间的儒略日,计算方法如下
设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;
计算任意时刻格林尼治视恒星时角,布布扣,bubuko.com
原文地址:http://blog.csdn.net/zhuimengshizhe87/article/details/26702997