%% The reactive distillation model
parameters_c
% All mole fractions in the sequence: MetOH/TG/DG/MG/GL/BD:
xx = [x,1-sum(x,2)];

% Vapour-liquid equilibria (multicomponent ideal VLE, Stichlmair-Fair, 'Distillation', p. 36, 1998)
y = (x*diag(alpha1(1:NC-1)))./((x*(alpha1(1:NC-1) - 1)' + 1)*ones(1,NC-1));

%% Density calculations
% Assuming constant temperature:
T = 91.27+273.15;           % [K]

% Molar density
rho = AR(:,1)./AR(:,2).^(1+(1-(T./AR(:,3))).^AR(:,4)); % [kmol/m^3]
N = (M*ones(1,NC)).*xx;     % molar holdup of components per tray [kmol]

% Volume fractions:
V = N./(ones(NT,1)*rho');   % volume of component per tray [m3]
Vtot = sum(V,2);            % total volume per tray [m3]
ve = V./(Vtot*ones(1,NC));  % volume fraction [m3/m3]
rho_l = sum((ones(NT,1)*rho').*ve,2);   % Common liquid density for the tray [kmol/m3]
 
% Weight fractions:
wt = xx.*(M*ones(1,NC)).*(ones(NT,1)*Mw);
wt_sum = sum(wt,2);
x_wt = wt./(wt_sum*ones(1,NC));

% Common liquid density for each tray [kg/m3]:
r_wt = (ones(NT,1)*rho').*x_wt.*(ones(NT,1)*Mw);
rho_L = sum(r_wt,2);


vol= volume_plate*ones(NT,1);


%% Dynamics
% Vapor flows assuming constant molar flows: [kmol/h]
V = VB*ones(NT-1,1); V(NF:NT-1) = V(NF:NT-1) + (1-qF)*F;

% Liquid flows are given by Franci's Weir Formula
Muw = volume_plate*rho_l;               % Molar holdup under weir [kmol]
Mw_avg = sum(xx.*(ones(NT,1)*Mw),2);    % Average molecular weight per tray [kg/kmol]

% Constants for the Francis weir equation:
Kbf = 60*rho_l*lw.*(10^3./(750*area_plate*rho_l)).^1.5;

L(2:NF,1)      = (Kbf(2:NF,1).*(max(M(2:NF) - Muw(2:NF,1),0)).^1.5);             % Liquid flows below feed (Kmol/min)
L(NT,1)        = LT;                               % Condenser's liquid flow (Kmol/min)

%Liquid control in the refulx drum and reboiler: flow equations for D and B
D = max(Ds + (-M(NT) + MDs)*KcD,0);      % Distillate flow
B = max(Bs + (-M(1) + MBs)*KcB,0);       % Bottoms flow


%% Kinetics

%          A          |   e^(-Ea/R/T)      | [min] |[x->kmol]          
k1 = 3.9      *(10^7)   *exp(-54998/(R*T))  *60  ./vol.*M.^2;
k2 = 5.78     *(10^5)   *exp(-41555/(R*T))  *60  ./vol.*M.^2;
k3 = 5.906    *(10^12)  *exp(-83094/(R*T))  *60  ./vol.*M.^2;
k4 = 9.888    *(10^9)   *exp(-61249/(R*T))  *60  ./vol.*M.^2;
k5 = 5.335    *(10^3)   *exp(-26865/(R*T))  *60  ./vol.*M.^2;
k6 = 2.1      *(10^4)   *exp(-40116/(R*T))  *60  ./vol.*M.^2;



r(:,2) = -(k1.*xx(:,2).*xx(:,1)) + (k2.*xx(:,3).*xx(:,6));
r(:,3) =  (k1.*xx(:,2).*xx(:,1)) - (k2.*xx(:,3).*xx(:,6)) - (k3.*xx(:,3).*xx(:,1)) + (k4.*xx(:,4).*xx(:,6));
r(:,4) =  (k3.*xx(:,3).*xx(:,1)) - (k4.*xx(:,4).*xx(:,6)) - (k5.*xx(:,4).*xx(:,1)) + (k6.*xx(:,5).*xx(:,6));
r(:,6) =  (k1.*xx(:,2).*xx(:,1)) - (k2.*xx(:,3).*xx(:,6)) + (k3.*xx(:,3).*xx(:,1)) - (k4.*xx(:,4).*xx(:,6)) + (k5.*xx(:,4).*xx(:,1)) - (k6.*xx(:,5).*xx(:,6));
r(:,5) =  (k5.*xx(:,4).*xx(:,1)) - (k6.*xx(:,5).*xx(:,6));
r(:,1) = -(k1.*xx(:,2).*xx(:,1)) + (k2.*xx(:,3).*xx(:,6)) - (k3.*xx(:,3).*xx(:,1)) + (k4.*xx(:,4).*xx(:,6)) - (k5.*xx(:,4).*xx(:,1)) + (k6.*xx(:,5).*xx(:,6));
r_tot = r*ones(6,1);

%% Molar balances:
% Time derivatives from material balances for: (1) total holdup and (2) component holdup
% Column
j = 2:NT-1; dMdt(j,1)  = L(j+1) - L(j) + V(j-1) - V(j) + r_tot(j);
i = 1:NC-1; dMxdt(j,i) = diag(L(j+1))*x(j+1,i) - diag(L(j))*x(j,i) + diag(V(j-1))*y(j-1,i) - diag(V(j))*y(j,i) + r(j,i);

% Correction for feed at the feed stage: The feed is assumed to be mixed into the feed stage
dMdt(NF)    = dMdt(NF) + F;
dMxdt(NF,:) = dMxdt(NF,:) + F*zF(1,1:NC-1);

% Reboiler (assumed to be an equilibrium stage)
dMdt(1)     = L(2) - V(1) - B;
dMxdt(1,i)  = L(2)*x(2,i) - V(1)*y(1,i) - B*x(1,i);

% Total condenser (no equilibrium stage)
dMdt(NT)    = V(NT-1) - L(NT) - D;
dMxdt(NT,i) = V(NT-1)*y(NT-1,i) - L(NT)*x(NT,i) - D*x(NT,i);

% Compute the derivative for the mole fractions from d(Mx) = x dM + M dx
dxdt = diag(M)\(dMxdt - diag(dMdt)*x);

% Rearrange elements of composition vector derivatives (dxdt)
dxdt = dxdt(:);

post_processing