clear; %clc; close all;

nx = 20;
L  = 0.238;  % length of domain, m


Dam = 8.48 / 1000^2;  % binary diffusivity, m^2/s
DaA = 13.72/ 1000^2;  % binary diffusivity, m^2/s
DmA = 19.91/ 1000^2;  % binary diffusivity, m^2/s

% ordering: acetone, methanol, air
Dbin = [ 0   Dam  DaA;
         Dam 0    DmA;
         DaA DmA  0 ];

x0 = [ 0.319; 0.528 ];  % mole fractions at bottom
xL = [ 0; 0; ];

ns = 3;

ct = 101325 / (8.314 * 328.5);

%% Step 1: guess Ni
Na = 1.3e-3;
Nm = 2.15e-3;
N = [ Na, Nm, 0 ];

%% Step 2: Calculate Phi, phi
Phi = zeros(ns-1);
phi = zeros(ns-1,1);
for i=1:ns-1
  for j=1:ns-1
    if i==j
      sumterm = 0;
      for k=1:ns
        if k~=i
          sumterm = sumterm + N(k)/(ct*Dbin(i,k)/L);
        end
      end
      Phi(i,i) = N(i) / (ct*Dbin(i,ns)/L) + sumterm;
    else
      Phi(i,j) = N(i) * ( 1/(ct*Dbin(i,ns)/L)-1/(ct*Dbin(i,j)/L) );
    end
  end
  phi(i) = - N(i)/( ct*Dbin(i,ns)/L );
end

%% Step 3: calculate compositions at the "top" of the tube
eta = linspace(0,1,nx);
xtop = zeros(ns,nx);
for i=1:nx
  xtop(1:ns-1,i) = expm( Phi*eta(i) ) * x0 + ( expm(Phi*eta(i)) - eye(ns-1) ) * Phi^-1 * phi;
  xtop(ns,i) = 1-sum(xtop(1:ns-1,i));
end

xerr = norm( xL - xtop(1:ns-1,nx) )

plot(eta,xtop);