%% zenith.m calculates the cos of the zenith angle at the known sub
%% ionospheric point. This is the value to multiple slant TEC to convert
%% it to vertical TEC.
%% 

function costh=zenith(long,latg,lons,lats,alts)

% long = ground station longitude (rads)
% latg = ground station latitude (rads)
% lons = sub ionospheric point longitude (rads)
% lats = sub ionospheric point latitude (rads)
% alts = sub ionospheric point altitude (km)
%
% costh = cos of the zenith angle

re=6371;                      %Earth radi (km)
sv=ones(size(long));

%Calculate seperation distance
cost=cos(lons-long).*cos(lats-latg);
dis=sqrt(re^2+(re+alts).^2-2*re*(re+alts).*cost);

%Calculate cos the the zenith angle
costh=(re^2-dis.^2-(re+alts).^2)./(-2*dis.*(re+alts));
