function [e,q,rh]=wetbulb(td,tw,p)    % assumes p in mb, t in C
  esw = 6.1121*exp(17.502*tw/(tw+240.97))*(1.0007+3.46e-6*p);% esat at wet bulb temp
  e = esw - 6.51e-4*p*(td-tw);   %psychrometric eqn for e (mb)
  es = 6.112*exp(17.502*td/(td+240.97))*(1.0007+3.46e-6*p);% esat at dry bulb temp
  q = e*621.97/(p-0.378*e);  % q g/kg
  rh = 100*e/es;    % RH%
  