function dxdt=f(t,x)
% This is file distpdyn.m (NOT USED IN BOOK)
% It simulates the open-loop pressure dynamics in a distillation column, which we see are very fast (no surprise!)
% Boundary conditions: Given flows in and out (given F0=boilup and FNT=condensation rate)
% Usage
% NT=40; x0 = 15*eye(1,NT); % initial states (not steady-state)
%[T,X] = ode15s('distpdyn',[0 100],x0) % First simulate to t=100 (steady-state) with ALT-expression for F(NT)
%  xss = X(length(X),:) % Save the steady-state mole fractions
  % Then change to F(NT)=FNT (SIM) and run simulations
%  [T,X] = ode15s('distpdyn',[0 0.1],xss) % Run new simulation (e.g change F0)
%  Vol=0.35; p = (8.31*300/Vol)/1.e5*X;
%  plot(T,p), figure(2), plot(T,p(:,NT))

% I. Data
A=0.7; % m2, column area 
h=0.5; % m, tray distance
V=A*h; % m3, per tray
R = 8.31; % J/mol,K
T = 300; % K

kpf=1; % Factor kpf for pressure drop (dp) in tray holes 
%       (not included liquid height dp here, which would be a much slower effect through change in L) 
kp=kpf*A; % mol/s / N/m2 
p0 = 1e5; % = 1 bar top pressure (initial steady-state)
NT=40;
% Flows in and out
F0ss = 40; % mol/s (boilup btm). 
%F0 = F0ss  % ALT
F0 = 1.1*40; % SIM
FNT = F0ss; % mol/s (condensing top)

% II. Extract present values of states
n = x;

% IIIa. Intermediate calculations
p=(R*T/V)*n;

i=1:NT-1;
F(i) = kp*(p(i)-p(i+1)) ;
%F(NT) = kp*(p(NT)-p0) ; % ALT. use this expression when finding initial steady state
F(NT)=FNT;  % SIM

% Derivatives
dndt(1)=F0  -F(1);
i=2:NT;
dndt(i)=F(i-1)-F(i);


dxdt=dndt';

