%CO2 regulated container with permeable barrier to open AIR
%2302026
%segallis


clear
perm = 1/2e5;           %permeability per unit time
or = 1.17*perm;         %O2 rate of diffusion
cr = 1*perm;            %CO2 rate of diffusion
nr = 1.25*perm;         %N2 rate of diffusion
o1 = 21e5;              %moles of O2 on left side
c1 = 0;                 %moles of CO2 on left side
n1 = 79e5;              %moles of N2 on left side
o2 = 0;                 %moles of O2 on right side
c2 = 200e5;             %moles of CO2 on right side
%c2 = 100e5;             %moles of CO2 on right side for 15 psi, i.e.  equal to ambient
t2 = c2;                %total moles on right side
n2 = 0;
k = 15/(o1+c1+n1);      %RT/V so we get 1 ATM (15 psi) on left at start of simulation

for t = 1:1e6;
    if mod(t,100e3) == 0
        fprintf("\r%d  ", t)
    endif
    olr = o1*or - o2*or;
    clr = c1*cr - c2*cr;
    nlr = n1*nr - n2*nr;
    o2 = o2 + olr;
    n2 = n2 + nlr;
    c2 = c2 + clr;
                               %cnew = -(clr + olr + nlr);                  %new CO2 provided by tank
    cnew = t2 - (c2+o2+n2);                     %new CO2 provided by tank
    c2 = c2 + cnew;                             %maintains exact pressure
    %c2 = c2 + max(0, cnew);                     %regulator allows pressure to rise above setpoint
    p1(t) = k*o1 + k*c1 + k*n1;
    p2(t) = k*o2 + k*c2 + k*n2;
    oo1(t) = o1;
    oo2(t) = o2;
    cc1(t) = c1;
    cc2(t) = c2;
    nn1(t) = n1;
    nn2(t) = n2;
    lr(t) = olr + clr + nlr;
    oolr(t) = olr;
    cclr(t) = clr;
    nnlr(t) = nlr;
    op(t) = o2/(o2+c2+n2)*100;
    cp(t) = c2/(o2+c2+n2)*100;
    np(t) = n2/(o2+c2+n2)*100;
    ccnew(t) = cnew;
end
fprintf("\n")

figure(1)
hold off
plot(p1)
hold on
plot(p2)
title("Pressure")

figure(2)
hold off
plot(oo1)
hold on
plot(oo2)
title("Oxygen")

figure(3)
hold off
plot(cc1)
hold on
plot(cc2)
title("CO2")

figure(32)
hold off
plot(nn1)
hold on
plot(nn2)
title("N2")

figure(4)
plot(lr)
title("Net Moles Left to Right")

figure(5)
hold off
plot(oolr)
hold on
plot(cclr)
plot(nnlr)
title("O2, CO2, N2 Moles Left to Right")
legend("O2", "CO2", "N2")

figure(41)
hold off
plot(op)
hold on
plot(cp)
hold on
plot(np)
title("Percent")
legend("O2", "CO2", "N2")

figure(6)
plot(ccnew)
title("Ctank")