function psi = psiu_36(zet) % v 3.6
% for zet<0, instability
x = (1-15*zet).^0.25;
psik = 2*log((1+x)/2) + log((1+x.*x)/2) - 2*atan(x) + 2*atan(1);
x = (1-10.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. 12
ii = find(zet>0);
if ~isempty(ii)
    am = 5; bm = am/6.5;
    BM = ((1-bm)/bm)^(1/3);
    y = (1+zet(ii)).^(1/3);
    psi(ii) = -(3*am/bm).*(y-1) + ...
              ((am*BM)./(2*bm))*(2*log((BM+y)./(BM+1)) -   ...
              log((BM^2 - BM*y+y.^2)./(BM^2 - BM+1)) + ...
              2*sqrt(3)*atan((2*y - BM)/(BM*sqrt(3))) - ...
              2*sqrt(3)*atan((2 - BM)/(BM*sqrt(3))));
end
end
