%% Analyze GPS Satellite Visibility
% This example shows how to simulate and analyze GPS satellite visibility at specified receiver positions and times. Use live script controls to set different parameters for the satellite simulation.
%% Specify Simulation Parameters
% Specify the start time, duration in hours, and time between samples in seconds of the simulation. Also specify the position of the receiver in geodetic coordinates and the mask angle, or minimum elevation angle, of the receiver.

dateString = "2021-04-08 00:00:00.000";
startTime = datetime(dateString,InputFormat="yyyy-MM-dd HH:mm:ss.SSS");
numHours = 24;
dt =60; % s
latitude =42.3013162; % deg
longitude = -71.3782972; % deg
altitude = 50; % m
recPos = [latitude longitude altitude]; % [deg deg m]
maskAngle = 5; % deg

%
%% Generate Satellite Visibilities
% Using the parameters, generate the satellite visibilities as a matrix of logical values. Each row in the matrix corresponds to a time step, and each column corresponds to a satellite. To plot visibilities, iterate through the time vector calculating the satellite positions and look angles based on GNSS constellation simulation.

secondsPerHour = 3600;
timeElapsed = 0:dt:(secondsPerHour*numHours);
t = startTime + seconds(timeElapsed);

satPos = gnssconstellation(t(1));
numSats = size(satPos,1);
numSamples = numel(t);
az = zeros(numSamples,numSats);
el = zeros(numSamples,numSats);
vis = false(numSamples,numSats);
satIDs = 1:numSats;

sp = skyplot([],[]);

for ii = 1:numel(t)
    satPos = gnssconstellation(t(ii));
    [az(ii,:),el(ii,:),vis(ii,:)] = lookangles(recPos,satPos,maskAngle);
    set(sp,AzimuthData=az(ii,vis(ii,:)), ...
        ElevationData=el(ii,vis(ii,:)), ...
        LabelData=satIDs(vis(ii,:)));
    drawnow limitrate
end

%
%% Plot Results
% Use the logical matrix to generate a satellite visibility chart and plot the total number of visible satellites at each time step. In general, at least four satellites must be visible to compute a positioning solution.

visPlotData = double(vis);
visPlotData(visPlotData == false) = NaN; % Hide invisible satellites.
visPlotData = visPlotData + (0:numSats-1); % Add space to satellites to be stacked.
colors = colororder;

figure
plot(t,visPlotData," .",Color=colors(1,:))
yticks(1:numSats)
yticklabels(string(satIDs))
grid on
ylabel("Satellite ID")
xlabel("Time")
title("Satellite Visibility Chart")
axis tight

%

numVis = sum(vis,2);
figure
area(t,numVis)
grid on
xlabel("Time")
ylabel("Number of satellites visible")
title("Number of GPS satellites visible")
axis tight

% Copyright 2021 The MathWorks, Inc.