% Porkchop plotter
%
% Requires data file named 'con'
% Format of con:
%
% JDi(1)  JDf(1) Vinf_earth Vinf_mars
% JDi(1)  JDf(2) Vinf_earth Vinf_mars
% JDi(1)  JDf(3) Vinf_earth Vinf_mars
% JDi(1)  JDf(4) Vinf_earth Vinf_mars
% JDi(2)  JDf(1) Vinf_earth Vinf_mars
% JDi(2)  JDf(2) Vinf_earth Vinf_mars
% JDi(2)  JDf(3) Vinf_earth Vinf_mars
% JDi(2)  JDf(4) Vinf_earth Vinf_mars
%
% (Vinf in m/s)
%
% Plots C3 at Earth, and Vinf at Mars
%
clear
clf
% Load data file
load con;

% Sort data into matrices required for contour plot
remem = con(1,1);
colu = 1;
rowu = 0;
for (i = 1:length(con(:,1)))
  if (con(i,1) ~= remem)
    remem = con(i,1);
    colu = colu + 1;
    rowu = 1;
    departmat(colu,rowu) = con(i,1);
    arrivemat(colu,rowu) = con(i,2);
    depvinf(colu,rowu) = con(i,3);
    arrvinf(colu,rowu) = con(i,4);
  else
    rowu = rowu + 1;
    departmat(colu,rowu) = con(i,1);
    arrivemat(colu,rowu) = con(i,2);
    depvinf(colu,rowu) = con(i,3);
    arrvinf(colu,rowu) = con(i,4);
  end 
end

%    Convert to Modified Julian Date (MJD)
% (MJD = 0.0 on 1858 November 17 at 0 hrs UTC)
departmat = departmat - 2400000.5;
arrivemat = arrivemat - 2400000.5;

% Plot C3L at Earth and Vinf at Mars
[C2,h2] = contour(departmat,arrivemat,(depvinf/1000).^2, [17.5,20,22.5,25,27.5],'r');
clabel(C2,h2);
hold on
[C3,h3] = contour(departmat,arrivemat,(arrvinf/1000), [2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0 5.0], 'b--');
clabel(C3,h3);

% Obtain tick mark labels
departmin = min(min(departmat));
departmax = max(max(departmat));
xtickmat = departmin:(departmax-departmin)/5:departmax;
for i = 1:length(xtickmat)
  xticklabel(i,:) = num2str(xtickmat(i),'%6.1f');
end
set(gca,'XTick',xtickmat)
set(gca,'XTickLabel',{xticklabel})

arrivemin = min(min(arrivemat));
arrivemax = max(max(arrivemat));
ytickmat = arrivemin:(arrivemax-arrivemin)/5:arrivemax;
for i = 1:length(ytickmat)
  yticklabel(i,:) = num2str(ytickmat(i),'%6.1f');
end
set(gca,'YTick',ytickmat)
set(gca,'YTickLabel',{yticklabel})

grid
hold off
xlabel('Earth Departure Date (MJD)')
ylabel('Mars Arrival Date (MJD)')
legend([h2(1),h3(1)],'C3L (km^2/s^2)','Arrival V_\infty (km/s)',2)

print -depsc2 e2m.eps
