function [theta, thetae, q, qsat, qsati] = thermo_rh(t,p,rh)
% thermo_rh generates thermodynamic variables from t(oC) p(mb) rh(%) vectors
% output [theta, thetae, q, qsat, qsati] = thermo_rh(t,p,rh)
% in K, K, g/kg, g/kg, g/kg

%convert t,p to SI units
 tk = t + 273.15; %K
 p = p*100; %Pa

%calculate the potential temperature from temperature array
 p0=1000*100;	%reference pressure in Pa
 R=287;		%gas constant
 cp=1004;	%specific heat wrt pressure
 K= R/cp;
 a = ((p./p0).^(-K));
 theta = tk.*a;

%calculate equivalent potential temperature 
 Lv = 2.5e6; 	%latent heat of vapourisation
 eps = 0.622; 	%Rd/Rv
 es0 = 0.61e3;	%reference saturation vapour pressure Pa
 tk0 = 273.15;	%reference temperature
% es_old = es0*exp(((eps*Lv)/R)*(1./tk0 - 1./tk)) %saturation vapour pressure
 [es esi] = thermo_es(t);   			%...wrt to water and ice 
 es = es*100; esi = esi*100; 		%convert to Pa
 rsat = eps*es./p;		 	%saturation mixing ratio
 rsati = eps*esi./p;		 	%sat mixing ratio wrt ice
 qsat = rsat./(rsat+1);				%sat specific humidity
 qsati = rsati./(rsati+1);			%sat specific humidity wrt ice
 thetae = theta.*exp((Lv.*rsat)./(cp.*tk));	%equivalent potential temp.

%calculate water vapour mixing ratio r and q specific humidity
 r = rsat.*(rh./100);
 q = qsat.*(rh./100);

 q = q*1e3;
 qsat = qsat*1e3;
 qsati = qsati*1e3;