lla2eci conversion for geo sat

4 次查看（过去 30 天）
Derrick Early
Derrick Early 2023-11-14
编辑： Derrick Early 2023-11-28
Ran in:
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)');

请先登录，再回答此问题。

采纳的回答

Paul
Paul 2023-11-26
Ran in:
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.
  3 个评论
显示 1更早的评论
Paul
Paul 2023-11-27
编辑：Paul 2023-11-28
Ran in:
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
Derrick Early 2023-11-28
编辑：Derrick Early 2023-11-28
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.

请先登录，再进行评论。

更多回答（1 个）

Ayush Anand
Ayush Anand 2023-11-23
编辑：Ayush Anand 2023-11-23
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!
  5 个评论
显示 3更早的评论
Paul
Paul 2023-11-27
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
Derrick Early 2023-11-27
编辑：Derrick Early 2023-11-27
Hi Paul, Thank you for your time. Your responses were very helpful.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Coordinate Systems 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by