lla2eci conversion for geo sat

When converting a geosynchronous orbit in lla to eci, I expected the eci z component to be zero. I used the IAU-2000 epoch of January 1, 2000, 12:00 TT (Terrestrial Time).

rg=42164172; % geo radius (m)
re=6378137;  % earth radius (m)

startTime = datetime(2000,1,1,12,0,0);
stopTime = startTime + days(1);
sampleTime = 60;                                      % seconds
t=(startTime:seconds(sampleTime):stopTime)';
nt=length(t);
lla=[zeros(nt,2) (rg-re)*ones(nt,1)];

utc=datevec(t);
p=lla2eci(lla,utc);

plot(t,p(:,3)); ylabel('eci z (m)');

Answer 1

[Paul]

Hi Derrick,

I think the effect we're seeing here is a slight angular offset betweent the ECI z-axis and the Earth's axis of rotation. After 24 hours, the sattelite comes back to its initial position in ECI (even thought the plot is labeled "ecef z", so I assume not exactly generated from the code as shown).

Also, keep in mind that Terrestrial Time (TT) is offset from UTC by a few 10's of seconds. Unfortunately, the Wikipedia pages I've been reading seem to have a discrepancy; I've come up with either a 64.184 or 69.184 offset depending on where I look. So don't know what to make of that.

I'm not quite sure what the Aerospace Toolbox is defining as "the" ECI frame. I didn't see anything on point tin the documentation. This doc page suggests that ECI might be J2000, but it doesn't say so definitively.

However, the ECI frame seems to be quite different than J2000. Assuming that TT is the same as UTC, we'd have the ECI to ECEF DCM at the epoch be

T0 = datetime(2000,1,1,12,0,0);
C_eci2ecef = dcmeci2ecef('IAU-2000/2006',T0); %using defaults for everything else

C_eci2ecef = 3×3

0.1816 -0.9834 -0.0000 0.9834 0.1816 0.0000 -0.0000 -0.0000 1.0000

[rotationAng1 rotationAng2 rotationAng3] = dcm2angle(C_eci2ecef);  % default ZYX sequence
[rotationAng1 rotationAng2 rotationAng3]*180/pi

ans = 1×3

-79.5394 0.0013 0.0018

So, at (or near) the (assumed) epoch, ECEF is off by nearly 80 degrees from being in alignment with ECI, which is a much bigger difference than would result by accounting for the difference between TT and UTC. Maybe I'm not using dcmeci2ecef correctly? Or maybe ECI is defined differently than assumed by T0. IDK.

The Aerospace Toolbox often confuses me, which is why I'd be hesitant to use it in anger should I ever have the need to do so. It would be nice if it offered a simple ECEF model where the Earth's angular velocity vector is constant in ECI, i.e., no precession or nutation, etc. which I think is modeled, and that ECEF can be defined as being exactly in aligment with ECI at some arbitrary time.

Disclaimer: I don't claim to have any expertise in the detailed earth rotation modeling as implemented in the Aerospace Toolbox.

[Derrick Early]

Looking at dcmeci2ecef, line 223 uses J2000. Apparently, the International Earth Rotation and Reference Systems Service (IERS) never intended for eci and ecef to be collinear at J2000. They start out with a earth rotation angle of about 280 deg.

IERS Technical Note No. 36 page 18.

[Paul]

Thank you for bringing that note to my attention.

According to Earth-centered inertial - Wikipedia

• J2000: One commonly used ECI frame is defined with the Earth's Mean Equator and Mean Equinox (MEME) at 12:00 Terrestrial Time on 1 January 2000. It can be referred to as J2K, J2000 or EME2000. The x-axis is aligned with the mean vernal equinox. The z-axis is aligned with the Earth's rotation axis (or equivalently, the celestial North Pole) as it was at that time. The y-axis is rotated by 90° East about the celestial equator.[5]"

Because the vernal equinox (March 20, 2000) is 79 days later than 1 Jan 2000, that would account for the -80 deg (or 280 deg) difference between ECI and ECEF at TT0.

A little trial and error shows that the time of alignment (i.e., rotation about z is 0 deg) is a little bit after 2000:03:20:12:06:40 UTC, and we can get the seconds more precisely from

T0 = fsolve(@(T0) dcm2angle(dcmeci2ecef('IAU-2000/2006',[2000 03 20 12 6 T0]))*180/pi,40,optimoptions('fsolve','TolFun',1e-10))

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

T0 = 41.1442

which yields a very small rotation about z (-1e-8 deg) from ECI to ECEF

format short e
startTime = datetime(2000,03,20,12,6,T0);
[rot1,rot2,rot3] = dcm2angle(dcmeci2ecef('IAU-2000/2006',datevec(startTime)));
[rot1,rot2,rot3]*180/pi

ans = 1×3

-1.1670e-08 -4.3063e-04 1.1134e-03

Recreating the plot from above ...

format long e
rg = 42164172; % geo radius (m)
re = 6378137;  % earth radius (m)
startTime        = datetime(2000,1,1,12,0,0);
stopTime   = startTime + days(1);
sampleTime = 60;                                      % seconds
t          = (startTime:seconds(sampleTime):stopTime)';
nt         = length(t);
lla        = [zeros(nt,2) (rg-re)*ones(nt,1)];
utc        = datevec(t);
p          = lla2eci(lla,utc);
figure
plot(t,p(:,3)); ylabel('eci z (m)'),grid

The other part of that statement is that the ECI z-axis is aligned with the earth's rotation axis "as it was" at TT0. However, the earth's rotation axis is not aligned with the ECEF z-axis (Earth-centered, Earth-fixed coordinate system - Wikipedia). That misalignment is apparently what's captured in rot2 and rot3

[rot1,rot2,rot3] = dcm2angle(dcmeci2ecef('IAU-2000/2006',datevec(startTime)));
[rot1,rot2,rot3]*180/pi

ans = 1×3

-7.953938102383586e+01 1.297575354190418e-03 1.809570371256625e-03

Assuming the misalignment is constant over the 24 hours, we'd expect the max and min deviation of the z position in ECI to be:

yline(rg*norm([rot2 rot3]))
yline(-rg*norm([rot2 rot3]))

That same ECEF wikipedia link says that this angular offset is due to "polar motion." However, the default for polarmotion in lla2eci and dcmeci2ecef is 0, so I'm still not sure what the source of that offset is in those functions as used above.

[Derrick Early]

Hi Paul,

This is excellent. Thank you.

Looking into dcmeci2ecef more, I found that the nutation rotation matrix, Q, is not identity at noon on 1 January 2000.

K>> [rot1,rot2,rot3]=dcm2angle(Q);

[rot1,rot2,rot3]*180/pi

ans =

0.0000 -0.0015 0.0016

In section 5.3.3, TN36 has the following: "In order to follow IAU 2000 Resolution B1.7, the realized celestial pole must be the [Celestial Intermediate Pole (CIP)]. This requires an offset at epoch in the conventional model for precession nutation as well as diurnal and higher frequency variations in the Earth's orientation. According to this resolution, the direction of the CIP at J2000.0 has to be offset from the pole of the GCRS in a manner consistent with the IAU 2006/2000A precession-nutation model."

Looks like the nutation offset at epoch is by design.

Answer 2

[Ayush Anand]

Hi Derrick,

I understand you are converting geosynchronous orbits from LLA to ECI coordinates, and expect the z-component of the ECI coordinates to be zero. However, this is not the case. The geosynchronous orbit is synchronized with the Earth's rotation and remains fixed relative to a point on the Earth's equator, but it does not mean that the orbit lies in the equatorial plane (which would give a zero ECI z-component).

The ECI coordinates are determined relative to the center of the Earth in a coordinate system which isn't rotating with the Earth. For a satellite to be geosynchronous, it needs to oribit the Earth at the same rate as the Earth's rotation, effectively staying in a fixed position above a point on the equator. This means the x and y components, representing its location along the equatorial plane, might not change much, but the z-component representing the altitude or the distance from the Earth's center wouldn't typically be zero as the satellite would be orbiting above the equator, maintaing a certain altitude to stay in sync with the Earth's rotation.

I hope this helps!

[Derrick Early]

Thank you for your response!

With the following line, I forced the orbit onto the equatorial plane.

lla=[zeros(nt,2) (rg-re)*ones(nt,1)];  

I thought that the ecef and eci frames were coincident at the j2k epoch of 1 January 2000 at 12 noon.

[Paul]

Why would the ECI z-component be altitude if the sattelite is basically over the equator? The ECI z-axis goes (nearly) through the North pole. Also, if the sattelite is orbiting above equator, wouldn't that also mean that its orbit lies in the equatorial plane?

[Derrick Early]

I was thinking that the ECI z-axis and the ECEF z-axis were collinear at January 1, 2000, 12:00 TT. Yes, I was attempting to model a sattelite on the equitorial plane. If the two z-axes were collinear, then the z positions in ECI should be zero.

[Paul]

Hi Derrick,

I understood what you were trying to do. I should have been clear that my comment was motivated by statements in the Answer by Ayush.

[Derrick Early]

Hi Paul, Thank you for your time. Your responses were very helpful.