%% orbcal2.m calculates the geographical longitude and
%% latitude points of interest for a satellite in a single
%% orbit. It works by first determining the points of interest
%% in the orbital plane, and then transforms these points into
%% geographical co-ordinates using orbtrsf.m.
%% orbcal2.m is different from orbcal.m in that it calculates the
%% orbit centred not at the equator (given by timut in UT), but 
%% at some specific reference time (in UT).
%%

function [lon,lat]=orbcal2(incl,elon,timut,timutr,orbp,timv,diro)

% incl = inclination of the orbit (degrees)
% elon = longitude at which orbit crosses equator (degrees)
% timut = univerisal time (UT) when orbit crosses equator
% timutr = univerisal time (UT) when orbit crosses geographic reference point
% orbp = orbital period (fractional hours)
% timv = UT orbital time vector (fractional hours)
% diro = direction of the orbit
% lon = geographical longitude of orbit (rads)
% lat = geographical latitude of orbit (rads)

th=diro*((timv-timut)/orbp)*(2*pi);      %required theta values in orbital plane
xcop=cos(th);                            %x co-odinates in orbital plane
ycop=sin(th);                            %y co-odinates in orbital plane
zcop=zeros(size(xcop));                  %z co-odinates in orbital plane

% Use orbtrsf to transform to the geographical co-ordinate system
[xcg,ycg,zcg]=orbtrsf(xcop,ycop,zcop,incl,elon,2);

% Convert from cartesian to spherical co-ordinates
[lon,lat,radg]=cart2sph(xcg,ycg,zcg);

we=2*pi/24;                              %angular velocity of Earths rotation (rads per hour)
lon=lon-we*(timv-timutr);                %correction to compenstate for rotating Earth
lon=longch(lon);                         %correct longitude values so all values are between -pi and pi

