1      SUBROUTINE CBESY(Z, FNU, KODE, N, CY, NZ, CWRK, IERR)
2C***BEGIN PROLOGUE  CBESY
3C***DATE WRITTEN   830501   (YYMMDD)
4C***REVISION DATE  890801   (YYMMDD)
5C***CATEGORY NO.  B5K
6C***KEYWORDS  Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT,
7C             BESSEL FUNCTION OF SECOND KIND
8C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
9C***PURPOSE  TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT
10C***DESCRIPTION
11C
12C         ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
13C         BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE
14C         ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
15C         -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESY RETURNS THE SCALED
16C         FUNCTIONS
17C
18C         CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z)   I = 1,...,N , Y=AIMAG(Z)
19C
20C         WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
21C         LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
22C         ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
23C         (REF. 1).
24C
25C         INPUT
26C           Z      - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI
27C           FNU    - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0E0
28C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
29C                    KODE= 1  RETURNS
30C                             CY(I)=Y(FNU+I-1,Z), I=1,...,N
31C                        = 2  RETURNS
32C                             CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N
33C                             WHERE Y=AIMAG(Z)
34C           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
35C           CWRK   - A COMPLEX WORK VECTOR OF DIMENSION AT LEAST N
36C
37C         OUTPUT
38C           CY     - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
39C                    VALUES FOR THE SEQUENCE
40C                    CY(I)=Y(FNU+I-1,Z)  OR
41C                    CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y))  I=1,...,N
42C                    DEPENDING ON KODE.
43C           NZ     - NZ=0 , A NORMAL RETURN
44C                    NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO
45C                    UNDERFLOW (GENERALLY ON KODE=2)
46C           IERR   - ERROR FLAG
47C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
48C                    IERR=1, INPUT ERROR   - NO COMPUTATION
49C                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU+N-1 IS
50C                            TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
51C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
52C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
53C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
54C                            ACCURACY
55C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
56C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
57C                            CANCE BY ARGUMENT REDUCTION
58C                    IERR=5, ERROR              - NO COMPUTATION,
59C                            ALGORITHM TERMINATION CONDITION NOT MET
60C
61C***LONG DESCRIPTION
62C
63C         THE COMPUTATION IS CARRIED OUT BY THE FORMULA
64C
65C         Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I
66C
67C         WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z)
68C         AND H(2,FNU,Z) ARE CALCULATED IN CBESH.
69C
70C         FOR NEGATIVE ORDERS,THE FORMULA
71C
72C              Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
73C
74C         CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD
75C         INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE
76C         POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)*
77C         SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS
78C         NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A
79C         LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM
80C         CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS,
81C         WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF
82C         ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z).
83C
84C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
85C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
86C         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
87C         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
88C         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
89C         IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO
90C         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
91C         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
92C         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
93C         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
94C         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
95C         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
96C         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
97C         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
98C         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
99C         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
100C         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
101C         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
102C
103C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
104C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
105C         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
106C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
107C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
108C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
109C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
110C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
111C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
112C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
113C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
114C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
115C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
116C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
117C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
118C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
119C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
120C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
121C         OR -PI/2+P.
122C
123C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
124C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
125C                 COMMERCE, 1955.
126C
127C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
128C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
129C
130C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
131C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
132C
133C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
134C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
135C                 1018, MAY, 1985
136C
137C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
138C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
139C                 MATH. SOFTWARE, 1986
140C
141C***ROUTINES CALLED  CBESH,I1MACH,R1MACH
142C***END PROLOGUE  CBESY
143C
144      COMPLEX CWRK, CY, C1, C2, EX, HCI, Z, ZU, ZV
145      REAL ELIM, EY, FNU, R1, R2, TAY, XX, YY, R1MACH, ASCLE, RTOL,
146     * ATOL, AA, BB
147      INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH
148      DIMENSION CY(N), CWRK(N)
149C***FIRST EXECUTABLE STATEMENT  CBESY
150      XX = REAL(Z)
151      YY = AIMAG(Z)
152      IERR = 0
153      NZ=0
154      IF (XX.EQ.0.0E0 .AND. YY.EQ.0.0E0) IERR=1
155      IF (FNU.LT.0.0E0) IERR=1
156      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
157      IF (N.LT.1) IERR=1
158      IF (IERR.NE.0) RETURN
159      HCI = CMPLX(0.0E0,0.5E0)
160      CALL CBESH(Z, FNU, KODE, 1, N, CY, NZ1, IERR)
161      IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
162      CALL CBESH(Z, FNU, KODE, 2, N, CWRK, NZ2, IERR)
163      IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
164      NZ = MIN0(NZ1,NZ2)
165      IF (KODE.EQ.2) GO TO 60
166      DO 50 I=1,N
167        CY(I) = HCI*(CWRK(I)-CY(I))
168   50 CONTINUE
169      RETURN
170   60 CONTINUE
171      TOL = AMAX1(R1MACH(4),1.0E-18)
172      K1 = I1MACH(12)
173      K2 = I1MACH(13)
174      K = MIN0(IABS(K1),IABS(K2))
175      R1M5 = R1MACH(5)
176C-----------------------------------------------------------------------
177C     ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
178C-----------------------------------------------------------------------
179      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
180      R1 = COS(XX)
181      R2 = SIN(XX)
182      EX = CMPLX(R1,R2)
183      EY = 0.0E0
184      TAY = ABS(YY+YY)
185      IF (TAY.LT.ELIM) EY = EXP(-TAY)
186      IF (YY.LT.0.0E0) GO TO 90
187      C1 = EX*CMPLX(EY,0.0E0)
188      C2 = CONJG(EX)
189   70 CONTINUE
190      NZ = 0
191      RTOL = 1.0E0/TOL
192      ASCLE = R1MACH(1)*RTOL*1.0E+3
193      DO 80 I=1,N
194C       CY(I) = HCI*(C2*CWRK(I)-C1*CY(I))
195        ZV = CWRK(I)
196        AA=REAL(ZV)
197        BB=AIMAG(ZV)
198        ATOL=1.0E0
199        IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 75
200          ZV = ZV*CMPLX(RTOL,0.0E0)
201          ATOL = TOL
202   75   CONTINUE
203        ZV = ZV*C2*HCI
204        ZV = ZV*CMPLX(ATOL,0.0E0)
205        ZU=CY(I)
206        AA=REAL(ZU)
207        BB=AIMAG(ZU)
208        ATOL=1.0E0
209        IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 85
210          ZU = ZU*CMPLX(RTOL,0.0E0)
211          ATOL = TOL
212   85   CONTINUE
213        ZU = ZU*C1*HCI
214        ZU = ZU*CMPLX(ATOL,0.0E0)
215        CY(I) = ZV - ZU
216        IF (CY(I).EQ.CMPLX(0.0E0,0.0E0) .AND. EY.EQ.0.0E0) NZ = NZ + 1
217   80 CONTINUE
218      RETURN
219   90 CONTINUE
220      C1 = EX
221      C2 = CONJG(EX)*CMPLX(EY,0.0E0)
222      GO TO 70
223  170 CONTINUE
224      NZ = 0
225      RETURN
226      END
227