function F = reactivedist_c_DAE(t,Y)
% Essential parameters:
parameters_c

% Inputs and disturbances:
LT = 226.8/60;          % Reflux [kmol/min]
VB = 309.4/60;          % Boilup [kmol/min]

% if t >= 2
%         F = 355/60;
%         zF = [300/355 55/350 0 0 0 0]; 
%     else
%         F = 350/60;
%         zF = [300/350 50/350 0 0 0 0]; 
% end
F  = (350/60);          % Feedrate [kmol/min]
qF = 1.064;             % Feed liquid fraction
zF = [300/350 50/350 0 0 0 0]; % Feed compositions for component x1 to xNC

% II. Extract present value of states
x0 = Y(1:NT*(NC-1),1);          % Compositions of compounds 1-5 on each tray
M0 = Y((NT*(NC-1))+1:NT*NC,1);  % Molar holdup on each tray [kmol]  
T0 = Y((NT*NC)+1:(NT*NC)+NT,1); % Temperature on each tray [K] 
x = reshape(x0,NT,NC-1);        % Reshaping the  compositions
M = reshape(M0,NT,1);           % Reshaping the molar holdups
TL = reshape(T0,NT,1);          % Reshaping the temperature

% 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


%% UNIQUAC model

% Vapour pressure equation where Psat [mmHg] and temperature T [K]
j = 1:NT;
Ps1(:,1) = (A(1,2)./TL) + (A(1,3)*log(TL)) + A(1,4)*TL.^A(1,5); % TG
Ps1(:,2) = (A(2,2)./TL) + (A(2,3)*log(TL)) + A(2,4)*TL.^A(2,5); % DG
Ps1(:,3) = (A(3,2)./TL) + (A(3,3)*log(TL)) + A(3,4)*TL.^A(3,5); % MG
Ps1(:,4) = (A(4,2)./TL) + (A(4,3)*log(TL)) + A(4,4)*TL.^A(4,5); % BD
Ps1(:,5) = (A(5,2)./TL) + (A(5,3)*log(TL)) + A(5,4)*TL.^A(5,5); % GL
Ps1(:,6) = (A(6,2)./TL) + (A(6,3)*log(TL)) + A(6,4)*TL.^A(6,5); % MetOH
Ps = 760*9.869*(10^-3)*exp((ones(NT,1)*A(:,1)') + Ps1);

%       TG     DG      MG      BD      GL      MetOH
XU = [xx(j,2) xx(j,3) xx(j,4) xx(j,6) xx(j,5) xx(j,1)];

% Calculating the surface area fractions:
Q_tot = (ones(NT,1)*Q).*XU;
Q_sum = sum(Q_tot,2);
F1 = (ones(NT,1)*Q)./(Q_sum*ones(1,6));

% Calculating the volume fraction per mixture mole fraction:
RPt = (ones(NT,1)*RP).*XU;
RP_sum = sum(RPt,2);
PHI = (ones(NT,1)*RP)./(RP_sum*ones(1,NC));

% Calculating matrix of UNIQUAC coefficients TAF_ij:
TAF1 = exp(-W/(Re*TL(1,1)));
TAF2 = exp(-W/(Re*TL(2,1)));
TAF3 = exp(-W/(Re*TL(3,1)));
TAF4 = exp(-W/(Re*TL(4,1)));
TAF5 = exp(-W/(Re*TL(5,1)));

% Calculating the combinatorial contribution:
lnGammaC = (1-PHI+log(PHI))-((z/2)*(ones(NT,1)*Q).*(1-(PHI./F1)+log(PHI./F1)));

% Calculating the residual contribution:
Q1 = ones(NC,1)*Q;          % Q as a (6,6) matrix
NOM(1,:) = sum(Q1.*(ones(NC,1)*XU(1,:)).*TAF1,2)';
DENOM(1,1) = sum(Q.*XU(1,:));
RATIO(1,1) = sum((Q.*XU(1,:).*TAF1(1,:))./NOM(1,1));
RATIO(1,2) = sum((Q.*XU(1,:).*TAF1(2,:))./NOM(1,2));
RATIO(1,3) = sum((Q.*XU(1,:).*TAF1(3,:))./NOM(1,3));
RATIO(1,4) = sum((Q.*XU(1,:).*TAF1(4,:))./NOM(1,4));
RATIO(1,5) = sum((Q.*XU(1,:).*TAF1(5,:))./NOM(1,5));
RATIO(1,6) = sum((Q.*XU(1,:).*TAF1(6,:))./NOM(1,6));
lnGammaR(1,:) = Q.*(1-log(NOM(1,:)./DENOM(1,1))-RATIO(1,:));

NOM(2,:) = sum(Q1.*(ones(NC,1)*XU(2,:)).*TAF2,2)';
DENOM(2,1) = sum(Q.*XU(2,:));
RATIO(2,1) = sum((Q.*XU(2,:).*TAF2(1,:))./NOM(2,1));
RATIO(2,2) = sum((Q.*XU(2,:).*TAF2(2,:))./NOM(2,2));
RATIO(2,3) = sum((Q.*XU(2,:).*TAF2(3,:))./NOM(2,3));
RATIO(2,4) = sum((Q.*XU(2,:).*TAF2(4,:))./NOM(2,4));
RATIO(2,5) = sum((Q.*XU(2,:).*TAF2(5,:))./NOM(2,5));
RATIO(2,6) = sum((Q.*XU(2,:).*TAF2(6,:))./NOM(2,6));
lnGammaR(2,:) = Q.*(1-log(NOM(2,:)./DENOM(2,1))-RATIO(2,:));

NOM(3,:) = sum(Q1.*(ones(NC,1)*XU(3,:)).*TAF3,2)';
DENOM(3,1) = sum(Q.*XU(3,:));
RATIO(3,1) = sum((Q.*XU(3,:).*TAF3(1,:))./NOM(3,1));
RATIO(3,2) = sum((Q.*XU(3,:).*TAF3(2,:))./NOM(3,2));
RATIO(3,3) = sum((Q.*XU(3,:).*TAF3(3,:))./NOM(3,3));
RATIO(3,4) = sum((Q.*XU(3,:).*TAF3(4,:))./NOM(3,4));
RATIO(3,5) = sum((Q.*XU(3,:).*TAF3(5,:))./NOM(3,5));
RATIO(3,6) = sum((Q.*XU(3,:).*TAF3(6,:))./NOM(3,6));
lnGammaR(3,:) = Q.*(1-log(NOM(3,:)./DENOM(3,1))-RATIO(3,:));

NOM(4,:) = sum(Q1.*(ones(NC,1)*XU(4,:)).*TAF4,2)';
DENOM(4,1) = sum(Q.*XU(4,:));
RATIO(4,1) = sum((Q.*XU(4,:).*TAF4(1,:))./NOM(4,1));
RATIO(4,2) = sum((Q.*XU(4,:).*TAF4(2,:))./NOM(4,2));
RATIO(4,3) = sum((Q.*XU(4,:).*TAF4(3,:))./NOM(4,3));
RATIO(4,4) = sum((Q.*XU(4,:).*TAF4(4,:))./NOM(4,4));
RATIO(4,5) = sum((Q.*XU(4,:).*TAF4(5,:))./NOM(4,5));
RATIO(4,6) = sum((Q.*XU(4,:).*TAF4(6,:))./NOM(4,6));
lnGammaR(4,:) = Q.*(1-log(NOM(4,:)./DENOM(4,1))-RATIO(4,:));

NOM(5,:) = sum(Q1.*(ones(NC,1)*XU(5,:)).*TAF5,2)';
DENOM(5,1) = sum(Q.*XU(5,:));
RATIO(5,1) = sum((Q.*XU(5,:).*TAF5(1,:))./NOM(5,1));
RATIO(5,2) = sum((Q.*XU(5,:).*TAF5(2,:))./NOM(5,2));
RATIO(5,3) = sum((Q.*XU(5,:).*TAF5(3,:))./NOM(5,3));
RATIO(5,4) = sum((Q.*XU(5,:).*TAF5(4,:))./NOM(5,4));
RATIO(5,5) = sum((Q.*XU(5,:).*TAF5(5,:))./NOM(5,5));
RATIO(5,6) = sum((Q.*XU(5,:).*TAF5(6,:))./NOM(5,6));
lnGammaR(5,:) = Q.*(1-log(NOM(5,:)./DENOM(5,1))-RATIO(5,:));

% Calculating the activity coefficents:
lnGamma = lnGammaC + lnGammaR;
Gamma = exp(lnGamma);

% Calculating the total pressure:
Ptot = sum(Gamma.*XU.*Ps,2);

%% 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(NC,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(:);

Temp = 760 - Ptot;

F = [dxdt; dMdt; Temp];