Skip to content

Commit 20ece60

Browse files
scivisionscivision
committed
ecef2eci, eci2ecef: much more accurate algorithms
Interface for these two functions is now like Matlab.
1 parent d2caeb2 commit 20ece60

27 files changed

Lines changed: 833 additions & 34 deletions

‎+matmap3d/aer2eci.m‎

Lines changed: 15 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
%
88
%%% Outputs
99
% * x, y, z: ECI x, y, z
10-
function [x, y, z] = aer2eci(utc, az, el, rng, lat, lon, alt)
10+
function [x, y, z] = aer2eci(utc, az, el, rng, lat, lon, alt, spheroid, angleUnit)
1111
arguments
1212
utc datetime
1313
az {mustBeReal}
@@ -16,10 +16,22 @@
1616
lat {mustBeReal}
1717
lon {mustBeReal}
1818
alt {mustBeReal}
19+
spheroid (1,1) matmap3d.referenceEllipsoid = matmap3d.wgs84Ellipsoid()
20+
angleUnit (1,1) string = "d"
1921
end
2022

21-
[x1, y1, z1] = matmap3d.aer2ecef(az, el, rng, lat, lon, alt);
23+
[x1, y1, z1] = matmap3d.aer2ecef(az, el, rng, lat, lon, alt, spheroid, angleUnit);
2224

23-
[x, y, z] = matmap3d.ecef2eci(utc, x1, y1, z1);
25+
x = nan(like=az);
26+
y = nan(like=az);
27+
z = nan(like=az);
28+
29+
for i = 1:numel(x)
30+
r_eci = matmap3d.ecef2eci(utc, [x1(i); y1(i); z1(i)]);
31+
32+
x(i) = r_eci(1);
33+
y(i) = r_eci(2);
34+
z(i) = r_eci(3);
35+
end
2436

2537
end

‎+matmap3d/ecef2eci.m‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
%% ECEF2ECI: Transforms Earth Centered Earth Fixed (ECEF) coordinates to
2+
% Earth Centered Inertial (ECI) coordinates
3+
%
4+
%%% Inputs
5+
% * utc: datetime UTC
6+
% * r_ecef: 3x1 vector of ECEF position [meters]
7+
% * v_ecef: 3x1 vector if ECEF velocity [m/s] (optional)
8+
%%% Outputs
9+
% * r_eci: 3x1 vector of ECI position [meters]
10+
% * v_eci; 3x1 vector of ECI velocity [m/s] (optional)
11+
%
12+
% Original by: Meysam Mahooti
13+
14+
function [r_eci, v_eci] = ecef2eci(utc, r_ecef, v_ecef, eopdata)
15+
arguments
16+
utc (1,1) datetime
17+
r_ecef (3,1) {mustBeReal}
18+
v_ecef {mustBeReal} = []
19+
eopdata (13, :) {mustBeReal} = h5read(fullfile(fileparts(mfilename('fullpath')), "private/EOP-All.h5"), '/eop')
20+
end
21+
22+
mjd_utc = Mjday(utc);
23+
24+
[x_pole,y_pole,UT1_UTC,LOD,~,~,~,~,TAI_UTC] = IERS(eopdata, mjd_utc, 'l');
25+
[~,~,~,TT_UTC] = timediff(UT1_UTC,TAI_UTC);
26+
27+
MJD_UT1 = mjd_utc + UT1_UTC/86400;
28+
MJD_TT = mjd_utc + TT_UTC/86400;
29+
30+
% ICRS to ITRS transformation matrix and its derivative
31+
P = PrecMatrix(51544.5, MJD_TT); % IAU 1976 Precession
32+
N = NutMatrix(MJD_TT); % IAU 1980 Nutation
33+
Theta = GHAMatrix(MJD_UT1,MJD_TT); % Earth rotation
34+
Po = PoleMatrix(x_pole,y_pole); % Polar motion
35+
U = Po*Theta*N*P; % ICRS to ITRS transformation
36+
37+
% Transformation from WGS to ICRS
38+
r_eci = U'*r_ecef;
39+
if nargout > 1
40+
mustBeEqualSize(r_ecef, v_ecef)
41+
S = zeros(3);
42+
S(1,2) = 1;
43+
S(2,1) = -1;
44+
Omega = 15.04106717866910/3600*pi/180 - 0.843994809*1e-9*LOD; % [rad/s] IERS
45+
dTheta = Omega*S*Theta; % Derivative of Earth rotation matrix [1/s]
46+
dU = Po*dTheta*N*P; % Derivative [1/s]
47+
v_eci = U'*v_ecef + dU'*r_ecef;
48+
end
49+
50+
end

‎+matmap3d/eci2aer.m‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -18,8 +18,14 @@
1818
alt {mustBeReal}
1919
end
2020

21-
[x, y, z] = matmap3d.eci2ecef(utc, x0, y0, z0);
21+
az = nan(like=x0);
22+
el = nan(like=x0);
23+
rng = nan(like=x0);
2224

23-
[az, el, rng] = matmap3d.ecef2aer(x, y, z, lat, lon, alt);
25+
for i = 1:numel(x0)
26+
r_ecef = matmap3d.eci2ecef(utc, [x0(i); y0(i); z0(i)]);
27+
28+
[az(i), el(i), rng(i)] = matmap3d.ecef2aer(r_ecef(1), r_ecef(2), r_ecef(3), lat, lon, alt);
29+
end
2430

2531
end

‎+matmap3d/eci2ecef.m‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
%% ECI2ECEF: Transforms Earth Centered Inertial (ECI) coordinates to Earth
2+
% Centered Earth Fixed (ECEF) coordinates
3+
%
4+
% original by: Meysam Mahooti
5+
%
6+
7+
function [r_ecef, v_ecef] = eci2ecef(utc, r_eci, v_eci, eopdata)
8+
arguments
9+
utc (1,1) datetime
10+
r_eci (3,1) {mustBeReal}
11+
v_eci {mustBeReal} = []
12+
eopdata (13, :) {mustBeReal} = h5read(fullfile(fileparts(mfilename('fullpath')), "private/EOP-All.h5"), '/eop')
13+
end
14+
15+
mjd_utc = Mjday(utc);
16+
17+
[x_pole,y_pole,UT1_UTC,LOD,~,~,~,~,TAI_UTC] = IERS(eopdata, mjd_utc,'l');
18+
[~,~,~,TT_UTC] = timediff(UT1_UTC,TAI_UTC);
19+
20+
MJD_UT1 = mjd_utc + UT1_UTC/86400;
21+
MJD_TT = mjd_utc + TT_UTC/86400;
22+
23+
% ICRS to ITRS transformation matrix and its derivative
24+
P = PrecMatrix(51544.5, MJD_TT); % IAU 1976 Precession
25+
N = NutMatrix(MJD_TT); % IAU 1980 Nutation
26+
Theta = GHAMatrix(MJD_UT1,MJD_TT); % Earth rotation
27+
Po = PoleMatrix(x_pole,y_pole); % Polar motion
28+
29+
S = zeros(3);
30+
S(1,2) = 1;
31+
S(2,1) = -1;
32+
Omega = 15.04106717866910/3600*pi/180 - 0.843994809*1e-9*LOD; % [rad/s] IERS
33+
dTheta = Omega*S*Theta; % Derivative of Earth rotation matrix [1/s]
34+
U = Po*Theta*N*P; % ICRS to ITRS transformation
35+
dU = Po*dTheta*N*P; % Derivative [1/s] % Derivative [1/s]
36+
37+
% Transformation from ICRS to ITRS
38+
r_ecef = U * r_eci;
39+
if nargout > 1
40+
mustBeEqualSize(r_eci, v_eci)
41+
v_ecef = U * v_eci + dU * r_eci;
42+
end
43+
44+
end

‎+matmap3d/private/EOP-All.h5‎

1.1 MB
Binary file not shown.

‎+matmap3d/private/EqnEquinox.m‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
%% EQNEQUINOX Computation of the equation of the equinoxes
2+
%
3+
% Input:
4+
% Mjd_TT Modified Julian Date (Terrestrial Time)
5+
%
6+
% Output:
7+
% EqE Equation of the equinoxes
8+
%
9+
% Notes:
10+
% The equation of the equinoxes dpsi*cos(eps) is the right ascension of
11+
% the mean equinox referred to the true equator and equinox and is equal
12+
% to the difference between apparent and mean sidereal time.
13+
%
14+
% Last modified: 2018/01/27 Meysam Mahooti
15+
%
16+
% Meysam Mahooti (2025).
17+
% ECI2ECEF & ECEF2ECI Transformations
18+
% https://www.mathworks.com/matlabcentral/fileexchange/61957-eci2ecef-ecef2eci-transformations)
19+
% MATLAB Central File Exchange.
20+
21+
function EqE = EqnEquinox (Mjd_TT)
22+
23+
% Nutation in longitude and obliquity
24+
dpsi = NutAngles (Mjd_TT);
25+
26+
% Equation of the equinoxes
27+
EqE = dpsi * cos(MeanObliquity(Mjd_TT));
28+
29+
end

‎+matmap3d/private/Frac.m‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
%% FRAC Fractional part of a number (y=x-[x])
2+
%
3+
% Last modified: 2018/01/27 Meysam Mahooti
4+
%
5+
% Meysam Mahooti (2025).
6+
% ECI2ECEF & ECEF2ECI Transformations
7+
% https://www.mathworks.com/matlabcentral/fileexchange/61957-eci2ecef-ecef2eci-transformations)
8+
% MATLAB Central File Exchange.
9+
10+
function [res] = Frac(x)
11+
12+
res = x-floor(x);
13+
14+
end

‎+matmap3d/private/GHAMatrix.m‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
%% GHAMATRIX Transformation from true equator and equinox to Earth equator and Greenwich meridian system
2+
%
3+
% Inputs:
4+
% Mjd_UT1 Modified Julian Date UT1
5+
% Mjd_TT Modified Julian Date TT
6+
%
7+
% Output:
8+
% GHAmat Greenwich Hour Angle matrix
9+
%
10+
% Last modified: 2022/06/16 Meysam Mahooti
11+
%
12+
% Meysam Mahooti (2025).
13+
% ECI2ECEF & ECEF2ECI Transformations
14+
% https://www.mathworks.com/matlabcentral/fileexchange/61957-eci2ecef-ecef2eci-transformations)
15+
% MATLAB Central File Exchange.
16+
17+
function GHAmat = GHAMatrix(Mjd_UT1,Mjd_TT)
18+
19+
GHAmat = R_z(gast(Mjd_UT1,Mjd_TT));
20+
21+
end

‎+matmap3d/private/IERS.m‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
%% IERS Management of IERS time and polar motion data
2+
%
3+
% Last modified: 2018/02/01 Meysam Mahooti
4+
%
5+
% Meysam Mahooti (2025).
6+
% ECI2ECEF & ECEF2ECI Transformations
7+
% https://www.mathworks.com/matlabcentral/fileexchange/61957-eci2ecef-ecef2eci-transformations)
8+
% MATLAB Central File Exchange.
9+
10+
function [x_pole,y_pole,UT1_UTC,LOD,dpsi,deps,dx_pole,dy_pole,TAI_UTC] = IERS(eop,Mjd_UTC,interp)
11+
arguments
12+
eop (13,:) {mustBeReal}
13+
Mjd_UTC (1,1) {mustBeReal, mustBePositive, mustBeFinite}
14+
interp {mustBeTextScalar} = 'n'
15+
end
16+
17+
Arcs = 3600*180/pi;
18+
19+
mjd = floor(Mjd_UTC);
20+
i = find(mjd==eop(4,:),1,'first');
21+
22+
switch interp
23+
case 'l'
24+
% linear interpolation
25+
preeop = eop(:,i);
26+
nexteop = eop(:,i+1);
27+
fixf = Mjd_UTC-floor(Mjd_UTC);
28+
% Setting of IERS Earth rotation parameters
29+
% (UT1-UTC [s], TAI-UTC [s], x ["], y ["])
30+
x_pole = preeop(5)+(nexteop(5)-preeop(5))*fixf;
31+
y_pole = preeop(6)+(nexteop(6)-preeop(6))*fixf;
32+
UT1_UTC = preeop(7)+(nexteop(7)-preeop(7))*fixf;
33+
LOD = preeop(8)+(nexteop(8)-preeop(8))*fixf;
34+
dpsi = preeop(9)+(nexteop(9)-preeop(9))*fixf;
35+
deps = preeop(10)+(nexteop(10)-preeop(10))*fixf;
36+
dx_pole = preeop(11)+(nexteop(11)-preeop(11))*fixf;
37+
dy_pole = preeop(12)+(nexteop(12)-preeop(12))*fixf;
38+
TAI_UTC = preeop(13);
39+
40+
x_pole = x_pole/Arcs; % Pole coordinate [rad]
41+
y_pole = y_pole/Arcs; % Pole coordinate [rad]
42+
dpsi = dpsi/Arcs;
43+
deps = deps/Arcs;
44+
dx_pole = dx_pole/Arcs; % Pole coordinate [rad]
45+
dy_pole = dy_pole/Arcs; % Pole coordinate [rad]
46+
case 'n'
47+
eop = eop(:,i);
48+
% Setting of IERS Earth rotation parameters
49+
% (UT1-UTC [s], TAI-UTC [s], x ["], y ["])
50+
x_pole = eop(5)/Arcs; % Pole coordinate [rad]
51+
y_pole = eop(6)/Arcs; % Pole coordinate [rad]
52+
UT1_UTC = eop(7); % UT1-UTC time difference [s]
53+
LOD = eop(8); % Length of day [s]
54+
dpsi = eop(9)/Arcs;
55+
deps = eop(10)/Arcs;
56+
dx_pole = eop(11)/Arcs; % Pole coordinate [rad]
57+
dy_pole = eop(12)/Arcs; % Pole coordinate [rad]
58+
TAI_UTC = eop(13); % TAI-UTC time difference [s]
59+
end
60+
61+
end

‎+matmap3d/private/MeanObliquity.m‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
%% MEANOBLIQUITY Computes the mean obliquity of the ecliptic
2+
%
3+
% Input:
4+
% Mjd_TT Modified Julian Date (Terrestrial Time)
5+
%
6+
% Output:
7+
% MOblq Mean obliquity of the ecliptic
8+
%
9+
% Last modified: 2018/01/27 Meysam Mahooti
10+
%
11+
% Meysam Mahooti (2025).
12+
% ECI2ECEF & ECEF2ECI Transformations
13+
% https://www.mathworks.com/matlabcentral/fileexchange/61957-eci2ecef-ecef2eci-transformations)
14+
% MATLAB Central File Exchange.
15+
16+
function MOblq = MeanObliquity(Mjd_TT)
17+
18+
T = (Mjd_TT - 51544.5)/36525;
19+
20+
MOblq = (pi/180)*(23.43929111-(46.8150+(0.00059-0.001813*T)*T)*T/3600);
21+
22+
end

0 commit comments

Comments
 (0)