%% satloc2.m calculates the azimuth, elevation and range of a point as 
%% observed from a second point on a sphere. For example the second point
%% could be a ground station on Earth or a LEO, and the first point could
%% be a GPS satellite.
%% The equations used are from the Satellite Times articles by Dr T.S.Kelso.

function [rg,azs,els]=satloc2(timmat,day,month,year,lont,latt,dist,lono,lato,diso)

% timmat - UT time (hours)
% day = day
% month = month
% year = year
% lont = longitude of target (rads)
% latt = latitude of target (rads)
% dist = radial distance of target (km)
% lono = longitude of obseving point (rads)
% lato = latitude of observing point (rads)
% diso = radial distance of observing point (km)
% 
% rg = range to satellite (km)
% azs = azimuth of the satellite observed from station (rads)
% els = elevation of the satellite observed from station (rads)

%Calculate ECI co-odinates of two points
[xit,yit,zit]=geo2eci(lont,latt,dist,timmat,day,month,year);
[xio,yio,zio]=geo2eci(lono,lato,diso,timmat,day,month,year);

we=7.29211510e-5;             %Earths rotation speed (rads/s)

%Calculate Greenwich Mean Sideral Time (GMST) 
[jd,dp2000]=jdate(day,month,year);
tu=dp2000/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 transformation co-ordinates
rx=xit-xio;
ry=yit-yio;
rz=zit-zio;

th=gmst+lono;
rs=sin(lato).*cos(th).*rx+sin(lato).*sin(th).*ry-cos(lato).*rz;
re=-sin(th).*rx+cos(th).*ry;
rz=cos(lato).*cos(th).*rx+cos(lato).*sin(th).*ry+sin(lato).*rz;

%Calculate required parameters
rg=sqrt(rs.^2+re.^2+rz.^2);
els=asin(rz./rg);
azs=atan(-re./rs);
azs(abs(azs)<1e-4)=0;
azs(rs>0)=azs(rs>0)+pi;
azs(azs<0)=azs(azs<0)+2*pi;
