%Code to solve PDE model where oxygen and antimicrobial are in
%a quasi-steady state


t_max =48; 
H = 20;
delta_t = .025; 
delta_z = .1; 
vb = 10;
tot_space_steps = (H+vb)/delta_z;

ew = 0.7;
ex0 = (1-ew)*.9;
A0 = 0;
O0 = [zeros(tot_space_steps+1,1)];
Ab0 = 0;

ex_bound=0;
ez_bound=0;
O_bound=1;

D_o = 5e6; 
D_a = 5e3; %5e3-5e6
delta_ox = 2e5;
delta_xa = 13.54;
delta_xphi = 3.8;
delta_ax = 320;
beta_a = 0.8125;
t_bar = 8;
epsilon = 1; %0.2-10
Q =0.2; %0.1-1
delta_a=0.0; %0 or 0.01 or 0.1
delta_ab=0; %0 or 5e3
eb=0.0; %0 or 0.01



[ex,ez,A,O,Ab,phi] = solve_biofilm_system(t_max,H,delta_t,delta_z,ex0,A0,O0,Ab0,O_bound,D_o,D_a,eb,delta_ox,delta_xa,delta_xphi,delta_ax,beta_a,t_bar,epsilon,Q,vb,ew,delta_a,delta_ab);
 
time_gap = 4;
no_of_iterations = t_max/time_gap + 1;
hold on 
for i = 1:(no_of_iterations)
    solution = ex(:,1+time_gap*(i-1)/delta_t);
    set(gca,'FontSize',20)
    plot(0:delta_z:H+vb,solution,'LineWidth', 2)
    xlabel('z')
    ylabel({'live cell fraction: \epsilon_x(z,t)'})
    ylim([0,inf])
    xlim([0,inf])
end 
hold off