%% eci2geo.m converts from Earth-Centered Inertial (ECI) 
%% co-ordinates to geographic co-ordinates. The equations
%% used are from the Satellite Times articles by 
%% Dr T.S.Kelso.

function [lonmat,latmat,dismat]=eci2geo(xi,yi,zi,timi,ts)

% xi - ECI x co-ordinate (km)
% yi - ECI y co-ordinate (km)
% zi - ECI z co-ordinate (km)
% timi - timi==1 time since JD 2451545.0 (1200 UT, 1 Jan 2000) (sec) 
%        otherwise Julian date
% ts - timi switch 
% lonmat - geographic longitude (rads)
% latmat - geographic longitude (rads)
% dismat - radial distance (km)

re=6371;                      %Earth radi (km)
we=7.29211510e-5;             %Earths rotation speed (rads/s)

if ts==1
 jd=2451545.0+timi/86400;     %required Julian date
else
 jd=timi;                     %required Julian date
end

%Calculate Greenwich Mean Sideral Time (GMST) 
timf=rem(jd+0.5,1);                                                   %UT of measurement (fraction of day)
timmat=24*timf;                                                       %UT of measurement (hours)
tu=(jd-timf-2451545.0)/36525;
gmst0=24110.54841+8640184.812866*tu+0.093104*tu.^2-6.2e-6*tu.^3;      %GMST at 0UT (sec)
gmst0r=mod(gmst0,86400)*(2*pi/86400);                                 %GMST at 0UT (rads)     
gmst=gmst0r+we*3600*timmat;                                           %GMST (rads)     

%Calculate the required geographic location
latmat=atan2(zi,sqrt(xi.^2+yi.^2));          %latitude (rads)
lonmat=atan2(yi,xi)-gmst;
lonmat=longch(lonmat);                       %longitude (rads)
dismat=sqrt(xi.^2+yi.^2+zi.^2);              %radial distance (km)
