% Script to plot Figure S5 in supplementary material

st_dir = '../'
file = '../Tables.xlsx';
T = readtable(file, "Sheet","sliprateDoubling_figS5");

sim = T.sim;
tstart = T.tstart_yr;
tend2 = T.tend_yr;
Seff0 = T.sigma_MPa;


mmyr2ms = 1000*365*24*3600;
vpl = 50/mmyr2ms; % plate rate in m/s
yr2sec = 365*24*3600;


for k = 1:length(sim)
data = load([st_dir 'Vmax_sim' num2str(sim(k))]);
t1 = data(:,1)/yr2sec; tend = t1(end);
vmax = data(:,2);
vmax = vmax/vpl;

data2 = load([st_dir 'porePressure_sim' num2str(sim(k))]);
Seff = data2(:,5)/1e6;
pressure = (Seff0(k) - Seff)/Seff0(k); % in m/s
t2=data2(:,1)/yr2sec; % in yrs


figure;
set(gcf,'color','w')
set(gca,'FontSize',14)
semilogy(t1,vmax,'b-');
xlabel('Time (yr)','fontName','Times');
ylabel('V_{max}/V_{pl}');
ax = gca; ax.FontSize = 14;
%ylim([1 10^4]);

yyaxis right
plot(t2,pressure,'-r');
ylim([0 1]);
ylabel('\Delta p_f at fault center (MPa)');
xlabel('time (years)');
ax = gca; ax.FontSize = 14;

xlim([tstart(k) tend2(k)]);
end