function psi = psit_36(zet) % v 3.6
% for zet<0, instability
x = (1-15*zet).^0.5;
psik = 2*log((1+x)/2);
x = (1-34.15*zet).^0.3333;
psic = 1.5*log((1+x+x.*x)/3) - ...
       sqrt(3)*atan((1+2*x)/sqrt(3)) + ...
       4*atan(1)/sqrt(3);
f = zet.*zet./(1+zet.*zet);
psi = (1-f).*psik+f.*psic;
% for z>0, stability, Grachev 2007, Eq. 13
ii = find(zet>0);
if ~isempty(ii)
    ah = 5; bh = 5; ch = 3;
    BH = sqrt(ch^2 - 4);
    psi(ii) = -(bh/2)*log(1+ch*zet(ii)+zet(ii).^2) + ...
              (((bh*ch)/(2*BH)) - ...
              (ah/BH))*(log((2*zet(ii)+ch-BH)./(2*zet(ii)+ch+BH)) - ...
              log((ch-BH)./(ch+BH)));
end
end

