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)');
采纳的回答
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
[rotationAng1 rotationAng2 rotationAng3] = dcm2angle(C_eci2ecef); % default ZYX sequence
[rotationAng1 rotationAng2 rotationAng3]*180/pi
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 个评论
Derrick Early
2023-11-27
编辑：Derrick Early
2023-11-27
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))
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
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
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
2023-11-28
编辑：Derrick Early
2023-11-28
更多回答（1 个）
Ayush Anand
2023-11-23
编辑：Ayush Anand
2023-11-23
1 个投票
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 个评论
Derrick Early
2023-11-24
编辑：Derrick Early
2023-11-24
Derrick Early
2023-11-27
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
2023-11-27
编辑：Derrick Early
2023-11-27
类别
在 帮助中心 和 File Exchange 中查找有关 Physical and Time Unit Conversions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
