function [u] = csrt(z)
% The following function returns a complex square root with a
% positive imaginary part.
u=sqrt(z);
J=find(imag(u)<0.0);
u(J)=(-1.0)*u(J);