FUNCTION CABS(X)
C RETURNS ABS(X), WHERE X IS COMPLEX
      DIMENSION X(2)
      CABS=SQRT(X(1)*X(1)+X(2)*X(2))
      RETURN
      END
C   
      SUBROUTINE CADD(X,Y,Z)
C RETURNS Z=X+Y, WHERE X, Y, AND Z ARE COMPLEX
      DIMENSION X(2),Y(2),Z(2)
      Z(1)=X(1)+Y(1)
      Z(2)=X(2)+Y(2)
      RETURN
      END
C   
      SUBROUTINE CSUB(X,Y,Z)
C RETURNS Z=X-Y, WHERE X, Y, AND Z ARE COMPLEX
      DIMENSION X(2),Y(2),Z(2)
      Z(1)=X(1)-Y(1)
      Z(2)=X(2)-Y(2)
      RETURN
      END
C   
      SUBROUTINE CMUL(X,Y,Z)
C RETURNS Z=X*Y, WHERE X, Y, AND Z ARE COMPLEX
      DIMENSION X(2),Y(2),Z(2)
      Z(1)=X(1)*Y(1)-X(2)*Y(2)
      Z(2)=X(1)*Y(2)+X(2)*Y(1)
      RETURN
      END
C   
      SUBROUTINE CDIV(X,Y,Z)
C RETURNS Z=X/Y, WHERE X, Y, AND Z ARE COMPLEX
      DIMENSION X(2),Y(2),Z(2)
      S=1./(Y(1)*Y(1)+Y(2)*Y(2))
      Z(1)=S*(X(1)*Y(1)+X(2)*Y(2))
      Z(2)=S*(X(2)*Y(1)-X(1)*Y(2))
      RETURN
      END
C   
      SUBROUTINE CSQRT(X,Z)
C RETURNS Z=SQRT(X), WHERE X AND Z ARE COMPLEX
      DIMENSION X(2),Z(2)
      CALL POLAR(X,R,THETA)
      CALL RECT(Z,SQRT(R),THETA/2.)
      RETURN
      END
C   
      SUBROUTINE POLAR(Z,R,THETA)
C RETURNS POLAR COORDINATES R AND THETA OF A COMPLEX NUMBER Z
      DIMENSION Z(2)
      DATA PI/3.1415927/
      THETA=0.
      R=CABS(Z)
      IF (R.EQ.0.) RETURN
      IF (ABS(Z(1)).LT.ABS(Z(2))) GOTO 10
      THETA=ATAN(Z(2)/Z(1))
      IF (Z(1).LT.0.) THETA=THETA+PI
      GOTO 20
   10 THETA=-ATAN(Z(1)/Z(2))+PI/2.
      IF (Z(2).LT.0.) THETA=THETA-PI
   20 IF (THETA.LT.0.) THETA=THETA+2.*PI
      RETURN
      END
C   
      SUBROUTINE RECT(Z,R,THETA)
C RETURNS THE COMPLEX NUMBER Z WHOSE POLAR COORDINATES
C ARE R AND THETA
      DIMENSION Z(2)
      Z(1)=R*COS(THETA)
      Z(2)=R*SIN(THETA)
      RETURN
      END