xzsqrt.f

Go to the documentation of this file.

      SUBROUTINE xzsqrt(AR, AI, BR, BI)
C***BEGIN PROLOGUE  XZSQRT
C***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY
C
C     DOUBLE PRECISION COMPLEX SQUARE ROOT, B=CSQRT(A)
C
C***ROUTINES CALLED  XZABS
C***END PROLOGUE  XZSQRT
      DOUBLE PRECISION ar, ai, br, bi, zm, dtheta, dpi, drt
      DOUBLE PRECISION xzabs
      DATA drt , dpi / 7.071067811865475244008443621d-1,
     1                 3.141592653589793238462643383d+0/
      zm = xzabs(ar,ai)
      zm = dsqrt(zm)
      IF (ar.EQ.0.0d+0) go to 10
      IF (ai.EQ.0.0d+0) go to 20
      dtheta = datan(ai/ar)
      IF (dtheta.LE.0.0d+0) go to 40
      IF (ar.LT.0.0d+0) dtheta = dtheta - dpi
      go to 50
   10 IF (ai.GT.0.0d+0) go to 60
      IF (ai.LT.0.0d+0) go to 70
      br = 0.0d+0
      bi = 0.0d+0
      RETURN
   20 IF (ar.GT.0.0d+0) go to 30
      br = 0.0d+0
      bi = dsqrt(dabs(ar))
      RETURN
   30 br = dsqrt(ar)
      bi = 0.0d+0
      RETURN
   40 IF (ar.LT.0.0d+0) dtheta = dtheta + dpi
   50 dtheta = dtheta*0.5d+0
      br = zm*dcos(dtheta)
      bi = zm*dsin(dtheta)
      RETURN
   60 br = zm*drt
      bi = zm*drt
      RETURN
   70 br = zm*drt
      bi = -zm*drt
      RETURN
      END