clearvars
close all
axisfontsize=16;
panelletterfontsize=18;
figname='Fig1_s0.2_a0.5_Kmn50_Kmx250_nth30_nu0.2_d0.6';
figfolder= 'C:\Users\matlh\OneDrive - University of Leeds\My Projects\In-Leeds\Cooperator-Treshold system\Final Figure 1\';
% Load figure data
load([figfolder 'Fig1Data.mat'])
% Plot figure
figure; set(gcf,'Position',FigData.position);
hold on; plot(FigData.Nth(:,1),FigData.Nth(:,2),'b--')
hold on; plot(FigData.K(:,1),FigData.K(:,2),'k-')
hold on; plot(FigData.Nc(:,1),FigData.Nc(:,2),'b-')
hold on; plot(FigData.Nd_r(:,1),FigData.Nd_r(:,2),'r-')
hold on; plot(FigData.Nd_g(:,1),FigData.Nd_g(:,2),'g-')
xlim(FigData.xlim); ylim(FigData.ylim);
% Modify the plot
f=gcf;
set(f,'PaperPositionMode','auto');         
set(f,'PaperOrientation','landscape');
set(f,'Position',[332.2 61 964.8 704],'Name','Fig1_s0.2_a0.5_Kmn50_Kmx250_nth30_nu0.2_d0.6')
Axes = gca;
sub(2)=subplot(2,2,2,Axes);
xLim=xlim; xlim auto
% One single red line for simulated Nc. 
o=findobj(gca,'Type','Line',{'Color',[0,1,0],'-or','Color',[1,0,0]});
o(1).YData=nanmean(cat(1,reshape(o(1).YData,1,[]),reshape(o(2).YData,1,[])));
set(o(1),'Color',[1,0,0]); delete(o(2));
ND=cat(1,reshape(o(1).XData,1,[]),reshape(o(1).YData,1,[]));
o=findobj(gca,'Type','Line','Color',[0,0,1],'LineStyle','-');
NC=cat(1,reshape(o(1).XData,1,[]),reshape(o(1).YData,1,[]));
% Thicken deterministic lines (K).
set(findobj(gca,'Type','Line','Color',[0,0,0]),'LineWidth',2);
% Compute deterministic Nc and Nd.
S=0.2;A=0.5;Kmn=50;Kmx=250;Nth=30;
XLim=xlim; dt=XLim(2)/1e6; tpoints=XLim(1):dt:XLim(2); clearvars XLim
K=findobj(gca,'Type','Line','Color',[0,0,0]);
k(:,1)=K.XData; k(:,2)=K.YData; NoNans=~isnan(k(:,1))&~isnan(k(:,2));
K=interp1(k(NoNans,1),k(NoNans,2),tpoints);
Nc=nan(size(tpoints)); Nd=Nc;
Nc(1)=Nth; Nd(1)=K(1)-Nth;
for i=2:numel(tpoints)
    nc=Nc(i-1);
    nd=Nd(i-1);
    h=Nc(i-1)<Nth;
    k=K(i);
    n=Nc(i-1)+Nd(i-1);
    denom=(1-S)*nc/n+(1-A*h)*(nd/n);
    dnc=nc*((1-S)-(1-A*h)*n/k-(A*h-S)*nc/k)/denom;
    dnd=nd*((1-A*h)-(1-S)*n/k+(A*h-S)*nd/k)/denom;
    Nc(i)=Nc(i-1)+dnc*dt;
    Nd(i)=Nd(i-1)+dnd*dt;
end; clearvars k
% Plot simulated and deterministic N
o=findobj(gca,'Type','Line',{'Color',[0,0,1],'-or','Color',[1,0,0]});
hold on; plot(o(1).XData,o(1).YData+o(2).YData,'-','Color',[0 .75 0],'LineWidth',0.5);
plot(tpoints,Nd+Nc,'-','Color',[0 .75 0],'LineWidth',2);
% Plot deterministic Nc and Nd.
plot(tpoints,Nc,'b-','LineWidth',2); plot(tpoints,Nd,'r-','LineWidth',2);
% Set x-axis to zero.
xlim(xLim); ylim([0 300]);
x0=min(xlim); o=findobj(gca); 
arrayfun(@(x)set(x,'XData',x.XData-x0),o(2:end),'un',0);
o(1).XLim=o(1).XLim-min(o(1).XLim); 
% Aestetics
yticks([0,50,250])
xticks(0:10:max(xlim)); xticklabels({})
text(0.035,0.935,'K','HorizontalAlignment','left','VerticalAlignment','middle','Color','k','FontWeight','Bold','FontSize',axisfontsize,'Units','Normalized')
text(0.15,0.9,'N','Color',[0 .75 0],'HorizontalAlignment','left','VerticalAlignment','middle','FontWeight','Bold','FontSize',axisfontsize,'Units','Normalized')
text(0.25,0.7,'N_S','Color',[1 0 0],'HorizontalAlignment','left','VerticalAlignment','middle','FontWeight','Bold','FontSize',axisfontsize,'Units','Normalized')
text(0.35,0.35,'N_R','Color',[0 0 1],'HorizontalAlignment','left','VerticalAlignment','middle','FontWeight','Bold','FontSize',axisfontsize,'Units','Normalized')
text(0.2,0.075,'N_{th}','Color',[0 0 1],'LineStyle',':','HorizontalAlignment','left','VerticalAlignment','middle','FontSize',axisfontsize,'Units','Normalized')
set(gca,'Fontsize',axisfontsize)
text(gca,-.05,1.05,'(b)','Units','Normalized','HorizontalAlignment','right',...
        'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
text(gca,-1.35,1.05,'(a)','Units','Normalized','HorizontalAlignment','right',...
        'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
%%
sub(4)=subplot(2,2,4);
plot(tpoints,Nth./K,'k--');
hold on; plot(NC(1,:),NC(2,:)./(NC(2,:)+ND(2,:)),'b-');
plot(tpoints,Nc./(Nc+Nd),'b-','LineWidth',2);
ylim([0,0.8]);
% Set x-axis to zero.
yticks([0,Nth/Kmx,Nth/Kmn,1])
xlim(xLim);
x0=min(xlim); o=findobj(gca); 
arrayfun(@(x)set(x,'XData',x.XData-x0),o(2:end),'un',0);
o(1).XLim=o(1).XLim-min(o(1).XLim);
% Aesthetics
box off
xticks(0:10:max(xlim)); xlabel('Time (a.u.)');
text(0.3,0.8,'\bf $x\equiv\frac{N_R}{N}$','Color',[0 0 1],'HorizontalAlignment','center','VerticalAlignment','middle','FontSize',1.2*axisfontsize,'Units','Normalized','Interpreter','Latex')
text(0.2,0.25,'$\frac{N_{th}}{K_+}$','Color',[0 0 0],'HorizontalAlignment','center','VerticalAlignment','middle','FontSize',axisfontsize,'Units','Normalized','Interpreter','Latex')
text(0.9,0.85,'$\frac{N_{th}}{K_-}$','Color',[0 0 0],'HorizontalAlignment','center','VerticalAlignment','middle','FontSize',axisfontsize,'Units','Normalized','Interpreter','Latex')
set(gca,'Fontsize',axisfontsize)
text(gca,-.05,1.05,'(c)','Units','Normalized','HorizontalAlignment','right',...
        'VerticalAlignment','bottom','Color',[0,0,0],'FontWeight','Bold','FontSize',panelletterfontsize)
%% Save Figure
savefig(f,[figfolder figname '.fig'])
saveas(f,[figfolder figname '.png'],'png')
saveas(f,[figfolder figname '.pdf'],'pdf')