function DAE = HEX_VL(t,w,param)
% April 2017. All given data are passed as parameters from the main script.
%Compares the liquid volume in each cell with a tolerance to select the
% valid set of eqautions
% I. Given data
M  = param(1);            % Cells numbers
NC = param(2);            % Number of components (1)Nitrogen - (2)Methane - (3)Ethane - (4)Propane 
V  = param(3);            % Total volume of HEX

% Data common for both phases
F0 = param(4);             % Feed rate [kmol/min]   
P0 = param(5);             % Inlet pressure, MPa
DP = param(6);             % Pressure drop along HEX, MPa
T0 = param(7);             % MR inlet temperature, K
z0 = param(15:15+NC-1);    % Specified composition [kmol1/kmol]
Qmax = param(8);           % Given heat, [MJ/min]
t1 = param(9);             % time when the heat starts increasing from 0, [min]
t2 = param(10);            % time when the heat stops increases and becomes constant, [min]
Tdew = param(end-1);       % Dew point, [K]
Tbub = param(end);         % Bubble point, [K]
Vcell = V/M;               % volume of each cell, [m3]

% Data for single phase (vapor or liquid)    
Z1P = param(11);                        % from FindZ.m
[~, h01P] = srkHEXV3(Z1P,z0,T0,P0);     % Specific enthalpy for vapor[MJ/kmol] 

% Data for vap-liquid phase
x0 = param(15+NC:15+2*NC-1);             % liquid composition for MR at -15 degC from Hysys
y0 = param(15+2*NC:15+3*NC-1);           % vapor composition from MR at -15 deg C from Hysys
ZG0 = param(12);                         % from FindZ
ZL0 = param(13);                         % from FindZ
v0 = param(14);                          % inlet vapor phase from Hysys                                            

[~, hfG] = srkHEXV3(ZG0,x0,T0,P0);        % Specific enthalpy for vapor[MJ/kmol]
[~, hfL] = srkHEXV3(ZL0,y0,T0,P0);        % Specific enthalpy for liquid[MJ/kmol]                                           
h02P = v0*hfG + (1-v0)*hfL;               % Specific feed enthalpy [MJ/mol]

% Given heat as a function of time, [[MJ/min]
   function Q = Q(t)
        if t < t1
            Q = 0;
        elseif t >= t1 && t < t2
            Q = Qmax/(t2-t1) * (t-t1);
        else
            Q = Qmax;
        end
    end

Qcell = Q(t)/M; % equal heat removed from each cell

% II. States. Same for both regions. 
  nT = w(1:M)';                                          % total holdup in  (both phases)       [kmol]  
  n(1:NC*M) = w(M+1:M+M*NC)';                            % component 1 holdup in tank (both phases) [kmol1]
  U = w(M+NC*M+1:2*M+NC*M)';                             % total internal energy (both phases)      [MJ]
  T = w(2*M+NC*M+1:3*M+M*NC)';                           % temperature (same in both phases)        [K]
  x(1:NC*M) = w(3*M+M*NC+1:3*M+2*M*NC)';                 % liquid composition
  y(1:NC*M) = w(3*M+2*M*NC+1:3*M+3*M*NC)';               % vapor composition
  VL = w(3*M+3*M*NC+1:4*M+3*M*NC)';                      % liquid volume [m3]
  Zg = w(4*M+3*M*NC+1:5*M+3*M*NC)';                      % compress  gas [-]
  Zl = w(5*M+3*M*NC+1:6*M+3*M*NC)';                      % compress liquid [-]
  P = w(6*M+3*M*NC+1:7*M+3*M*NC)';                       % pressure, [MPa] 

% III. System of equations 

% IIIa.The outlet flow is calculated with a linear valve equation
F = zeros(1,M);      % [kmol/min],preallocating memory        
k = F0/DP*(M+2);     % [kmol/min*MPa], valve constant, considering equal pressure drop in each cell
Pout = P0-DP;        % [MPa], outlet pressure

for j=1:M-1
    F(j) = k*(P(j)-P(j+1)); % Outlet flow
end

for j=M
    F(j) = k*(P(j)-Pout);  
end

% IIIb. SRK. Compute fugacity, enthalpy, molar volumes
% preallocating memory
phil = zeros(M,NC); phig = zeros(M,NC); hl = zeros(1,M); hg = zeros(1,M);
Vlm = zeros(1,M); Vgm = zeros(1,M); Al = zeros(1,M); Ag = zeros(1,M); Bl = zeros(1,M);Bg = zeros(1,M);

for j=1:M
    
    [phil_temp,hl_temp,Vlm_temp,Al_temp,Bl_temp]= srkHEXV3(Zl(j), [x(j*NC-3) x(j*NC-2) x(j*NC-1)  x(j*NC)], T(j),P(j));    % liquid
    [phig_temp,hg_temp,Vgm_temp,Ag_temp,Bg_temp]= srkHEXV3(Zg(j), [y(j*NC-3) y(j*NC-2) y(j*NC-1) y(j*NC)], T(j), P(j));    % vapor
    
    phil(j,1)= phil_temp(1);  phil(j,2)= phil_temp(2); phil(j,3)= phil_temp(3);  phil(j,4)= phil_temp(4);
    phig(j,1)= phig_temp(1);  phig(j,2)= phig_temp(2); phig(j,3)= phig_temp(3);  phig(j,4)= phig_temp(4);
    hl(j)= hl_temp; Vlm(j)= Vlm_temp; Al(j)= Al_temp; Bl(j)= Bl_temp;
    hg(j)= hg_temp; Vgm(j)= Vgm_temp; Ag(j)= Ag_temp; Bg(j)= Bg_temp; 
    
end

K = phil./phig;  % VLE constant

% Phase distribution 
VG = zeros(1,M); nG = zeros(1,M);
f1 = zeros(1,M); f2 = zeros(1,M*(NC)); f3 = zeros(1,M); f4 = zeros(1,M); f5 = zeros(1,M*NC); 
f6 = zeros(1,M*NC); f7 = zeros(1,M); f8 = zeros(1,M); f9 = zeros(1,M); f10 = zeros(1,M); % preallocation memory

% Checking of location in the phase region, e.g. 1 phase (vapor or liquid)
% or 2 phase (vapor-liquid)

for j=1  % cell # 1. 
   
 % Compare the liquid volume in each cell with a given tolerance
    if  VL(j) < 1e-5                   
        VG(j) = Vcell;               % gas volume only for vapor phase [m3]
        nG(j) = VG(j)/Vgm(j);        % gas holdup                            [kmol]
    
 % IIIc. System of equations for single phase   
    f1(j) = F0-F(j);                                                             
    f2((NC)*j-3:(NC)*j) = F0*z0 - F(j)*y(NC*j-3:NC*j);
    f3(j) = F0*h01P - F(j)*hg(j)- Qcell;                              
    f4(j) = U(j) + P(j)*Vcell - hg(j)*nG(j);    
    f5((NC)*j-3:(NC)*j) = y(NC*j-3:NC*j) - n(NC*j-3:NC*j)./nT(j);    
    f6((NC)*j-3:(NC)*j) = x(NC*j-3:NC*j) - x0(1:4);
    f7(j) = VL(j);
    f8(j) = Zg(j)^3-Zg(j)^2+Zg(j)*(Ag(j)-Bg(j)-Bg(j)^2)-Ag(j)*Bg(j);  
    f9(j) = Zl(j) - ZL0;
    f10(j) = nT(j)-nG(j);   
    
    else % system is in 2 phase region
    % phase distribution for 2-phase region
        VG = Vcell - VL;      % gas volume for 2 phase region
        nL = VL./Vlm;         % liquid holdup                         [kmol]
        nG = VG./Vgm;         % gas holdup                            [kmol]
        v = nG./nT;           % vapor fraction
    
    % IIId. System of equations for 2 phase region    
    
    K1=K(j,:); % K-values
        
    f1(j) = F0-F(j);                                                             
    f2(NC*j-3:NC*j) = F0*v0*y0+ F0*(1-v0)*x0 - F(j)*v(j)*y(NC*j-3:NC*j)-F(j)*(1-v(j))*x(NC*j-3:NC*j);  
    f3(j) = F0*h02P - F(j)*v(j)*hg(j) - F(j)*(1-v(j))*hl(j) - Qcell;  
    f4(j) = U(j) + P(j)*Vcell - hl(j).*nL(j) - hg(j).*nG(j);                     
    f5(NC*j-3:NC*j) = n(NC*j-3:NC*j) - x(NC*j-3:NC*j).*nL(j) - y(NC*j-3:NC*j).*nG(j);     
    f6(NC*j-3:NC*j) = y(NC*j-3:NC*j) - K1(1:NC).*x(NC*j-3:NC*j);                           
    f7(j) = nT(j)-nG(j)-nL(j); 
    f8(j) = (Zg(j)^3-Zg(j)^2+Zg(j)*(Ag(j)-Bg(j)-Bg(j)^2)-Ag(j)*Bg(j));                                       
    f9(j) = (Zl(j)^3-Zl(j)^2+Zl(j)*(Al(j)-Bl(j)-Bl(j)^2)-Al(j)*Bl(j));
    f10(j)= sum(x(j*NC-3:j*NC)) - sum(y(j*NC-3:j*NC));
    
    end % if statement
     G = roots([1 -1 Ag(j)-Bg(j)-Bg(j)^2 -Ag(j)*Bg(j)]);
    Zg(j);
    L = roots([1 -1 Al(j)-Bl(j)-Bl(j)^2 -Al(j)*Bl(j)]);
    Zl(j);
end % for loop

for j=2:M   % Cells 2 to N   
    
    % check if liquid volume is below a certain tolerance
    if  VL(j)< 1e-5              
        VG(j) = Vcell;           % gas ocupies entire cell volume [m3]
        nG(j) = VG(j)/Vgm(j);    % gas holdup                     [kmol]
        
    f1(j) = F(j-1)-F(j);                                                  % overall mass balance           
    f2((NC)*j-3:(NC)*j) = F(j-1)*y(NC*j-7:NC*j-NC)- F(j)*y(NC*j-3:NC*j);  % component mass balance
    f3(j) = F(j-1)*hg(j-1) - F(j)*hg(j)- Qcell;                           % energy balance 
    f4(j) = U(j) + P(j)*Vcell - hg(j)*nG(j);                              % internal energy
    f5((NC)*j-3:(NC)*j) = y(NC*j-3:NC*j) - n(NC*j-3:NC*j)./nT(j);         % vapor composition
    f6((NC)*j-3:(NC)*j) = x(NC*j-3:NC*j) - x0(1:4);                       % "liquid composition", dummy eq
    f7(j) = VL(j);                                                        % "liquid volume", dummy eq
    f8(j) = Zg(j)^3-Zg(j)^2+Zg(j)*(Ag(j)-Bg(j)-Bg(j)^2)-Ag(j)*Bg(j);      % gas comrpesibility
    f9(j) = Zl(j) - ZL0;                                                  % "liquid comrpesibility", dummy eq
    f10(j) = nT(j)-nG(j);                                                 % total holdup
    
    else % system is in 2 phase region
    % phase distribution for 2-phase region
        VG = Vcell - VL;      % gas volume for 2 phase region
        nL = VL./Vlm;         % liquid holdup                         [kmol]
        nG = VG./Vgm;         % gas holdup                            [kmol]
        v = nG./nT;           % vapor fraction

% IIId. System of equations for 2 phase region      

    K1=K(j,:);  % K-values for all components
            
    f1(j) = F(j-1)-F(j);                                                                                        % overall mass balance
    f2(NC*j-3:NC*j) = F(j-1)*v(j-1)*y(NC*j-7:NC*j-NC)+F(j-1)*(1-v(j-1))*x(NC*j-7:NC*j-NC) - ...                 % component mass balance
                      F(j)  *v(j)  *y(NC*j-3:NC*j)   -F(j)  *(1-v(j))  *x(NC*j-3:NC*j);                         
    f3(j) = F(j-1)*(1-v(j-1))*hl(j-1) + F(j-1)*v(j-1)*hg(j-1) - F(j)*(1-v(j))*hl(j) - F(j)*v(j)*hg(j) - Qcell;  % energy balance                            
    f4(j) = (U(j) + P(j)*Vcell- hl(j)*nL(j) - hg(j)*nG(j));                                                     % internal energy
    f5(NC*j-3:NC*j) = n(NC*j-3:NC*j) - x(NC*j-3:NC*j).*nL(j) - y(NC*j-3:NC*j).*nG(j);                           % component holdup
    f6(NC*j-3:NC*j) = y(NC*j-3:NC*j) - K1(1:NC).*x(NC*j-3:NC*j);                                                % VLE
    f7(j)= nT(j)-nG(j)-nL(j);                                                                                   % total holdup
    f8(j) = (Zg(j)^3-Zg(j)^2+Zg(j)*(Ag(j)-Bg(j)-Bg(j)^2)-Ag(j)*Bg(j));                                          % gas compresibility
    f9(j) = (Zl(j)^3-Zl(j)^2+Zl(j)*(Al(j)-Bl(j)-Bl(j)^2)-Al(j)*Bl(j));                                          % liquid compresibilty
    f10(j)= sum(x(j*NC-3:j*NC)) - sum(y(j*NC-3:j*NC));                                                          % 1-sum(x) - (1-sum(y))
    
    end % if statement
    G = roots([1 -1 Ag(j)-Bg(j)-Bg(j)^2 -Ag(j)*Bg(j)]);
    Zg(j);
    L = roots([1 -1 Al(j)-Bl(j)-Bl(j)^2 -Al(j)*Bl(j)]);
    Zl(j);
end % for statemant

DAE  = [f1'; f2';f3'; f4'; f5'; f6'; f7';f8';f9';f10']; % System of DAE 

end % end function structure
        