%sealed container with higher pressure on side 2 and a permeable barrier to open AIR (same scenario as a balloon)
%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 = 106.7e5;           %moles of CO2 on right side for 16 psi, i.e.  +1 psi over ambient
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;
    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;
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")