function dhr=dh_dT(T,vf)
% Compute derivative of enthalpy with respect to T
% T : Temperature [K]
% vf: Vapour fraction (1=vapour and 0=liquid)
% Note: Pentane currently left out as a component

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)]';
dhlCmp=CplCmp;                                          % [J/mol]    Specific enthalpy of liquid
 
% Heat of vaporization at is constant
% dhlv_dt=0

% Gas phase 
Tp    =[T^5 T^4 T^3 T^2 T^1 1]';

dhv   = Cpv*Tp;
dhvCmp=[dhv(1) (dhv(2)+dhv(3)+dhv(4))/3 dhv(6) dhv(7) dhv(8)]';

if T >= Tref
  dhr=[(1-vf)*dhlCmp + vf*dhvCmp];                   % [J/mol/K]
else
  dhr=[(1-vf)*dhlCmp + vf*dh_dT(Tref,vf)];           % [J/mol/K]
end
return

