varargout=GenInputData(Type)

function varargout=GenInputData(Type)
% Generate input data for DESOS
% Possible data types are:
%
% 'contacts':     once defined the scenario hereafter in the code,
%                 output the vector containing contacts with selected
%                 ground stations
%
% 'acquisitions': output acquisition data rates according to the number of
%                 active instrument polarizations, working time period of
%                 the modes, and the geo-localization of the S/C (over
%                 land, overseas)
%
% 'both':         output both contacts and acquisitions vectors
%
% Note that S/C state vector is supposed to be stored in the file
% 'OrbitalData.mat'
%
% Example: [contact,acquisitions]=GenInputData('both'); 
%
% Author: Ugo Pattacini

cont=false;
acq=false;
switch lower(Type)
    case 'contacts'
        cont=true;
    case 'acquisitions'
        acq=true;
    case 'both'
        cont=true;
        acq=true;
    otherwise
        error('Unknown Type');
end

%%%%% ORBITAL PARAMETERS
NOrbit=175;                     % orbits of repeat cycle
TOrbit=5913.8606;               % orbital period [s]
H=692844.666;                   % satellite altitude [m]
EOC=5;                          % GND End Of Coverage EOC []
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 ACQ SCENARIO
% acquisitions time [s] and rate [bps] over land
AcqTbl_Land{1}.time=2.5*60; AcqTbl_Land{1}.rate=352e6; AcqTbl_Land{1}.numpol=1;
AcqTbl_Land{2}.time=2.5*60; AcqTbl_Land{2}.rate=352e6; AcqTbl_Land{2}.numpol=2;
AcqTbl_Land{3}.time=7.5*60; AcqTbl_Land{3}.rate=301e6; AcqTbl_Land{3}.numpol=1;
AcqTbl_Land{4}.time=7.5*60; AcqTbl_Land{4}.rate=301e6; AcqTbl_Land{4}.numpol=2;

% acquisitions time [s] and rate [bps] overseas
AcqTbl_Water{1}.time=40*60; AcqTbl_Water{1}.rate=6.7e6; AcqTbl_Water{1}.numpol=1;

AcqTblLen_Land=length(AcqTbl_Land);
AcqTblLen_Water=length(AcqTbl_Water);

% calculate total acquisition time per orbit [s]
AcqT_Land=0;
for j=1:AcqTblLen_Land
    AcqT_Land=AcqT_Land+AcqTbl_Land{j}.time;
end

AcqT_Water=0;
for j=1:AcqTblLen_Water
    AcqT_Water=AcqT_Water+AcqTbl_Water{j}.time;
end

%%%%% 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
t=0;
Contacts=[];
Acquisitions=[];

% main loop
for i=1:L
    title(sprintf('Orbit #%d',i));
    
    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;      
    
    orbitTime=(t:Ts:t+(k2-k1)*Ts)';

    % update Contacts
    if cont
        idx=isInside(GNDData,lat,lon);    
        tmpContacts=-ones(k2-k1+1,2);
        tmpContacts(:,1)=orbitTime;
        tmpContacts(idx,2)=1;
        
        % refine
        d=diff(tmpContacts(:,2));
        tmpContacts(find(~d)+1,:)=[];
        Contacts=cat(1,Contacts,tmpContacts);

        % plot
        if ~acq
            if exist('hdl1','var')
                delete([hdl1 hdl2]);
            end
            hdl1=plotm(lat,lon,'b.');
            hdl2=plotm(lat(idx),lon(idx),'r.');
        end
        drawnow;
    end   
    
    % update Acquisitions
    if acq
        [idx_land,idx_water]=isInside(OrbitalData,lat,lon);
        idx1=RefineAcqIdx(idx_land,idx_water,Ts,AcqT_Land);
        idx_land=setdiff(idx_land,idx1);
        idx_water=setdiff(idx_water,idx1);
        idx2=RefineAcqIdx(idx_water,idx_land,Ts,AcqT_Water);        
        tmpAcquisitions=-ones(k2-k1+1,2);
        tmpAcquisitions(:,1)=orbitTime;
        
        % fill with the data rates over land
        l=round(rand*(AcqTblLen_Land-1))+1;
        cumT=0;
        for j=1:length(idx1)
            tmpAcquisitions(idx1(j),2)=AcqTbl_Land{l}.rate;
            tmpAcquisitions(idx1(j),3)=AcqTbl_Land{l}.numpol;
            cumT=cumT+Ts;
            if cumT>=AcqTbl_Land{l}.time
                cumT=0;
                l=l+1;
                if l>AcqTblLen_Land
                    l=1;
                end
            end
        end
        
        % fill with the data rates overseas
        l=round(rand*(AcqTblLen_Water-1))+1;
        cumT=0;
        for j=1:length(idx2)
            tmpAcquisitions(idx2(j),2)=AcqTbl_Water{l}.rate;
            tmpAcquisitions(idx2(j),3)=AcqTbl_Water{l}.numpol;
            cumT=cumT+Ts;
            if cumT>=AcqTbl_Water{l}.time
                cumT=0;
                l=l+1;
                if l>AcqTblLen_Water
                    l=1;
                end
            end
        end

        % refine
        d=diff(tmpAcquisitions(:,2).*tmpAcquisitions(:,3));
        tmpAcquisitions(find(~d)+1,:)=[];
        Acquisitions=cat(1,Acquisitions,tmpAcquisitions);
        
        % plot
        if exist('hdl1','var')
            delete([hdl1 hdl2 hdl3]);
        end
        hdl1=plotm(lat,lon,'b.');
        hdl2=plotm(lat(idx1),lon(idx1),'k.');
        hdl3=plotm(lat(idx2),lon(idx2),'c.');
        drawnow;
    end
    
    t=SC_N*i*Ts;
end

% final refining
if cont
    d=diff(Contacts(:,2));
    Contacts(find(~d)+1,:)=[];
end
if acq
    d=diff(Acquisitions(:,2).*Acquisitions(:,3));
    Acquisitions(find(~d)+1,:)=[];
end

% set output
if cont && acq
    varargout{1}=Contacts;
    varargout{2}=Acquisitions;
elseif cont
    varargout{1}=Contacts;
else
    varargout{1}=Acquisitions;
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 varargout=isInside(VectData,lat,lon)

d=size(lat);
y=reshape(zeros(length(lat),1),d);

L=length(VectData.Lat);
for i=1:L
    yi=inpolygon(lat,lon,VectData.Lat{i},VectData.Lon{i});
    y=y+double(yi);
end

varargout{1}=find(y);
if nargout>1
    varargout{2}=find(~y);
end



function idx=RefineAcqIdx(idx1,idx2,Ts,AcqT)

l=length(idx1);
m=floor(AcqT/Ts);
n=m-l;

if n>0
    i=round(rand*(length(idx2)-n))+1;
    idx=union(idx1,idx2(i:i+n-1));
else
    i=round(-rand*n)+1;
    idx=idx1(i:i+m-1);
end