example3.m

Sky lab calculate position of satellite in GNSS systems and calculate elevation and azimuth them and

%EXAMPLE: Plot world-map with areas of visibility around major cities
%
% This script demonstrates some of the capabilities of Skylab, the
% Matlab version of SkyPlot. This demo reads an ephemeris-file and
% computes satellite-positions for a given epoch. For a certain
% satellite, the area where this satellite is visible on the world is
% shown on a world-map.

% ----------------------------------------------------------------------
% File.....: example3.m
% Date.....: 17-AUG-2000
% Version..: 1.0
% Author...: Peter Joosten
%            Mathematical Geodesy and Positioning
%            Delft University of Technology
% ----------------------------------------------------------------------

% ----------------------------------------------
% --- Initialize, close all previous figures ---
% ----------------------------------------------

close all;
clear all;

% -------------------------------------------------------------
% --- Read ephemeris (RINEX nav-file or YUMA almanac)       ---
% --- and delete double occurences of SVN's                 ---
% -------------------------------------------------------------

eph = sortrows(rdyuma('../data/yuma59.txt'));

% -----------------------------------------
% --- Station coordinates (WGS84 / XYZ) ---
% -----------------------------------------

station1(1) = 3923625; % Delft (Geodesy-building)
station1(2) =  298462;
station1(3) = 5002803;
cutoff     =      15; % Cutoff elevation in degrees

% -----------------
% --- Time/span ---
% -----------------

tsat = mktsat ('11-OCT-2000 12:00:00','11-OCT-2000 14:00:00',300);

% -------------------------------------------------------
% --- Compute satellite positions (XYZ WGS84 & AZ/EL) ---
% -------------------------------------------------------

[xsat,ysat,zsat,azim,elev,rotmat]    = cpaziele (tsat,eph,station1);

% --------------------------------------
% --- Create the figure, set options ---
% --------------------------------------

h = axesm ('MapProjection','Mercator', ...
	   'AngleUnits','degrees', ...
	   'MapLatLimit',[-80 80], ...
	   'MapLonLimit',[-180 180], ...
	   'grid','on', ...
	   'parallellabel','on', ...
	   'meridianlabel','on', ...
	   'mlabelparallel',-85, ...
	   'fontsize',7, ...
	   'frame','on');

set (gca,'position',[0 0 1 1]);

% ------------------------------------------
% --- Plot visibility area for station 1 ---
% ------------------------------------------

Rsat=(398600.44e9*(43200/(2*pi))^2)^(1/3);
Rsta=sqrt(sum(station1.^2));
station1=xyz2plh(station1);

viewangle1 = 180 - (90+cutoff) - rad2deg(asin((Rsta/Rsat)*sin(deg2rad(90+cutoff))));
[lat,lon]  = scircle1(rad2deg(station1(1)),rad2deg(station1(2)),viewangle1);

h = fillm (lat, lon, [0.9 0.9 0.9]);
set (h,'EraseMode','none');
h = plotm (lat,lon,'k-');

% --------------------------------
% --- Load and plot coastlines ---
% --------------------------------
  
load coast;
h = patchm (lat,long,[0.8 1.0 0.8]);

% ------------------------------
% --- Plot satellites tracks ---
% ------------------------------

lcol = get(gca,'ColorOrder');

for i = 1:size(eph,1);
  
  plh = xyz2plh ([xsat(i,:); ysat(i,:); zsat(i,:)]);
  idx = [0 find(abs(plh(2,1:size(plh,2)-1)-plh(2,2:size(plh,2)))>deg2rad(180)) size(plh,2)];
  
  for j = 1:length(idx)-1;

    h = plotm (rad2deg(plh(1,idx(j)+1:idx(j+1))),rad2deg(plh(2,idx(j)+1:idx(j+1))));
    set (h,'LineStyle','-');
    set (h,'LineWidth',2);
    set (h,'Color',lcol(mod(i,size(lcol,1))+1,:));
    
  end;
  
end;