complex.f90

! ---------------------------------------------------------------------
!      Utility subroutines used by any program from Numath library
!      with (not intrinsic) complex numbers z = (zr,zi).
! ---------------------------------------------------------------------
! Reference: From Numath Library By Tuan Dang Trong in Fortran 77
!            [BIBLI 18].
!
!                               F90 Release 1.0 By J-P Moreau, Paris
!                                         (www.jpmoreau.fr)
! ---------------------------------------------------------------------
MODULE complex

CONTAINS

SUBROUTINE zsqrt(AR, AI, BR, BI)
!***BEGIN PROLOGUE  ZSQRT
!***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY
!
!     DOUBLE PRECISION COMPLEX SQUARE ROOT, B=CSQRT(A)
!
!***ROUTINES CALLED  ZABS
!***END PROLOGUE  ZSQRT
      DOUBLE PRECISION ar, ai, br, bi, zm, dtheta, dpi, drt
      DATA drt , dpi / 7.071067811865475244008443621d-1, &
                       3.141592653589793238462643383d+0/
      zm = zabs(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 SUBROUTINE zsqrt

SUBROUTINE zexp(AR, AI, BR, BI)
!***BEGIN PROLOGUE  ZEXP
!***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY
!
!     DOUBLE PRECISION COMPLEX EXPONENTIAL FUNCTION B=EXP(A)
!
!***ROUTINES CALLED  (NONE)
!***END PROLOGUE  ZEXP
      DOUBLE PRECISION ar, ai, br, bi, zm, ca, cb
      zm = dexp(ar)
      ca = zm*dcos(ai)
      cb = zm*dsin(ai)
      br = ca
      bi = cb
RETURN
END SUBROUTINE zexp

SUBROUTINE zmlt(AR, AI, BR, BI, CR, CI)
!***BEGIN PROLOGUE  ZMLT
!***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY
!
!     DOUBLE PRECISION COMPLEX MULTIPLY, C=A*B.
!
!***ROUTINES CALLED  (NONE)
!***END PROLOGUE  ZMLT
DOUBLE PRECISION ar, ai, br, bi, cr, ci, ca, cb
      ca = ar*br - ai*bi
      cb = ar*bi + ai*br
      cr = ca
      ci = cb
RETURN
END SUBROUTINE zmlt

SUBROUTINE zdiv(AR, AI, BR, BI, CR, CI)
!***BEGIN PROLOGUE  ZDIV
!***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY
!
!     DOUBLE PRECISION COMPLEX DIVIDE C=A/B.

!***ROUTINES CALLED  ZABS
!***END PROLOGUE  ZDIV
      DOUBLE PRECISION ar, ai, br, bi, cr, ci, bm, ca, cb, cc, cd
      bm = 1.0d0/zabs(br,bi)
      cc = br*bm
      cd = bi*bm
      ca = (ar*cc+ai*cd)*bm
      cb = (ai*cc-ar*cd)*bm
      cr = ca
      ci = cb
RETURN
END SUBROUTINE zdiv

SUBROUTINE zlog(AR, AI, BR, BI, IERR)
!***BEGIN PROLOGUE  ZLOG
!***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY

!     DOUBLE PRECISION COMPLEX LOGARITHM B=CLOG(A)
!     IERR=0,NORMAL RETURN      IERR=1, Z=CMPLX(0.0,0.0)
!***ROUTINES CALLED  ZABS
!***END PROLOGUE  ZLOG
      DOUBLE PRECISION ar, ai, br, bi, zm, dtheta, dpi, dhpi
      DATA dpi , dhpi  / 3.141592653589793238462643383d+0, &
                         1.570796326794896619231321696d+0/

      ierr=0
      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.EQ.0.0d+0) GO TO 60
      bi = dhpi
      br = dlog(dabs(ai))
      IF (ai.LT.0.0d+0) bi = -bi
      RETURN
   20 IF (ar.GT.0.0d+0) GO TO 30
      br = dlog(dabs(ar))
      bi = dpi
      RETURN
   30 br = dlog(ar)
      bi = 0.0d+0
      RETURN
   40 IF (ar.LT.0.0d+0) dtheta = dtheta + dpi
   50 zm = zabs(ar,ai)
      br = dlog(zm)
      bi = dtheta
      RETURN
   60 CONTINUE
      ierr=1
RETURN
END SUBROUTINE zlog

DOUBLE PRECISION FUNCTION zabs(ZR, ZI)
!***BEGIN PROLOGUE  ZABS
!***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY

!     ZABS COMPUTES THE ABSOLUTE VALUE OR MAGNITUDE OF A DOUBLE
!     PRECISION COMPLEX VARIABLE CMPLX(ZR,ZI)

!***ROUTINES CALLED  (NONE)
!***END PROLOGUE  ZABS
      DOUBLE PRECISION zr, zi, u, v, q, s
      u = dabs(zr)
      v = dabs(zi)
      s = u + v
!-----------------------------------------------------------------------
!     S*1.0D0 MAKES AN UNNORMALIZED UNDERFLOW ON CD! MACHINES INTO A
!     TRUE FLOATING ZERO
!-----------------------------------------------------------------------
      s = s*1.0d+0
      IF (s.EQ.0.0d+0) GO TO 20
      IF (u.GT.v) GO TO 10
      q = u/v
      zabs = v*dsqrt(1.d+0+q*q)
      RETURN
   10 q = v/u
      zabs = u*dsqrt(1.d+0+q*q)
      RETURN
   20 zabs = 0.0d+0
      RETURN
END FUNCTION zabs

END MODULE complex
! end of file Complex.f90