
function psi=psit(zet)
	x=(1-15*zet).^.5;
	psik=2*log((1+x)/2);
	x=(1-10.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);
    a=1.007;b=7;c=1.5;ac=a^3;x=zet(ii(jj)).^.333;
    psi(ii(jj))=log(1+zet(ii(jj))/ac)-3*b*x+(a^2*b+c)/2/a*log((a+x).^3./(ac+zet(ii(jj))))+sqrt(3)*(a^2*b-c)/a*(atan((2*x-a)/sqrt(3)/a)+pi/6);
end;

