function [c,ceq]=f(y)
%function [c,ceq]=f(y)
% This is file flashss.m
% INPUT: Initial state vector y
% OUTPUT: c is not used, ceq = value of equality constraints
 
% All in SI units 
% I. Data (parameters and independent variabless)
% Data for: 1-methanol, 2-ethanol
A1=8.08097; B1=1582.271; C1=239.726;   % Antoine psat1 [mmHg] w/ T [C] (T-range: 15C - 84C)
A2=8.11220; B2=1592.864; C2=226.184;   % Antoine psat2 [mmHg] w/ T [C] (T-range: 20C - 93C)
cpl1=80; cpl2=131; cpv1=44; cpv2=65;   % heat capacity        [J/K mol]
T0=298.15; hvap01=38000; hvap02=43000; % heat of vap. at T0   [J/mol]
Vl1= 40.7e-6; Vl2=58.7e-6;             % liquid molar volumes [m3/mol]
R=8.13;                                % J/mol K

% Feed data
q=0.1;                                 % m3/s
z1=0.5;                                % mol1/mol
Tf=1.1*400;                            % K  (increase from 400K to 440K)
p0=1e5;                                % N/m2  (=1 bar downstream pressure)
Vf = z1*Vl1 + (1-z1)*Vl2;              % m3/mol (molar volume feed)
F = q/Vf;                              % mol/s (feed rate)

% Total tank volume 
Vtot = 10;                             % m3

% Valve constant and controller gain
kg=0.1; kl=100e3;                      % 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]
VL=y(5);                               % liquid volume                            [m3]

% IIIa. Intermediate calculations
% VLE
x1 = n1/n;
p1s=10.^(A1-B1/(T-273.15+C1))/750e-5;  % psat1 from Antoine                    [N/m2]
p2s=10.^(A2-B2/(T-273.15+C2))/750e-5;  % psat2 from Antoine                    [N/m2]
p1 = x1*p1s;                           % partial pressure component 1          [N/m2]
p2 = (1-x1)*p2s;                       % partial pressure component 2          [N/m2]
p = p1+p2;                             % pressure = sum of partial pressures   [N/m2]
y1 =  p1/p;                            % vapor fraction component 1            [mol1/mol]
Vm  = x1*Vl1 + (1-x1)*Vl2;             % molar volume (liquid phase)           [m3/mol]

% Phase distibution  
VG = Vtot - VL;                        % gas volume                            [m3]
nL = VL / Vm;                          % liquid holdup                         [mol]
nG = p*VG / (R*T);                     % gas holdup  (ideal gas law)           [mol]

% Enthalpies [J/mol] (Ref.state: pure liquid at T0)
hF = [z1*cpl1 + (1-z1)*cpl2] * (Tf-T0);
hL = [x1*cpl1 + (1-x1)*cpl2] * (T -T0);  
hG = [y1*cpv1 + (1-y1)*cpv2] * (T -T0) + y1*hvap01 + (1-y1)*hvap02;

% Vapor and liquid flow 
G = kg*(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: Overall mass balance
f2 = F*z1 - G*y1 - L*x1;         % =dn1/dt Dynamic: Component 1 mass balance
f3 = F*hF - G*hG - L*hL;         % =dU/dt  Dynamic: Energy balance
f4 = U + p*Vtot - hL*nL - hG*nG; % = 0     Algebraic: Internal energy U (can find T from this)
f5 = n - nG - nL;                % = 0     Algebraic: Total holdup n    (can find VL from this)
F  = [f1; f2; f3; f4; f5];  
c=[]; ceq=F;
