function F=f(t,y)
% This is file flashsrk.m 
% INPUT: Time t and state vector y
% OUTPUT: Right hand side of DAE set: M dy/dt = f(y)
% States n=y(1); n1=y(2); U=y(3); T=y(4); x1=y(5), y1=y(6); p=y(7):
% Usage with ode15s:   
%   options = odeset('mass',diag([1 1 1 0 0 0 0 0]));
%   y0 = [100.4e3 48.1e3 5.02e8 344.55 50.17 0.4793 0.6 1.1e5]
%   [t,y]= ode15s(@flashsrk,[0 100],y0,options)   
% Plot temperature: 	 plot(t,y(:,4))
 
% All in SI units 
liquid=1;vapor=2;
% I. Data (parameters and independent variables)
R=8.13;                                % J/mol K

% Feed data: flow F (from q), composition z1 and enthalpy (hf from Tf assuming liquid)
F = 1000                               % mol/s
z1=0.5;                                % mol1/mol
Tf=1.0*300;                            % K  (increase from 300K to 340K)
pf=2.e5  % N/m2
[Zfl,phifl,hfl,Vflm]=srksis([z1 1-z1 0 0],Tf,pf,liquid)      
[Zfg,phifg,hfg,Vfgm]=srksis([z1 1-z1 0 0],Tf,pf,vapor)       
qlf = 0.9;                                % Fraction liquid
hf = qlf*hfl + (1-qlf)*hfg

% 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=y(1)                                % total holdup in tank (both phases)       [mol]  
n1=y(2)                               % component 1 holdup in tank (both phases) [mol1]
U=y(3)                                % total internal energy (both phases)      [J]
T=y(4)                                % temperature (same in both phases)        [K]
x1=y(5)                               % liquid composition
y1=y(6)                               % vapor composition
p =y(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 = max(Vtot/1000 , (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];  



