function H2O = H2O_psych(dbC, dpC, stPa)
% H2O_psych.m edition 30th January 2017 given drybulb, dewpoint, and pressure;
% example dbC=25; dpC=20; stPa=101325;
%% psychrometrics from ASHRAE Fundamentals 1993
C(1) = -5.6745359E+03;
% C(2) = -5.1523058E-01;
C(2) = +6.3925247;
C(3) = -9.6778430E-03;
C(4) = +6.2215701E-07;
C(5) = +2.0747825E-09;
C(6) = -9.4840240E-13;
C(7) = +4.1635019;
C(8) = -5.8002206E+03;
% C(9) = +5.5162560;
C(9) = 1.3914993;
C(10)= -4.8640239E-02;
C(11)= +4.1764768E-05;
C(12)= -1.4452093E-08;
C(13)= +6.5459673;
C(14)= 6.54;
C(15)=14.526;
C(16)= 0.7389;
C(17)= 0.09486;
C(18)= 0.04569;

%% 24th March 2016 amendment to avoid supersaturation
if dbC<dpC;     dpC=dbC; end
%% applies to MAX, mean, and MIN (constant dewpoint)
% ASHRAE 1993 Fundamentals Chapter 6 numerical calculation
% given drybulb (dbC), dewpoint (dpC), and pressure (kPa)
% water vapour partial pressure pw = pws(dpC) by equations 3,4
% then humidity ratio W from equation 20
DPA=dpC +273.15; % absolute temperature Kelvins of dew point condition
if dpC<0 % below freezing
    lnpwsdp = C(1)/DPA ...
        + C(2) ...
        + C(3)*DPA ...
        + C(4)*DPA^2 ...
        + C(5)*DPA^3 ...
        + C(6)*DPA^4 ...
        + C(7)*log(DPA);
    pwsxdp=exp(lnpwsdp); % Pascals vapour pressure at saturation with DP
else % above freezing
    lnpwsdp = C(8)/DPA ...
        + C(9) ...
        + C(10)*DPA ...
        + C(11)*DPA^2 ...
        + C(12)*DPA^3 ...
        + C(13)*log(DPA);
    pwsxdp=exp(lnpwsdp); % Pascals vapour pressure at saturation with DP
end % now we have vapour pressure at DP saturation
H2O=0.62198*pwsxdp/(stPa-pwsxdp); % equation 20
H2O=1000*H2O; % Humidity expressed in g H20 per kg dry air
