forked from poliastro/validation
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ecef2teme.m
92 lines (78 loc) · 3.23 KB
/
ecef2teme.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
% ----------------------------------------------------------------------------
%
% function ecef2teme
%
% this function trsnforms a vector from the earth fixed (ITRF) frame to the
% true equator mean equniox frame (teme). the results take into account
% the effects of sidereal time, and polar motion.
%
% author : david vallado 719-573-2600 30 oct 2017
%
% revisions
%
% inputs description range / units
% recef - position vector earth fixed km
% vecef - velocity vector earth fixed km/s
% aecef - acceleration vector earth fixedkm/s2
% ttt - julian centuries of tt centuries
% jdut1 - julian date of ut1 days from 4713 bc
% lod - excess length of day sec
% xp - polar motion coefficient arc sec
% yp - polar motion coefficient arc sec
% eqeterms - use extra two terms (kinematic) after 1997 0, 2
%
% outputs :
% rteme - position vector teme km
% vteme - velocity vector teme km/s
% ateme - acceleration vector teme km/s2
%
% locals :
% st - matrix for pef - tod
% pm - matrix for ecef - pef
%
% coupling :
% gstime - greenwich mean sidereal time rad
% polarm - rotation for polar motion pef - ecef
%
% references :
% vallado 2013, 231-233
%
% [rteme, vteme, ateme] = ecef2teme( recef, vecef, aecef, ttt, jdut1, lod, xp, yp, eqeterms );
% ----------------------------------------------------------------------------
function [rteme, vteme, ateme] = ecef2teme( recef, vecef, aecef, ttt, jdut1, lod, xp, yp, eqeterms );
deg2rad = pi/180.0;
% ------------------------ find gmst --------------------------
gmst= gstime( jdut1 );
% find omeage from nutation theory
omega= 125.04452222 + ( -6962890.5390 *ttt + ...
7.455 *ttt*ttt + 0.008 *ttt*ttt*ttt ) / 3600.0;
omega= rem( omega, 360.0 ) * deg2rad;
% teme does not include the geometric terms here
% after 1997, kinematic terms apply
if (jdut1 > 2450449.5 ) && (eqeterms > 0)
gmstg = gmst ...
+ 0.00264*pi /(3600*180)*sin(omega) ...
+ 0.000063*pi /(3600*180)*sin(2.0 *omega);
else
gmstg = gmst;
end
gmstg = rem (gmstg, 2.0*pi);
st(1,1) = cos(gmstg);
st(1,2) = -sin(gmstg);
st(1,3) = 0.0;
st(2,1) = sin(gmstg);
st(2,2) = cos(gmstg);
st(2,3) = 0.0;
st(3,1) = 0.0;
st(3,2) = 0.0;
st(3,3) = 1.0;
[pm] = polarm(xp, yp, ttt, '80');
rpef = pm*recef;
rteme = st*rpef;
thetasa= 7.29211514670698e-05 * (1.0 - lod/86400.0 );
omegaearth = [0; 0; thetasa;];
vpef = pm*vecef;
vteme = st*(vpef + cross(omegaearth,rpef));
temp = cross(omegaearth,rpef);
ateme = st*( pm*aecef + cross(omegaearth,temp) ...
+ 2.0*cross(omegaearth,vpef) );