      SUBROUTINE A02AAF(XXR, XXI, YR, YI)
C     
C     COMPUTES THE SQUARE ROOT OF A COMPLEX NUMBER
C
      DOUBLE PRECISION ZERO, HALF, XI, XXI, H, XXR, YR, YI, A02ABF
      DATA ZERO /0.0D0/, HALF /0.5D0/
C
      XI = XXI
      H = DSQRT((DABS(XXR)+A02ABF(XXR,XI))*HALF)
      IF (XI.NE.ZERO) XI = XI/(H+H)
      IF (XXR.LT.ZERO) GO TO 20
      YR = H
      YI = XI
      RETURN
   20 IF (XI.LT.ZERO) GO TO 40
      YR = XI
      YI = H
      RETURN
   40 YR = -XI
      YI = -H
      RETURN
      END

      DOUBLE PRECISION FUNCTION A02ABF(XXR, XXI)
C
C     RETURNS THE ABSOLUTE VALUE OF A COMPLEX NUMBER VIA ROUTINE
C
      DOUBLE PRECISION ZERO, ONE, XR, XXR, XI, XXI, H
      DATA ZERO /0.0D0/, ONE /1.0D0/
C
      XR = DABS(XXR)
      XI = DABS(XXI)
      IF (XI.LE.XR) GO TO 20
      H = XR
      XR = XI
      XI = H
   20 IF (XI.NE.ZERO) GO TO 40
      A02ABF = XR
      RETURN
   40 H = XR*DSQRT(ONE+(XI/XR)**2)
      A02ABF = H
      RETURN
      END

      DOUBLE PRECISION FUNCTION A02ABS(XX)
C           SUM OF ABSOLUTE VALUES
      COMPLEX*16 XX
      DOUBLE PRECISION XXR,XXI
      CALL A02CD(XX,XXR,XXI)
      A02ABS = DABS(XXR) + DABS(XXI)
      RETURN
      END


      SUBROUTINE A02ACF(XXR, XXI, YYR, YYI, ZR, ZI)
C
C     DIVIDES ONE COMPLEX NUMBER BY A SECOND  Z = X / Y
C
      DOUBLE PRECISION ONE, YYR, YYI, H, A, ZR, XXR, XXI, ZI
      DATA ONE /1.0D0/
C
      IF (DABS(YYR).LE.DABS(YYI)) GO TO 20
      H = YYI/YYR
      A = ONE/(H*YYI+YYR)
      ZR = (XXR+H*XXI)*A
      ZI = (XXI-H*XXR)*A
      RETURN
   20 H = YYR/YYI
      A = ONE/(H*YYR+YYI)
      ZR = (H*XXR+XXI)*A
      ZI = (H*XXI-XXR)*A
      RETURN
      END

      SUBROUTINE A02DC(XR,XI,X)
      DOUBLE PRECISION XR,XI
      COMPLEX*16 X
      X = CMPLX(XR, XI)
      RETURN
      END

      SUBROUTINE A02CD(X,XR,XI)
      DOUBLE PRECISION XR,XI
      COMPLEX*16 X
      XR = REAL(X)
      XI = AIMAG(X)
      RETURN
      END

