ecef <=> eci: keep _naive functions

Author: scivision

+matmap3d/aer2eci.m

@@ -1,8 +1,5 @@
 %% AER2ECI convert AER (azimuth, elevation, slant range) to ECI
 %
-% NOTE: because underlying ecef2eci() is rotation only, error can be order
-% 1..10%
-%
 %%% Inputs
 % * utc: datetime UTC
 % * az,el,rng: (degrees, meters)

+matmap3d/ecef2eci.m +matmap3d/ecef2eci_naive.m+matmap3d/ecef2eci.m renamed to +matmap3d/ecef2eci_naive.m

@@ -1,13 +1,13 @@
-%% ECEF2ECI rotate ECEF coordinates to ECI
-% because this doesn't account for nutation, etc. error is often > 1%
+%% ECEF2ECI_NAIVE rotate ECEF coordinates to ECI
+% because this doesn't account for nutation, etc. error is generally significant
 %
 %%% Inputs
 % x0, y0, z0:  ECEF position (meters)
 % utc: time UTC
 %%% Outputs
 % * x,y,z:  ECI position (meters)
 
-function [x,y,z] = ecef2eci(utc, x0, y0, z0)
+function [x,y,z] = ecef2eci_naive(utc, x0, y0, z0)
 arguments
   utc (:,1) datetime
   x0 (:,1) {mustBeReal,mustBeEqualSize(utc,x0)}
@@ -19,9 +19,9 @@
 gst = matmap3d.greenwichsrt(juliandate(utc));
 
 % Convert into ECEF
-x = nan(size(gst));
-y = nan(size(x));
-z = nan(size(x));
+x = nan(like=gst);
+y = nan(like=gst);
+z = nan(like=gst);
 
 for j = 1:length(x)
   eci = matmap3d.R3(gst(j)).' * [x0(j), y0(j), z0(j)].';

+matmap3d/eci2ecef.m +matmap3d/eci2ecef_naive.m+matmap3d/eci2ecef.m renamed to +matmap3d/eci2ecef_naive.m

@@ -1,11 +1,12 @@
-%% ECI2ECEF rotate ECI coordinates to ECEF
-% because this doesn't account for nutation, etc. error is often > 1%
+%% ECI2ECEF_NAIVE rotate ECI coordinates to ECEF
+% because this doesn't account for nutation, etc. error is generally significant
 %
 % x_eci, y_eci, z_eci:  eci position vectors
 % utc: Matlab datetime UTC
 %
 % x,y,z:  ECEF position (meters)
-function [x,y,z] = eci2ecef(utc, x_eci, y_eci, z_eci)
+
+function [x,y,z] = eci2ecef_naive(utc, x_eci, y_eci, z_eci)
 arguments
   utc (:,1) datetime
   x_eci (:,1) {mustBeReal,mustBeEqualSize(utc,x_eci)}
@@ -16,9 +17,9 @@
 % Greenwich hour angles (radians)
 gst = matmap3d.greenwichsrt(juliandate(utc));
 
-x = nan(size(x_eci));
-y = nan(size(x));
-z = nan(size(x));
+x = nan(like=x_eci);
+y = nan(like=y_eci);
+z = nan(like=z_eci);
 
 for j = 1:length(x)
   ecef = matmap3d.R3(gst(j)) * [x_eci(j), y_eci(j), z_eci(j)].';

.archive/juliantime.m

@@ -0,0 +1,26 @@
+%% juliandate  datetime to Julian time
+%
+% from D.Vallado Fundamentals of Astrodynamics and Applications p.187
+% and J. Meeus Astronomical Algorithms 1991 Eqn. 7.1 pg. 61
+%
+% parameters:
+% t: datetime vector
+%
+% results:
+% jd: julian date (days since Jan 1, 4173 BCE
+
+function jd = juliandate(t)
+
+[y, mon, d, h, m, s] = datevec(t);
+
+i = mon < 3;
+y(i) = y(i) - 1;
+mon(i) = mon(i) + 12;
+
+A = fix(y / 100.);
+B = 2 - A + fix(A / 4.);
+C = ((s / 60. + m) / 60. + h) / 24.;
+
+jd = fix(365.25 * (y + 4716)) + fix(30.6001 * (mon + 1)) + d + B - 1524.5 + C;
+
+end