function plot_figures()

unzip('fig_data.zip')

set(0,'defaultTextInterpreter','latex')

plot_fig3('fig3a.mat','a')
plot_fig3('fig3b.mat','b')

plot_figs4to6('fig4.mat',10,4)
plot_figs4to6('fig5.mat',30,5)
plot_figs4to6('fig6.mat',10,6)

plot_fig7('fig7a.mat','a')
plot_fig7('fig7b.mat','b')

plot_fig8('fig8a.mat','a')
plot_fig8('fig8b.mat','b')

plot_figE1('figE1a.mat','a')
plot_figE1('figE1b.mat','b')

plot_figE2a('figE2a.mat')
plot_figE2b('figE2b.mat')

plot_figsE3toE4('figE3a.mat',3,'a')
plot_figsE3toE4('figE3b.mat',3,'b')

plot_figsE3toE4('figE4a.mat',4,'a')
plot_figsE3toE4('figE4b.mat',4,'b')

function plot_fig3(data,subfig)

load(data)

figure; surf(Plout,Qin,dl*mean_Rc_final,'EdgeColor','none'); view([0 90]); shading interp; colormap jet
ylabel({'$Q_{l,in}$','(m$^3$s$^{-1}$)'}); xlabel('$P_{l,out}$ (Pa)'); ylim([Qin(1) Qin(end)]); xlim([Plout(1) Plout(end)])
set(0,'defaultTextInterpreter','none')
cb=colorbar; title(cb,'$\overline{R}_c$ ($\mu$m)','Interpreter','latex')
set(0,'defaultTextInterpreter','latex')
saveas(gcf,['fig3' subfig '.fig'],'fig');

function plot_figs4to6(data,Ntime_pts,fig)

load(data)

figure; set(gcf,'Units','Normalized','OuterPosition',[0 0 0.5 1]);
for i=1:Ntime_pts
% Plot ox concn
subplot(2,2,1); plot(L*z,Cin*cO(:,i),'k-'); hold on 
xlabel('$z$ (m)'); ylabel({'$c_c^{\mathrm{O}}|_{r=R_c}$','(mol m$^{-3}$)'},'rot',0);

% Plot lactate concn
subplot(2,2,2); plot(L*z,CL*cL(:,i),'k-'); hold on
xlabel('$z$ (m)'); ylabel({'$c_c^{\mathrm{L}}|_{r=R_c}$','(mol m$^{-3}$)'},'rot',0);

% Plot fluid shear stress
subplot(2,2,3); plot(L*z,Pc*sigma_p(:,i),'k-',L*z,Pc*sigma_d*ones(Nz,1),'k--'); hold on
xlabel('$z$ (m)'); ylabel('$\sigma$ (Pa)','rot',0)

% Plot cell layer outer surface
subplot(2,2,4); 
if i==1
    plot(L*z,1e6*dl*Rc(:,i),'k--')
else
    plot(L*z,1e6*dl*Rc(:,i),'k-')
end
xlim([0 L]);
xlabel('$z$ (m)'); ylabel({'$R_c(z,t)$','($\mu$m)'},'rot',0); hold on
end
saveas(gcf,['fig' num2str(fig) '.fig'],'fig')

function plot_fig7(data,subfig)

load(data)

figure; [ax,p1,p2]=plotyy(Qin1,1e6*dl*mean_Rc_final1,Qin1,1e6*dl*std_dev_Rc_final_smthd);
set(ax,{'ycolor'},{'b';'r'})
xlabel('$Q_{l,in}$ (m$^3$s$^{-1}$)'); 
ylabel(ax(1),{'$\overline{R}_c$','($\mu$m)'},'rot',0);
ylabel(ax(2),{'SD','($\mu$m)'},'rot',0);
set(p1,'LineWidth',2);
set(p2,'LineWidth',2,'LineStyle','--','Color','r');
xlim(ax(1),[Qin1(1) Qin1(end)]); ylim(ax(1),[1e6*dl*mean_Rc_final1(1)-0.01 1e6*dl*mean_Rc_final1(end)+0.01])
xlim(ax(2),[Qin1(1) Qin1(end)]); ylim(ax(2),[1e6*dl*std_dev_Rc_final_smthd(end)-0.005 1.501]);
saveas(gcf,['fig7' subfig '.fig'],'fig')

function plot_fig8(data,subfig)

load(data)

figure; 
if strcmp(subfig,'a')
    i=find(abs(dl*mean_Rc_final1-(dl+dm+dc))==min(abs(dl*mean_Rc_final1-(dl+dm+dc))));
    plot(CLmax1,1e6*(dl*mean_Rc_final1-(dl+dm+dc)),'k',[CLmax1(i) CLmax1(i)],[-50 20],'k','LineWidth',0.5)
elseif strcmp(subfig,'b')
    plot(CLmax1,1e6*(dl*mean_Rc_final1-(dl+dm+dc)),'k')
end
xlabel('$C_{max}^{\mathrm{L}}$ (mol m$^{-3}$)'); ylabel({'$\overline{\Delta R_c}$','($\mu$m)'},'rot',0)
saveas(gcf,['fig8' subfig '.fig'],'fig')

function plot_figE1(data,subfig)

load(data)

% Plot fluid pressure
figure; surf(L*z,dl*1e6*[rl;rm;re],Patm+[Pl*pl_mat;Pm*pm;Pc*p_e],'EdgeColor','None'); zlim([1e5 1.43e5]); view([128 28])
xlabel('$z$ (m)'); ylabel('$r$ ($\mu$m)'); zlabel('$p$ (Pa)','rot',0); hold on
plot3(L*z,dl*1e6*ones(Nz,1),Patm+Pl*pl,'k--',L*z,dl*1e6*Rm*ones(Nz,1),Patm+Pm*pm(Nr,:),'k--',L*z,dl*1e6*Rc_vec,Patm+Pc*p_e(round((Rc_vec(1)-Rm)/(Re-Rm)*2*Nr),:),'w-.','Linewidth',2);
saveas(gcf,['figE1' subfig '.fig'],'fig')

function plot_figE2a(data)

load(data)

% Plot lumen axial fluid velocity
figure; surf(L*z,dl*1e6*rl,U*ulz,'EdgeColor','None'); zlim([0 0.5375]); view([128 28]) %view([0 90])
xlabel('$z$ (m)'); ylabel('$r$ ($\mu$m)'); zlabel('$u_{l,z}$ (ms$^{-1}$)','rot',0)
saveas(gcf,'figE2a.fig','fig')

function plot_figE2b(data)

load(data)

% Plot radial fluid velocity
figure; surf(L*z,dl*1e6*[rl;rm;re],epsilon*U*[ulr;umr;uer],'EdgeColor','None'); zlim([0 0.8e-4]); view([128 28]) %view([0 90])
xlabel('$z$ (m)'); ylabel('$r$ ($\mu$m)'); zlabel('$u_r$ (ms$^{-1}$)','rot',0); hold on
plot3(L*z,dl*1e6*ones(Nz,1),epsilon*U*ulr(Nr,:),'k--',L*z,dl*1e6*Rm*ones(Nz,1),epsilon*U*umr(Nr,:),'k--',L*z,dl*1e6*Rc_vec,epsilon*U*lambda.^2*1./Rc_vec'.*pl/16,'w-.','Linewidth',2);
saveas(gcf,'figE2b.fig','fig')

function plot_figsE3toE4(data,fig,subfig)

load(data)

% Plot concentration throughout bioreactor
figure; surf(L*z,dl*1e6*[rl;rm;re],C*[cl;cm;ce],'EdgeColor','None'); shading interp; 
if fig==3
    view([128 40])
elseif fig==4
    view([42 14])
end
hold on
xlabel('$z$ (m)'); ylabel('$r$ ($\mu$m)'); zlabel({str,'(mol m$^{-3}$)'},'rot',0)
plot3(L*z,dl*1e6*ones(Nz,1),C*cl(Nr,:),'k--',L*z,dl*1e6*Rm*ones(Nz,1),C*cm(Nr,:),'k--',L*z,dl*1e6*Rc,C*(Ac+Bc.*Rc.^nc+m*M*Rc.^2./(2*(nc-2))),'w-.','LineWidth',2); hold off
saveas(gcf,['figE' num2str(fig) subfig '.fig'],'fig')