function  [xprime,B,D,xB,xD] = column(X,u,Dd,param)

% Simplified version of the multicomA model from skoge.ntnu.no.
% This version made March 2011 by Magnus G. Jacobsen.

%The model has NC*NT states: For each tray the total molar holdup and the
%fraction of the NC-1 lightest components are included.

column_parameters % Call the parameters that are specific to the distillation column

NC = param.NC ;
NT = param.NT ;
NF = param.NF ;
alpha = param.alpha ; % Relative volatilities - lowest for B, the desired product.

% Splitting the states
x = X(1:(NC-1)*NT) ;                   % Liquid compositions from btm to top
M = X((NC-1)*NT+1:NC*NT)  ;             % Liquid hold up from btm to top

% Rearrange elements of composition vector (x) for later use
Iu = reshape(1:(NC-1)*NT,NT,NC-1);
x  = x(Iu)    ;

LT = u(1);                      % Reflux
VB = u(2);                      % Boilup

% Disturbances
F  = Dd(1);                             % Feedrate
qF = Dd(2);                             % Feed liquid fraction
zF = Dd(3:end)  ;                       % Feed compositions

% The distillation model

% Level controllers for column and ocndenser
D = Ds + (M(NT) - MDs)*KcD;      % Distillate flow
B = Bs + (M(1) - MBs)*KcB;       % Bottoms flow

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

% Vapor flows assuming constant molar flows
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 L(i)=K*Mow(i)^1.5 
% Liquid flow L(i) dependent only on the holdup over the weir Mow(i)
% M(i)= Mow(i) + Muw(i) (Total holdup = holdup over weir + holdup below weir) 
L(2:NF,1)      = Kbf*(M(2:NF) - Muw).^1.5;         % Liquid flows below feed (Kmol/min)
L(NF+1:NT-1,1) = Kuf*(M(NF+1:NT-1) - Muw).^1.5;    % Liquid flows above feed (Kmol/min)
L(NT,1)        = LT;                               % Condenser's liquid flow (Kmol/min)

% 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);
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);

% 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) - LT - 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 = inv(diag(M))*(dMxdt - diag(dMdt)*x);
dxdt =(diag(M))\(dMxdt - diag(dMdt)*x);

% Rearrange elements of composition vector derivatives (dxdt)
dxdt = dxdt(:);
xprime = [dxdt ; dMdt] ;

% Isolate top and bottom compositions for later use
xB = [x(1,:) 1-sum(x(1,:))]' ;
xD = [x(end,:) 1-sum(x(end,:))]' ;
    