function [wbC H20]= psych(dbC, dpC, stPa)
% psych.m edition 30th January 2017 given drybulb, dewpoint, and pressure;
% example dbC=25; dpC=20; stPa=101325; [wbC H20]= psych(dbC, dpC, stPa)
%% 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;

%% initial estimate of webbulb is needed for iteration to convergence
% Ben Liley's formula developed after Australian BoM
% temperature dbC in degrees Celsius
SatVP_a = 6108 * exp(17.2694 * dbC/(dbC + 237.3)); % Pascals
SatVP_b =  61064 * exp(21.88 * dbC/(dbC + 265.5)) ; % vapour pressure over ice correction for dbC < 0
SatVP =min(SatVP_a,SatVP_b); % minimum of the two
dP = 1013.25 -stPa/100; % hPa pressure deficit from standard pressure
wbC = 0.07 + (0.34902 + 1.62e-4*dP) * dpC ...
    + (0.63730 - 1.88e-4*dP) * dbC ...
    + (8.028e-3 + 2.35e-6*dP) * dpC^2 ...
    - (4.176e-3 + 9.3e-7*dP) * dbC*dpC ...
    - 4.468e-3 * dbC^2 ...
    + (6.035e-5 + 1.34e-8*dP + 3.9e-10*dP^2) * dpC^3 ...
    - (5.522e-5 + 2.6e-8*dP) * dbC*dpC^2 ...
    + 1.32e-5 * dpC*dbC^2 ...
    + 1.9e-5 * dbC^3;
wbC=min(wbC,dbC); % Liley's wetbulb estimate degrees C
% ASHRAE 1993 Fundamentals Chapter 6 numerical calculation
% given drybulb (dbC), wetbulb (wbC), and pressure (kPa)
% water vapour partial pressure pw* at wetbulb by equations 3,4
% then humidity ratio W* at saturation from equation 21
% then humidity ratio W from equation 33

WBA=wbC +273.15; % absolute temperature Kelvins of wetbulb condition
if dpC<0 % corresponding drybulb below freezing point
    lnpwswb = C(1)/WBA ...
        + C(2) ...
        + C(3)*WBA ...
        + C(4)*WBA^2 ...
        + C(5)*WBA^3 ...
        + C(6)*WBA^4 ...
        + C(7)*log(WBA);
    pwsxdp=exp(lnpwswb); % Pascals vapour pressure at saturation with DP
else % above freezing
    lnpwswb = C(8)/WBA ...
        + C(9) ...
        + C(10)*WBA ...
        + C(11)*WBA^2 ...
        + C(12)*WBA^3 ...
        + C(13)*log(WBA);
    pwsxdp=exp(lnpwswb); % Pascals vapour pressure at saturation with DP
end % now we have vapour pressure at DP saturation
H2O_star=0.62198*pwsxdp/(stPa-pwsxdp); % equation 21
H2O_liley=((2501-2.326*wbC)*H2O_star-1.006*(dbC-wbC))/(2501+1.805*dbC-4.186*wbC); % equation 33

H2O_liley=1000*H2O_liley'; % vertical array of Humidity expressed in g H20 per kg dry air
wbLiley=wbC; % backup copy of Liley's estimate of wetbulb tmeperature

%% 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


%%
% comparison with initial wetbulb estimate

discrep=0; % zero assuming no discrepencies between initial estimate of WBT and that calculated below
WB(1)=wbC; % initial guess of WB is taken from wbC %%%
go=1; % go loop flag =0 to skip
limit=1000; % limit on kounter kr (don't need to increase this much more than maximum of "enough"
convergence=1/1000000; % convergence criteria (one part per million)
for kr=1:limit; % kounter to refine estimate of WB 
    if go==1; % go through all these calcs
        WBA=WB(kr) +273.15; % absolute temperature Kelvins of wet bulb condition
        if WB(kr)<0 % below freezing
            lnpwsx = C(1)/WBA + C(2) + C(3)*WBA + C(4)*WBA^2 + C(5)*WBA^3 + C(6) * WBA^4 + C(7)*log(WBA);
            pwsx=exp(lnpwsx); % Pascals vapour pressure at saturation with WB
        else % above freezing
            lnpwsx = C(8)/WBA + C(9) + C(10)*WBA + C(11)*WBA^2 + C(12)*WBA^3 + C(13)*log(WBA);
            pwsx=exp(lnpwsx); % Pascals vapour pressure at saturation with WB
        end % now we have vapour pressure at WB saturation
        Wsatx=0.62198*pwsx/(stPa-pwsx); % g-H2O/g-air at saturation of WB
        W = ((2501-2.326*WB(kr))*Wsatx-1.006*(dbC-WB(kr)))/(2501+1.86*dbC-4.186*WB(kr)); 
        % estimated W (g-H2O/g-air)
        % error(kr)=((W*1000)-H2O)/(W*1000)      % delta:ashrae
        H2Oerror(kr)=(1-H2O/(W*1000)); % overestimate of H20 on psychrometric y-axis
        H2Odiff(kr)=W*1000-H2O_liley; % difference of H20 on psychrometric y-axis
        WBerror(kr) = WB(kr)- wbLiley; % wetbulb error = WBerror(kr) = WB(kr)- wbC where the latter is Liley's estimate
        if kr<limit; % next estimate of wet bulb temperature
            WB(kr+1)=WB(kr)-H2Oerror(kr)/log10((3162/kr)^2); % scaled by product of error and degree adjustment
        end;  % there you have next interation guess for WB
        if abs(H2Oerror(kr))< convergence;
           % if abs(WBerror(kr))<0.01 ; % previously abs(H2Oerror(kr))< convergence;
           enough=kr; % needed interations for the day j
           % substitute wbC the latest estimate of wetbulb, overwriting liley's estimate
           wbC=WB(kr);
           go=0; % skip out of calculations after convergence satisfied
        else
           if kr<limit;
           else % at the limit then accept enough as the limit
               enough=kr;
           end
        end;
    end;
end
errWB=WBerror(enough); % error of wetbulb 
H20=W*1000; % absolute humidity (g/kg)
errH2O=H2Oerror(enough); % error of absolute humidity
