GNU Octave: libcruft/amos/xzsqrt.f Source File

Go to the documentation of this file.
00001       SUBROUTINE XZSQRT(AR, AI, BR, BI)
00002 C***BEGIN PROLOGUE  XZSQRT
00003 C***REFER TO  ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY
00004 C
00005 C     DOUBLE PRECISION COMPLEX SQUARE ROOT, B=CSQRT(A)
00006 C
00007 C***ROUTINES CALLED  XZABS
00008 C***END PROLOGUE  XZSQRT
00009       DOUBLE PRECISION AR, AI, BR, BI, ZM, DTHETA, DPI, DRT
00010       DOUBLE PRECISION XZABS
00011       DATA DRT , DPI / 7.071067811865475244008443621D-1,
00012      1                 3.141592653589793238462643383D+0/
00013       ZM = XZABS(AR,AI)
00014       ZM = DSQRT(ZM)
00015       IF (AR.EQ.0.0D+0) GO TO 10
00016       IF (AI.EQ.0.0D+0) GO TO 20
00017       DTHETA = DATAN(AI/AR)
00018       IF (DTHETA.LE.0.0D+0) GO TO 40
00019       IF (AR.LT.0.0D+0) DTHETA = DTHETA - DPI
00020       GO TO 50
00021    10 IF (AI.GT.0.0D+0) GO TO 60
00022       IF (AI.LT.0.0D+0) GO TO 70
00023       BR = 0.0D+0
00024       BI = 0.0D+0
00025       RETURN
00026    20 IF (AR.GT.0.0D+0) GO TO 30
00027       BR = 0.0D+0
00028       BI = DSQRT(DABS(AR))
00029       RETURN
00030    30 BR = DSQRT(AR)
00031       BI = 0.0D+0
00032       RETURN
00033    40 IF (AR.LT.0.0D+0) DTHETA = DTHETA + DPI
00034    50 DTHETA = DTHETA*0.5D+0
00035       BR = ZM*DCOS(DTHETA)
00036       BI = ZM*DSIN(DTHETA)
00037       RETURN
00038    60 BR = ZM*DRT
00039       BI = ZM*DRT
00040       RETURN
00041    70 BR = ZM*DRT
00042       BI = -ZM*DRT
00043       RETURN
00044       END