一、文本格式
通常给出的地球指向参数文本格式为:
1962 01 01 37665 -0.012700 0.213000 0.0326338 0.0017230 0.064261 0.006067 0.000000 0.000000 2
1962 01 01:公历日期
37665 :儒略日
-0.0127 :x极移,单位:角秒
0.213 :y极移,单位:角秒
0.0326338:世界时UT1与协调世界时UTC的差
0.0017230:LOD 地球日长变化 单位:s
0.064261 :x向的章动值 单位:角秒
0.006067 :y向的章动值 单位:角秒
0.000000 :x向极移的变化量 单位:角秒
0.000000 :y向极移的变化量 单位:角秒
2 :国际原子时TAI与协调世界时UTC之间相差整数秒 单位:s
角秒即为1°/3600,又已知pi=180°,求得1角秒约等于4.848136811e-6
二、应用
格林尼治地方平时M称为世界时,记为UT。天文经度处的地方平时与世界时的关系为
时间的单位是小时(h),经度的单位是度(°)。通过测站直接观测天体获得的世界时记为UT0,对应于瞬时极子午圈。修正极移引起子午圈变化后得到的世界时记为UT1,对应于平均极子午圈。
根据瞬时真天球坐标系与瞬时地球坐标系的转换关系,两者的差别仅是X轴的指向差一个格林尼治真恒星时。真恒星时的计算由平恒星时和赤经岁差引起。平恒星时的计算与UT1有关
其中是该天UT1时刻的简略儒略日的整数部分与2000年1月1日12:00的简略儒略日之差除以100年即36525天,为该天的简略儒略日与2000年1月1日12:00的简略儒略日之差除以100年即36525天
则格林尼治地方平时的求解即写为
代码如下:
%--------------------------------------------------------------------------
%
% gmst: Greenwich Mean Sidereal Time
%
% Input:
% Mjd_UT1 Modified Julian Date UT1
%
% Output:
% gmstime GMST in [rad]
%
% Last modified: 2015/08/12 M. Mahooti
%
%--------------------------------------------------------------------------
function gmstime = gmst(Mjd_UT1)
Secs = 86400; % Seconds per day
MJD_J2000 = 51544.5;
Mjd_0 = floor(Mjd_UT1);
UT1 = Secs*(Mjd_UT1-Mjd_0); % [s]
T_0 = (Mjd_0 -MJD_J2000)/36525;
T = (Mjd_UT1-MJD_J2000)/36525;
gmst = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1...
+ (0.093104-6.2e-6*T)*T*T; % [s]
gmstime = 2*pi*Frac(gmst/Secs); % [rad], 0..2pi
地球自转轴的方向还受到被称为章动的小周期扰动的影响。它们是由于在进差处理中平均的月球和太阳扭矩的月和年变化。对章动的主要贡献来自于月球轨道相对于地球赤道的不同方向,这一点由月球的升交点黄经表示,它引起了周期性的变化。在本文给出IAU1980章动模型,其表达式如下式所示
共106项,其常数计入如下形式。
每个术语都描述了月球和太阳轨道的平均元素的周期函数
整数系数为。其它参数是月球的平近地点角,太阳平近地点角,月球的升交点赤经,太阳和月球的平均经度之间的差异,和月球轨道的升交点赤经。计算公式如下:
式中:T为该天的简略儒略日与2000年1月1日12:00的简略儒略日之差除以100年即36525天。其计算代码如下:
%--------------------------------------------------------------------------
%
% NutAngles: Nutation in longitude and obliquity
%
% Input:
% Mjd_TT Modified Julian Date (Terrestrial Time)
%
% Output:
% dpsi, deps
%
% Last modified: 2018/01/04 M. Mahooti
%
%--------------------------------------------------------------------------
function [dpsi, deps] = NutAngles(Mjd_TT)
global const
T = (Mjd_TT-const.MJD_J2000)/36525;
T2 = T*T;
T3 = T2*T;
rev = 360*3600; % arcsec/revolution
N_coeff = 106;
C = [
% l l' F D Om dpsi *T deps *T
0, 0, 0, 0, 1,-1719960,-1742, 920250, 89 ; % 1
0, 0, 0, 0, 2, 20620, 2, -8950, 5 ; % 2
-2, 0, 2, 0, 1, 460, 0, -240, 0 ; % 3
2, 0,-2, 0, 0, 110, 0, 0, 0 ; % 4
-2, 0, 2, 0, 2, -30, 0, 10, 0 ; % 5
1,-1, 0,-1, 0, -30, 0, 0, 0 ; % 6
0,-2, 2,-2, 1, -20, 0, 10, 0 ; % 7
2, 0,-2, 0, 1, 10, 0, 0, 0 ; % 8
0, 0, 2,-2, 2, -131870, -16, 57360, -31 ; % 9
0, 1, 0, 0, 0, 14260, -34, 540, -1 ; % 10
0, 1, 2,-2, 2, -5170, 12, 2240, -6 ; % 11
0,-1, 2,-2, 2, 2170, -5, -950, 3 ; % 12
0, 0, 2,-2, 1, 1290, 1, -700, 0 ; % 13
2, 0, 0,-2, 0, 480, 0, 10, 0 ; % 14
0, 0, 2,-2, 0, -220, 0, 0, 0 ; % 15
0, 2, 0, 0, 0, 170, -1, 0, 0 ; % 16
0, 1, 0, 0, 1, -150, 0, 90, 0 ; % 17
0, 2, 2,-2, 2, -160, 1, 70, 0 ; % 18
0,-1, 0, 0, 1, -120, 0, 60, 0 ; % 19
-2, 0, 0, 2, 1, -60, 0, 30, 0 ; % 20
0,-1, 2,-2, 1, -50, 0, 30, 0 ; % 21
2, 0, 0,-2, 1, 40, 0, -20, 0 ; % 22
0, 1, 2,-2, 1, 40, 0, -20, 0 ; % 23
1, 0, 0,-1, 0, -40, 0, 0, 0 ; % 24
2, 1, 0,-2, 0, 10, 0, 0, 0 ; % 25
0, 0,-2, 2, 1, 10, 0, 0, 0 ; % 26
0, 1,-2, 2, 0, -10, 0, 0, 0 ; % 27
0, 1, 0, 0, 2, 10, 0, 0, 0 ; % 28
-1, 0, 0, 1, 1, 10, 0, 0, 0 ; % 29
0, 1, 2,-2, 0, -10, 0, 0, 0 ; % 30
0, 0, 2, 0, 2, -22740, -2, 9770, -5 ; % 31
1, 0, 0, 0, 0, 7120, 1, -70, 0 ; % 32
0, 0, 2, 0, 1, -3860, -4, 2000, 0 ; % 33
1, 0, 2, 0, 2, -3010, 0, 1290, -1 ; % 34
1, 0, 0,-2, 0, -1580, 0, -10, 0 ; % 35
-1, 0, 2, 0, 2, 1230, 0, -530, 0 ; % 36
0, 0, 0, 2, 0, 630, 0, -20, 0 ; % 37
1, 0, 0, 0, 1, 630, 1, -330, 0 ; % 38
-1, 0, 0, 0, 1, -580, -1, 320, 0 ; % 39
-1, 0, 2, 2, 2, -590, 0, 260, 0 ; % 40
1, 0, 2, 0, 1, -510, 0, 270, 0 ; % 41
0, 0, 2, 2, 2, -380, 0, 160, 0 ; % 42
2, 0, 0, 0, 0, 290, 0, -10, 0 ; % 43
1, 0, 2,-2, 2, 290, 0, -120, 0 ; % 44
2, 0, 2, 0, 2, -310, 0, 130, 0 ; % 45
0, 0, 2, 0, 0, 260, 0, -10, 0 ; % 46
-1, 0, 2, 0, 1, 210, 0, -100, 0 ; % 47
-1, 0, 0, 2, 1, 160, 0, -80, 0 ; % 48
1, 0, 0,-2, 1, -130, 0, 70, 0 ; % 49
-1, 0, 2, 2, 1, -100, 0, 50, 0 ; % 50
1, 1, 0,-2, 0, -70, 0, 0, 0 ; % 51
0, 1, 2, 0, 2, 70, 0, -30, 0 ; % 52
0,-1, 2, 0, 2, -70, 0, 30, 0 ; % 53
1, 0, 2, 2, 2, -80, 0, 30, 0 ; % 54
1, 0, 0, 2, 0, 60, 0, 0, 0 ; % 55
2, 0, 2,-2, 2, 60, 0, -30, 0 ; % 56
0, 0, 0, 2, 1, -60, 0, 30, 0 ; % 57
0, 0, 2, 2, 1, -70, 0, 30, 0 ; % 58
1, 0, 2,-2, 1, 60, 0, -30, 0 ; % 59
0, 0, 0,-2, 1, -50, 0, 30, 0 ; % 60
1,-1, 0, 0, 0, 50, 0, 0, 0 ; % 61
2, 0, 2, 0, 1, -50, 0, 30, 0 ; % 62
0, 1, 0,-2, 0, -40, 0, 0, 0 ; % 63
1, 0,-2, 0, 0, 40, 0, 0, 0 ; % 64
0, 0, 0, 1, 0, -40, 0, 0, 0 ; % 65
1, 1, 0, 0, 0, -30, 0, 0, 0 ; % 66
1, 0, 2, 0, 0, 30, 0, 0, 0 ; % 67
1,-1, 2, 0, 2, -30, 0, 10, 0 ; % 68
-1,-1, 2, 2, 2, -30, 0, 10, 0 ; % 69
-2, 0, 0, 0, 1, -20, 0, 10, 0 ; % 70
3, 0, 2, 0, 2, -30, 0, 10, 0 ; % 71
0,-1, 2, 2, 2, -30, 0, 10, 0 ; % 72
1, 1, 2, 0, 2, 20, 0, -10, 0 ; % 73
-1, 0, 2,-2, 1, -20, 0, 10, 0 ; % 74
2, 0, 0, 0, 1, 20, 0, -10, 0 ; % 75
1, 0, 0, 0, 2, -20, 0, 10, 0 ; % 76
3, 0, 0, 0, 0, 20, 0, 0, 0 ; % 77
0, 0, 2, 1, 2, 20, 0, -10, 0 ; % 78
-1, 0, 0, 0, 2, 10, 0, -10, 0 ; % 79
1, 0, 0,-4, 0, -10, 0, 0, 0 ; % 80
-2, 0, 2, 2, 2, 10, 0, -10, 0 ; % 81
-1, 0, 2, 4, 2, -20, 0, 10, 0 ; % 82
2, 0, 0,-4, 0, -10, 0, 0, 0 ; % 83
1, 1, 2,-2, 2, 10, 0, -10, 0 ; % 84
1, 0, 2, 2, 1, -10, 0, 10, 0 ; % 85
-2, 0, 2, 4, 2, -10, 0, 10, 0 ; % 86
-1, 0, 4, 0, 2, 10, 0, 0, 0 ; % 87
1,-1, 0,-2, 0, 10, 0, 0, 0 ; % 88
2, 0, 2,-2, 1, 10, 0, -10, 0 ; % 89
2, 0, 2, 2, 2, -10, 0, 0, 0 ; % 90
1, 0, 0, 2, 1, -10, 0, 0, 0 ; % 91
0, 0, 4,-2, 2, 10, 0, 0, 0 ; % 92
3, 0, 2,-2, 2, 10, 0, 0, 0 ; % 93
1, 0, 2,-2, 0, -10, 0, 0, 0 ; % 94
0, 1, 2, 0, 1, 10, 0, 0, 0 ; % 95
-1,-1, 0, 2, 1, 10, 0, 0, 0 ; % 96
0, 0,-2, 0, 1, -10, 0, 0, 0 ; % 97
0, 0, 2,-1, 2, -10, 0, 0, 0 ; % 98
0, 1, 0, 2, 0, -10, 0, 0, 0 ; % 99
1, 0,-2,-2, 0, -10, 0, 0, 0 ; % 100
0,-1, 2, 0, 1, -10, 0, 0, 0 ; % 101
1, 1, 0,-2, 1, -10, 0, 0, 0 ; % 102
1, 0,-2, 2, 0, -10, 0, 0, 0 ; % 103
2, 0, 0, 2, 0, 10, 0, 0, 0 ; % 104
0, 0, 2, 4, 2, -10, 0, 0, 0 ; % 105
0, 1, 0, 1, 0, 10, 0, 0, 0 % 106
];
% Mean arguments of luni-solar motion
% l mean anomaly of the Moon
% l' mean anomaly of the Sun
% F mean argument of latitude
% D mean longitude elongation of the Moon from the Sun
% Om mean longitude of the ascending node
l = mod ( 485866.733 + (1325.0*rev + 715922.633)*T ...
+ 31.310*T2 + 0.064*T3, rev );
lp = mod ( 1287099.804 + ( 99.0*rev + 1292581.224)*T ...
- 0.577*T2 - 0.012*T3, rev );
F = mod ( 335778.877 + (1342.0*rev + 295263.137)*T ...
- 13.257*T2 + 0.011*T3, rev );
D = mod ( 1072261.307 + (1236.0*rev + 1105601.328)*T ...
- 6.891*T2 + 0.019*T3, rev );
Om = mod ( 450160.280 - ( 5.0*rev + 482890.539)*T ...
+ 7.455*T2 + 0.008*T3, rev );
% Nutation in longitude and obliquity [rad]
dpsi = 0;
deps = 0;
for i=1:N_coeff
arg = ( C(i,1)*l+C(i,2)*lp+C(i,3)*F+C(i,4)*D+C(i,5)*Om )/const.Arcs;
dpsi = dpsi + ( C(i,6)+C(i,7)*T ) * sin(arg);
deps = deps + ( C(i,8)+C(i,9)*T ) * cos(arg);
end
dpsi = 1e-5 * dpsi/const.Arcs;
deps = 1e-5 * deps/const.Arcs;
格林尼治真地方时(GAST)与格林尼治平地方时(GMST)的关系为
%--------------------------------------------------------------------------
%
% gast: Greenwich Apparent Sidereal Time
%
% Input:
% Mjd_UT1 Modified Julian Date UT1
%
% Output:
% gstime GAST in [rad]
%
% Last modified: 2018/01/04 M. Mahooti
%
%--------------------------------------------------------------------------
function gstime = gast(Mjd_UT1)
global const
gstime = mod(gmst(Mjd_UT1) + EqnEquinox(Mjd_UT1), const.pi2);
%--------------------------------------------------------------------------
%
% EqnEquinox: Computation of the equation of the equinoxes
%
% Input:
% Mjd_TT Modified Julian Date (Terrestrial Time)
%
% Output:
% EqE Equation of the equinoxes
%
% Notes:
% The equation of the equinoxes dpsi*cos(eps) is the right ascension of
% the mean equinox referred to the true equator and equinox and is equal
% to the difference between apparent and mean sidereal time.
%
% Last modified: 2015/08/12 M. Mahooti
%
%--------------------------------------------------------------------------
function EqE = EqnEquinox (Mjd_TT)
% Nutation in longitude and obliquity
[dpsi, deps] = NutAngles (Mjd_TT);
% Equation of the equinoxes
EqE = dpsi * cos(MeanObliquity(Mjd_TT));
由于行星的进动,黄道的倾角略有下降,相当于
代码如下:
%--------------------------------------------------------------------------
%
% MeanObliquity: Computes the mean obliquity of the ecliptic
%
% Input:
% Mjd_TT Modified Julian Date (Terrestrial Time)
%
% Output:
% MOblq Mean obliquity of the ecliptic
%
% Last modified: 2018/01/04 M. Mahooti
%
%--------------------------------------------------------------------------
function MOblq = MeanObliquity(Mjd_TT)
global const
T = (Mjd_TT-const.MJD_J2000)/36525;
MOblq = const.Rad*(23.43929111-(46.8150+(0.00059-0.001813*T)*T)*T/3600);
最后得到了真实日期的坐标系(采用的进动章动理论定义)与地球赤道和格林尼治子午线对齐的系统之间的转换矩阵
代码如下:
%--------------------------------------------------------------------------
%
% GHAMatrix: Transformation from true equator and equinox to Earth equator
% and Greenwich meridian system
%
% Input:
% Mjd_UT1 Modified Julian Date UT1
%
% Output:
% GHAmat Greenwich Hour Angle matrix
%
% Last modified: 2015/08/12 M. Mahooti
%
%--------------------------------------------------------------------------
function GHAmat = GHAMatrix(Mjd_UT1)
GHAmat = R_z(gast(Mjd_UT1));
%--------------------------------------------------------------------------
% Input:
% angle angle of rotation [rad]
%
% Output:
% rotmat rotation matrix
%--------------------------------------------------------------------------
function [rotmat] = R_z(angle)
C = cos(angle);
S = sin(angle);
rotmat = zeros(3,3);
rotmat(1,1) = C; rotmat(1,2) = S; rotmat(1,3) = 0.0;
rotmat(2,1) = -1.0*S; rotmat(2,2) = C; rotmat(2,3) = 0.0;
rotmat(3,1) = 0.0; rotmat(3,2) = 0.0; rotmat(3,3) = 1.0;