function [e,a,q,td]=wilg(t,f,p)

%calculate humidity parameters
%input:
% t- temperature in C
% f- relative humidity
% p- air pressure
% output:
% e- pressure of water vapor
% a- absolute humidity in g/m^3
% q- specific humidity in g/kg
% td- dew point in C

load cis_pary.dat;

T=cis_pary(:,1);
Ew=cis_pary(:,2);
El=cis_pary(:,3);

if t>=0
  E=interp1(T,Ew,t);
 elseif (t<0)&(t>-50)
  E=interp1(T,El,t);
 else
  E=El(1)*exp(-2.5*1.e6/461.5*(1./(273.16+t)-1./(273.16-50)));
end 

e=E.*f./100;
a=217*e./(t+273.16);

if t>=0
 td=interp1(Ew,T,e);
elseif (t<0)&(t>-50);
 td=interp1(El,T,e);
else
 td=nan;
end

q=622*e./p;  
