function F=f(t,w)
% This is file flashsrk.m (DYNAMIC FLASH with given energy and volume)
% INPUT: Time t and state vector w
% OUTPUT: Right hand side of DAE set: M dw/dt = f(w)
% States n=w(1); n1=w(2); U=w(3); T=w(4); x1=w(5), y1=w(6); p=w(7):
% Usage with ode15s:   
%   options = odeset('mass',diag([1 1 1 0 0 0 0]));
%   w0 = [52.5e3 15.7e3 -1.315e9 250 0.297 0.581 1.07e5]    % initial state
%   [t,w]= ode15s(@flashsrk,[0 100],w0,options)             % integrate
%   plot(t,w(:,4))                                          % plot temperature
 
% All in SI units 
liquid=1;vapor=2;
% I. Data (parameters and independent variables)
R=8.13;                                % J/mol K

% Given feed data: flow F, composition z1 and enthalpy (hf from Tf assuming liquid)
F = 1000;                              % Specified feed rate [mol/s]
z1=0.5;                                % Specified composition [mol1/mol]
hf = -10.e3;        			       % Specified enthalpy [J/mol] (ref: h=0 for ideal gas at 298K)

% Other data 
Vtot = 10;                             % m3    (total tank volume)
p0=1*1e5;                                % N/m2  (=1 bar downstream pressure)

% Valve constant and controller gain
kg=0.1; kl=10000;                      % note that the P-controller gain kl is large

% II. Extract present value of states
n=w(1);                                % total holdup in tank (both phases)       [mol]  
n1=w(2);                               % component 1 holdup in tank (both phases) [mol1]
U=w(3);                                % total internal energy (both phases)      [J]
T=w(4);                                % temperature (same in both phases)        [K]
x1=w(5);                               % liquid composition
y1=w(6);                               % vapor composition
p =w(7);                               % pressure [N/m2]

% IIIa. Intermediate calculations

% Fugacity, enthalpy and volumes from SRK
[Zl,phil,hl,Vlm]=srksis([x1 1-x1 0 0],T,p,liquid);        % liquid
[Zg,phig,hg,Vgm]=srksis([y1 1-y1 0 0],T,p,vapor);        % vapor

K = phil./phig;

% Phase distribution  
VG = (Vtot-n*Vlm) * Vgm/(Vgm-Vlm);  % liquid volume                         [m3]
VL = Vtot - VG;                       % gas volume                            [m3]
nL = VL / Vlm;                          % liquid holdup                         [mol]
nG = VG / Vgm;                          % gas holdup                            [mol]

% Vapor and liquid flow 
G = kg*max(0,(p-p0));   % simple valve equation for outflow of gas              [mol/s]
VLs = Vtot/2;           % Setpoint level (volume): keep 50% liquid in tank      [m3]
L = kl*(VL-VLs);        % Level controller with proportional gain kl            [mol/s]

% IIIb. Evaluate right hand side of DAE-set: M dy/dt = f(y)
f1 = F - G - L;                  % =dn/dt  Dynamic (1): Overall mass balance
f2 = F*z1 - G*y1 - L*x1;         % =dn1/dt Dynamic (NC-1): Component 1 mass balance
f3 = F*hf - G*hg - L*hl;         % =dU/dt  Dynamic (1): Energy balance
f4 = U + p*Vtot - hl*nL - hg*nG; % = 0     Algebraic (1): Internal energy U (can find T from this)
f5 = n1 - x1*nL - y1*nG;         % = 0     Algebraic (NC-1): Component holdup n1   (can find x1 from this)
f6 = y1 - K(1)*x1;               % = 0	   Algebraic (NC-1): VLE for comp. 1 (can find y1 from this)
f7 = (1-y1) - K(2)*(1-x1);       % = 0	   Algebraic (1): VLE for last component (can find p from this)
F  = [f1; f2; f3; f4; f5; f6; f7];  


