%function skyplot(orbfile,sv,latR,lonR,altR,cutoff);

% SKYPLOT	Produces satellite skyplot from sp3 file.
%           Input:
%             orbfile = an orbit file in sp3 format
%             sv = vector of PRN numbers
%             latR, lonR, altR = ground site position (WGS)
%             cutoff = elevation cutoff angle in degrees
%
%           Call: skyplot(orbfile, sv, latR, lonR, altR, cutoff)

polar(0,90);
hold on;

% convert reference site to ECEF
[XR,YR,ZR] = wgs2xyz(lonR,latR,altR);

disp(['-------------']);
disp(['Producing sky plot...']);

type = 'nav';

for (i=1:length(sv))

  if (type == 'sp3')

    % read sp3 file for satellite sv and convert to meters
    [Xs,Ys,Zs,dT,Ts] = read_sp3(orbfile,sv(i));

    % matrix of satellite positions (ECEF, meters)
    S = [Xs.*1000 Ys.*1000 Zs.*1000];

  elseif (type == 'nav')

    % figure out day of week and calculate time offset in seconds
    [week,dow,mjd] = gweek(2015,278);
    DT = dow*24*3600;

    % read rinex navigation file
    eph = read_rinexn(orbfile);

    % time vector of observations in seconds
    Ts = 0:30:24*60*2;
    S = [];
    for (j=1:length(Ts))
      satpos = get_satpos(Ts(j)+DT,sv(i),eph,0);
      if (~isempty(satpos)) S_tmp = [satpos(1) satpos(2) satpos(3)]; end
      S = [S_tmp ; S];
    end
  end

  % compute elevation angle, range, and azimuth
  [azim,elev,hlen] = azelle(S,[XR,YR,ZR]);

  % convert elevation to zenith angle
  zen_ang = (pi/2) - elev;

  % cutoff observations below cutoff angle
  cutoff = cutoff * pi/180;
  I = find(zen_ang<((pi/2)-cutoff));

  % sky plot
  polar(azim(I),zen_ang(I).*180/pi,'*b');
  %h = mmpolar(azim(I),zen_ang(I).*180/pi,'*b','TZeroDirection','North','RLimit',[0 90]);

end

