[xr,xi,yr,yi]=csroot(xr,xi,yr,yi);

function [xr,xi,yr,yi]=csroot(xr,xi,yr,yi);
%***BEGIN PROLOGUE  CSROOT
%***SUBSIDIARY
%***PURPOSE  Compute the complex square root of a complex number.
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (CSROOT-S)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     (YR,YI) = complex sqrt(XR,XI)
%
%***SEE ALSO  EISDOC
%***ROUTINES CALLED  PYTHAG
%***REVISION HISTORY  (YYMMDD)
%   811101  DATE WRITTEN
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  CSROOT
persistent s ti tr ; 

if isempty(s), s=0; end;
if isempty(tr), tr=0; end;
if isempty(ti), ti=0; end;
%
%     BRANCH CHOSEN SO THAT YR .GE. 0.0 AND SIGN(YI) .EQ. SIGN(XI)
%***FIRST EXECUTABLE STATEMENT  CSROOT
tr = xr;
ti = xi;
s = sqrt(0.5e0.*(pythag(tr,ti)+abs(tr)));
if( tr>=0.0e0 )
yr = s;
end;
if( ti<0.0e0 )
s = -s;
end;
if( tr<=0.0e0 )
yi = s;
end;
if( tr<0.0e0 )
yr = 0.5e0.*(ti./yi);
end;
if( tr>0.0e0 )
yi = 0.5e0.*(ti./yr);
end;
end
%DECK CSROT