Plot2DTrajectory(GND1,EOC1,GND2,EOC2,H,logsout,RepT)

function Plot2DTrajectory(GND1,EOC1,GND2,EOC2,H,logsout,RepT)
% Script used within Repointing model
%
% Author: Ugo Pattacini

% plot 2D globe
figure;
landareas=shaperead('landareas.shp','UseGeoCoords',true);
axesm('flatplrq','Frame','on','Grid','on');
geoshow(landareas,'FaceColor',[1 1 .5],'EdgeColor',[.6 .6 .6]);

% plot Ground Stations locations
GND1=GND1*180/pi;
GND2=GND2*180/pi;
h(1)=plotm(GND1(1),GND1(2),'g^','LineWidth',1.5);
M{1}='GND1';
h(2)=plotm(GND2(1),GND2(2),'b^','LineWidth',1.5);
M{2}='GND2';
n=3;

% plot trajectory
h(n)=plotm(logsout.mu_l.Data(:,1),logsout.mu_l.Data(:,2),'k');
M{n}='Out of Coverage';
n=n+1;
    
% plot trajectory under GND1 visibility
c1=find(logsout.coverage.gnd.Data(:,1)>0);
if ~isempty(c1)
    tc1_1=logsout.coverage.gnd.Time(c1(1));
    tc1_2=logsout.coverage.gnd.Time(c1(end));
    i=find(logsout.mu_l.Time>tc1_1 & logsout.mu_l.Time<tc1_2);
    h(n)=plotm(logsout.mu_l.Data(i,1),logsout.mu_l.Data(i,2),'g');
    M{n}='GND1 Coverage';
    n=n+1;
    
    % plot GND1 visibility circle
    [vis_lat,vis_long]=calcVisibility(GND1,EOC1,H);
    plotm(vis_lat,vis_long,'g--');
end

% plot trajectory under GND2 visibility
c2=find(logsout.coverage.gnd.Data(:,2)>0);
if ~isempty(c2)
    tc2_1=logsout.coverage.gnd.Time(c2(1));
    tc2_2=logsout.coverage.gnd.Time(c2(end));
    j=find(logsout.mu_l.Time>tc2_1 & logsout.mu_l.Time<tc2_2);
    h(n)=plotm(logsout.mu_l.Data(j,1),logsout.mu_l.Data(j,2),'b');    
    M{n}='GND2 Coverage';    
    n=n+1;

    % plot GND2 visibility circle
    [vis_lat,vis_long]=calcVisibility(GND2,EOC2,H);
    plotm(vis_lat,vis_long,'b--');    
end

% plot trajectory under GND1 and GND2 visibility overlap
if exist('i','var') && exist('j','var')
    k=intersect(i,j);
    if ~isempty(k)
        h(n)=plotm(logsout.mu_l.Data(k,1),logsout.mu_l.Data(k,2),'c');
        M{n}='Coverage Overlap';    
        n=n+1;
    end
end

% plot repointing time
Data=reshape(logsout.coverage.repointing.Data,size(logsout.coverage.repointing.Time));
rep=find(Data(:,1)>0);
if ~isempty(rep)
    trep_1=logsout.coverage.repointing.Time(rep(1));
    trep_2=logsout.coverage.repointing.Time(rep(end));
    l=find(logsout.mu_l.Time>trep_1 & logsout.mu_l.Time<trep_2);
    t=logsout.mu_l.Time(l);
    lat=logsout.mu_l.Data(l,1);
    long=logsout.mu_l.Data(l,2);
    ll=[find(t>tc1_1 & t<tc1_2) find(t>tc2_1 & t<tc2_2)];
    if ~isempty(ll)
        lat=lat(ll);
        long=long(ll);
        h(n)=plotm(lat,long,'r','LineWidth',1.5);
        M{n}='Repointing Zone';
        lat_m=(min(lat)+max(lat))/2;
        long_m=max(long);
        textm(lat_m,long_m+10*pi/180,sprintf('\\leftarrow Repointing time = %d s',round(RepT)),...
              'FontSize',13,'Color','r');
    end
end

title('Satellite Trajectory');
legend(h,M);



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

global ellipsoid;
global lat;
global long;
global elev;

ellipsoid=almanac('earth','ellipsoid','meters');
lat=GND(1)*pi/180;
long=GND(2)*pi/180;
elev=EOC*pi/180;

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

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

vis_lat=vis_lat*180/pi;
vis_long=vis_long*180/pi;



function varargout=calcAltitude(flag,theta,r)

global ellipsoid;
global lat;
global long;
global elev;

[xl,yl,zl]=sph2cart(theta,elev,r);
[x,y,z]=lv2ecef(xl,yl,zl,lat,long,0,ellipsoid);
[phi,lambda,h]=ecef2geodetic(x,y,z,ellipsoid);

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