N=NumOfContacts

function N=NumOfContacts
% Compute over the whole repeat-cycle a statistics on the number of
% contacts, given the GS scenario (defined hereafter)
%
% Author: Ugo Pattacini

NOrbit=175;                     % orbits of repeat cycle
TOrbit=5913.8606;               % orbital period [s]
H=692844.666;                   % satellite altitude [m]
EOC=5;                          % GND end of coverage []
Ts=5;                           % orbit trajectory sample time [s]

GND{1}=[78.22 15.38];           % Svalbard
GND{2}=[67.88 20.25];           % Kiruna
GND{3}=[40.3852 16.4210];       % Matera
GND{4}=[41.978 13.604];         % Fucino
GND{5}=[53.2126 -105.9353];     % Prince Albert
GND{6}=[-31.55241 -64.4636];    % Cordoba
GND{7}=[-72 2.53];              % Troll

%%%%% DEFINE HERE THE GS SCENARIO
% create visibility circles
[vis_lat,vis_lon]=calcVisibility(GND{1},EOC,H); % Svalbard
GNDData.Lat{1}=vis_lat;GNDData.Lon{1}=vis_lon;
[vis_lat,vis_lon]=calcVisibility(GND{3},EOC,H); % Matera
GNDData.Lat{2}=vis_lat;GNDData.Lon{2}=vis_lon;
[vis_lat,vis_lon]=calcVisibility(GND{5},EOC,H); % Prince Albert
GNDData.Lat{3}=vis_lat;GNDData.Lon{3}=vis_lon;

OrbitalData=load('OrbitalData.mat');

ellipsoid=almanac('earth','ellipsoid','meters');

axesm mercator;

% plot lands
for i=1:length(OrbitalData.Lat)
    plotm(OrbitalData.Lat{i},OrbitalData.Lon{i},'k');
end

% plot visibility circles
for i=1:length(GNDData.Lat)
    plotm(GNDData.Lat{i},GNDData.Lon{i},'r');
end

SC_N=round(TOrbit/Ts);
M=length(OrbitalData.StateVector);
L=min(NOrbit,floor(M/SC_N));

% init
N=zeros(L,4);

% main loop
for i=1:L
    if exist('hdl1','var')
        delete([hdl1 hdl2]);
    end
    
    k1=1+SC_N*(i-1);
    k2=min(SC_N*i,M);
    
    % get orbital trajectory
    x=OrbitalData.StateVector(k1:k2,1);
    y=OrbitalData.StateVector(k1:k2,2);
    z=OrbitalData.StateVector(k1:k2,3);

    % transformation
    [lat,lon]=ecef2geodetic(x,y,z,ellipsoid);
    lat=lat*180/pi;
    lon=lon*180/pi;

    idx=isInside(GNDData,lat,lon);
    N(i,1)=Ts*length(idx{1});
    N(i,2)=Ts*length(idx{2});
    N(i,3)=Ts*length(idx{3});
    N(i,4)=Ts*length(idx{4});
        
    title(sprintf('Orbit #%d',i));

    % plot
    hdl1=plotm(lat,lon,'b.');
    hdl2=plotm(lat(idx{4}),lon(idx{4}),'r.');
    drawnow;  
end



function [vis_lat,vis_lon]=calcVisibility(GND,EOC,H)

global g_ellipsoid;
global g_lat;
global g_lon;
global g_elev;

g_ellipsoid=almanac('earth','ellipsoid','meters');
g_lat=GND(1)*pi/180;
g_lon=GND(2)*pi/180;
g_elev=EOC*pi/180;

theta=(0:1:360)*pi/180;
L=length(theta);
vis_lat=zeros(L,1);
vis_lon=zeros(L,1);

for i=1:L
    f=@(r)(calcAltitude(false,theta(i),r)-H);
    r=fzero(f,1e5);
    [vis_lat(i) vis_lon(i)]=calcAltitude(true,theta(i),r);
end

vis_lat=vis_lat*180/pi;
vis_lon=vis_lon*180/pi;

% refining needed for near-polar ground stations
i=find(180-vis_lon<5,1,'last');
j=find(180+vis_lon<5,1,'first');
if ~isempty(i) && ~isempty(j)
    if i<length(vis_lon)
        k=i+1;
    else
        k=1;
    end
    s=sign(vis_lat(i));
    tmp_lon=[vis_lon(1:i);180;180;-180;-180];
    tmp_lat=[vis_lat(1:i);vis_lat(i);s*90;s*90;vis_lat(k)];
    if i<length(vis_lon)
        vis_lon=[tmp_lon;vis_lon(i+1:end)];
        vis_lat=[tmp_lat;vis_lat(i+1:end)];
    else
        vis_lon=tmp_lon;
        vis_lat=tmp_lat;
    end
end



function varargout=calcAltitude(Type,theta,r)

global g_ellipsoid;
global g_lat;
global g_lon;
global g_elev;

[xl,yl,zl]=sph2cart(theta,g_elev,r);
[x,y,z]=lv2ecef(xl,yl,zl,g_lat,g_lon,0,g_ellipsoid);
[phi,lambda,h]=ecef2geodetic(x,y,z,g_ellipsoid);

d=size(r);
if Type
    varargout{1}=reshape(phi,d);
    varargout{2}=reshape(lambda,d);
    varargout{3}=reshape(h,d);    
else
    varargout{1}=reshape(h,d);    
end



function idx=isInside(VectData,lat,lon)

L=length(VectData.Lat);
idx=cell(L+1,1);

for i=1:L
    y=inpolygon(lat,lon,VectData.Lat{i},VectData.Lon{i});
    idx{i}=find(y);
    
    idx{L+1}=union(idx{L+1},idx{i});    
end