%% satloc.m is design to estimate the location of a satellite based
%% on its observed azimuth and elevation from a ground station for
%% the case of a spherical Earth. The user inputs the geographic
%% longitude and latitude of the observing ground station, the azimuth
%% and elevation of the observed satellite and the known altitude of
%% the satellite. satloc.m then returns the geographic longitude and
%% latitude of the satellite. Derived with help from Orbital Motion by
%% A. E. Roy p31 521.3 R888o.

function [lons,lats]=satloc(long,latg,azs,els,alt)

% long = geographic longitude of the ground station (degrees)
% latg = geographic latitude of the ground station (degrees)
% azs = azimuth of the satellite observed from station (degrees)
% els = elevation of the satellite observed from station (degrees)
% alt = altitude of the satellite (km)
% lons = geographic longitude of the satellite (degrees)
% lats = geographic latitude of the satellite (degrees)

re=6371;                           %Earth radi (km)

%Convert angles to rads
long=deg2rad(long);               
latg=deg2rad(latg);               
azs=deg2rad(azs);
els=deg2rad(els);

%Calculate the great circle angle between station and satellite
gsth=-els+acos((re./(re+alt)).*cos(els));

%Calculate the longitude and latitude of the satellite
lond=asin(sin(azs).*sin(gsth));
lons=long+lond;
latd=acos(cos(gsth)./cos(lond));
poic=abs(azs)>(pi/2) & abs(azs)<((3/2)*pi);
latd(poic)=-latd(poic);
lats=latg+latd;

%Convert result angles to degrees
lons=rad2deg(lons);
lats=rad2deg(lats);
