
function psi=psit(zet)
	x=(1-15*zet).^.5;
	psik=2*log((1+x)/2);
	x=(1-34.15*zet).^.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;  
   
   ii=find(zet>0);
if ~isempty(ii);
	%psi=-4.7*zet;
    psi(ii)=-5*zet(ii);
    jj=find(zet(ii)>.03);
    b=22;c=-2.5;ac=7;a=ac^0.333;
    x=zet(ii(jj)).^.333;
    psi(ii(jj))=log(1+zet(ii(jj))/ac)+(a*b-c)/2/a^2*log((a+x).^3./(ac+zet(ii(jj))))-sqrt(3)*(a*b+c)/a^2*(atan((2*x-a)/sqrt(3)/a)+pi/6);
end;

