GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
xzsqrt.f
Go to the documentation of this file.
1  SUBROUTINE xzsqrt(AR, AI, BR, BI)
2 C***BEGIN PROLOGUE XZSQRT
3 C***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZBESY,ZAIRY,ZBIRY
4 C
5 C DOUBLE PRECISION COMPLEX SQUARE ROOT, B=CSQRT(A)
6 C
7 C***ROUTINES CALLED XZABS
8 C***END PROLOGUE XZSQRT
9  DOUBLE PRECISION ar, ai, br, bi, zm, dtheta, dpi, drt
10  DOUBLE PRECISION xzabs
11  DATA drt , dpi / 7.071067811865475244008443621d-1,
12  1 3.141592653589793238462643383d+0/
13  zm = xzabs(ar,ai)
14  zm = dsqrt(zm)
15  IF (ar.EQ.0.0d+0) go to 10
16  IF (ai.EQ.0.0d+0) go to 20
17  dtheta = datan(ai/ar)
18  IF (dtheta.LE.0.0d+0) go to 40
19  IF (ar.LT.0.0d+0) dtheta = dtheta - dpi
20  go to 50
21  10 IF (ai.GT.0.0d+0) go to 60
22  IF (ai.LT.0.0d+0) go to 70
23  br = 0.0d+0
24  bi = 0.0d+0
25  RETURN
26  20 IF (ar.GT.0.0d+0) go to 30
27  br = 0.0d+0
28  bi = dsqrt(dabs(ar))
29  RETURN
30  30 br = dsqrt(ar)
31  bi = 0.0d+0
32  RETURN
33  40 IF (ar.LT.0.0d+0) dtheta = dtheta + dpi
34  50 dtheta = dtheta*0.5d+0
35  br = zm*dcos(dtheta)
36  bi = zm*dsin(dtheta)
37  RETURN
38  60 br = zm*drt
39  bi = zm*drt
40  RETURN
41  70 br = zm*drt
42  bi = -zm*drt
43  RETURN
44  END