function hr=h(T,vf)
% Compute enthalpy
% T : Temperature [K]
% vf: Vapour fraction (1=vapour and 0=liquid)
Tref=273.15;      % Reference temperature [K]

HeatPropertyConstants;

%Components (Cmp): H G P N A where
% H: Hydrogen
% G: Equal mix of ethane, propane and butane
% P: Hepthane
% N: Cycloheptane
% A: Toluene

% Liquid phase
CplCmp=[Cpl(1) (Cpl(2)+Cpl(3)+Cpl(4))/3 Cpl(6) Cpl(7) Cpl(8)]';
hlCmp =CplCmp.*(T-Tref);                                          % [J/mol]

% Enthalpy of vaporization
hlvCmp=[hlv(1) (hlv(2)+hlv(3)+hlv(4))/3 hlv(6) hlv(7) hlv(8)]';   % [J/mol]

% Gas properties 
% hv=Integral from Tref to T of Cp(T)dT
Tp    =[T^6 T^5 T^4 T^3 T^2 T]';
Trefp =[Tref^6 Tref^5 Tref^4 Tref^3 Tref^2 Tref]';
Icons =[1/6 1/5 1/4 1/3 1/2 1];
I     =ones(size(Cpv(:,1)));
hv    =I*Icons.*Cpv*(Tp-Trefp);
hvCmp =[hv(1) (hv(2)+hv(3)+hv(4))/3 hv(6) hv(7) hv(8)]';          % [J/mol]

% Compute enthalpy
if T >= Tref
  hr=[(1-vf)*hlCmp + vf*(hlvCmp+hvCmp)];                          % [J/mol]
else % if T<Tref use constant Cp for vapour phase
  hr=[(1-vf)*hlCmp + vf*(hlvCmp+dh_dT(Tref,vf)*(T-Tref))];        % [J/mol]
%  warning(['Enthalpy calculation h(T): T<Tref T=' num2str(T)])
end
return
