1      SUBROUTINE ZBESH(ZR, ZI, FNU, KODE, M, N, CYR, CYI, NZ, IERR)
2C***BEGIN PROLOGUE  ZBESH
3C***DATE WRITTEN   830501   (YYMMDD)
4C***REVISION DATE  890801   (YYMMDD)
5C***CATEGORY NO.  B5K
6C***KEYWORDS  H-BESSEL FUNCTIONS,BESSEL FUNCTIONS OF COMPLEX ARGUMENT,
7C             BESSEL FUNCTIONS OF THIRD KIND,HANKEL FUNCTIONS
8C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
9C***PURPOSE  TO COMPUTE THE H-BESSEL FUNCTIONS OF A COMPLEX ARGUMENT
10C***DESCRIPTION
11C
12C                      ***A DOUBLE PRECISION ROUTINE***
13C         ON KODE=1, ZBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
14C         HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1
15C         OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX
16C         Z.NE.CMPLX(0.0,0.0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI.
17C         ON KODE=2, ZBESH RETURNS THE SCALED HANKEL FUNCTIONS
18C
19C         CY(I)=EXP(-MM*Z*I)*H(M,FNU+J-1,Z)       MM=3-2*M,   I**2=-1.
20C
21C         WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER AND
22C         LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN THE
23C         NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1).
24C
25C         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
26C           ZR,ZI  - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0),
27C                    -PT.LT.ARG(Z).LE.PI
28C           FNU    - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0D0
29C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
30C                    KODE= 1  RETURNS
31C                             CY(J)=H(M,FNU+J-1,Z),   J=1,...,N
32C                        = 2  RETURNS
33C                             CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M))
34C                                  J=1,...,N  ,  I**2=-1
35C           M      - KIND OF HANKEL FUNCTION, M=1 OR 2
36C           N      - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1
37C
38C         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
39C           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
40C                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
41C                    CY(J)=H(M,FNU+J-1,Z)  OR
42C                    CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M))  J=1,...,N
43C                    DEPENDING ON KODE, I**2=-1.
44C           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
45C                    NZ= 0   , NORMAL RETURN
46C                    NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE
47C                              TO UNDERFLOW, CY(J)=CMPLX(0.0D0,0.0D0)
48C                              J=1,...,NZ WHEN Y.GT.0.0 AND M=1 OR
49C                              Y.LT.0.0 AND M=2. FOR THE COMPLMENTARY
50C                              HALF PLANES, NZ STATES ONLY THE NUMBER
51C                              OF UNDERFLOWS.
52C           IERR   - ERROR FLAG
53C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
54C                    IERR=1, INPUT ERROR   - NO COMPUTATION
55C                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU TOO
56C                            LARGE OR CABS(Z) TOO SMALL OR BOTH
57C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
58C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
59C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
60C                            ACCURACY
61C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
62C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
63C                            CANCE BY ARGUMENT REDUCTION
64C                    IERR=5, ERROR              - NO COMPUTATION,
65C                            ALGORITHM TERMINATION CONDITION NOT MET
66C
67C***LONG DESCRIPTION
68C
69C         THE COMPUTATION IS CARRIED OUT BY THE RELATION
70C
71C         H(M,FNU,Z)=(1/MP)*EXP(-MP*FNU)*K(FNU,Z*EXP(-MP))
72C             MP=MM*HPI*I,  MM=3-2*M,  HPI=PI/2,  I**2=-1
73C
74C         FOR M=1 OR 2 WHERE THE K BESSEL FUNCTION IS COMPUTED FOR THE
75C         RIGHT HALF PLANE RE(Z).GE.0.0. THE K FUNCTION IS CONTINUED
76C         TO THE LEFT HALF PLANE BY THE RELATION
77C
78C         K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z)
79C         MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1
80C
81C         WHERE I(FNU,Z) IS THE I BESSEL FUNCTION.
82C
83C         EXPONENTIAL DECAY OF H(M,FNU,Z) OCCURS IN THE UPPER HALF Z
84C         PLANE FOR M=1 AND THE LOWER HALF Z PLANE FOR M=2.  EXPONENTIAL
85C         GROWTH OCCURS IN THE COMPLEMENTARY HALF PLANES.  SCALING
86C         BY EXP(-MM*Z*I) REMOVES THE EXPONENTIAL BEHAVIOR IN THE
87C         WHOLE Z PLANE FOR Z TO INFINITY.
88C
89C         FOR NEGATIVE ORDERS,THE FORMULAE
90C
91C               H(1,-FNU,Z) = H(1,FNU,Z)*CEXP( PI*FNU*I)
92C               H(2,-FNU,Z) = H(2,FNU,Z)*CEXP(-PI*FNU*I)
93C                         I**2=-1
94C
95C         CAN BE USED.
96C
97C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
98C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
99C         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
100C         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
101C         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
102C         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
103C         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
104C         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
105C         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
106C         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
107C         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
108C         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
109C         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
110C         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
111C         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
112C         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
113C         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
114C         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
115C         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
116C
117C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
118C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
119C         ROUNDOFF,1.0D-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
120C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
121C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
122C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
123C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
124C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
125C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
126C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
127C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
128C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
129C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
130C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
131C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
132C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
133C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
134C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
135C         OR -PI/2+P.
136C
137C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
138C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
139C                 COMMERCE, 1955.
140C
141C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
142C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
143C
144C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
145C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
146C
147C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
148C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
149C                 1018, MAY, 1985
150C
151C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
152C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
153C                 MATH. SOFTWARE, 1986
154C
155C***ROUTINES CALLED  ZACON,ZBKNU,ZBUNK,ZUOIK,ZABS,I1MACH,D1MACH
156C***END PROLOGUE  ZBESH
157C
158C     COMPLEX CY,Z,ZN,ZT,CSGN
159      DOUBLE PRECISION AA, ALIM, ALN, ARG, AZ, CYI, CYR, DIG, ELIM,
160     * FMM, FN, FNU, FNUL, HPI, RHPI, RL, R1M5, SGN, STR, TOL, UFL, ZI,
161     * ZNI, ZNR, ZR, ZTI, D1MACH, ZABS, BB, ASCLE, RTOL, ATOL, STI,
162     * CSGNR, CSGNI
163      INTEGER I, IERR, INU, INUH, IR, K, KODE, K1, K2, M,
164     * MM, MR, N, NN, NUF, NW, NZ, I1MACH
165      DIMENSION CYR(N), CYI(N)
166C
167      DATA HPI /1.57079632679489662D0/
168C
169C***FIRST EXECUTABLE STATEMENT  ZBESH
170      IERR = 0
171      NZ=0
172      IF (ZR.EQ.0.0D0 .AND. ZI.EQ.0.0D0) IERR=1
173      IF (FNU.LT.0.0D0) IERR=1
174      IF (M.LT.1 .OR. M.GT.2) IERR=1
175      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
176      IF (N.LT.1) IERR=1
177      IF (IERR.NE.0) RETURN
178      NN = N
179C-----------------------------------------------------------------------
180C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
181C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
182C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
183C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
184C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
185C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
186C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
187C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
188C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
189C-----------------------------------------------------------------------
190      TOL = DMAX1(D1MACH(4),1.0D-18)
191      K1 = I1MACH(15)
192      K2 = I1MACH(16)
193      R1M5 = D1MACH(5)
194      K = MIN0(IABS(K1),IABS(K2))
195      ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
196      K1 = I1MACH(14) - 1
197      AA = R1M5*DBLE(FLOAT(K1))
198      DIG = DMIN1(AA,18.0D0)
199      AA = AA*2.303D0
200      ALIM = ELIM + DMAX1(-AA,-41.45D0)
201      FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
202      RL = 1.2D0*DIG + 3.0D0
203      FN = FNU + DBLE(FLOAT(NN-1))
204      MM = 3 - M - M
205      FMM = DBLE(FLOAT(MM))
206      ZNR = FMM*ZI
207      ZNI = -FMM*ZR
208C-----------------------------------------------------------------------
209C     TEST FOR PROPER RANGE
210C-----------------------------------------------------------------------
211      AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0)))
212      AA = 0.5D0/TOL
213      BB=DBLE(FLOAT(I1MACH(9)))*0.5D0
214      AA = DMIN1(AA,BB)
215      IF (AZ.GT.AA) GO TO 260
216      IF (FN.GT.AA) GO TO 260
217      AA = DSQRT(AA)
218      IF (AZ.GT.AA) IERR=3
219      IF (FN.GT.AA) IERR=3
220C-----------------------------------------------------------------------
221C     OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
222C-----------------------------------------------------------------------
223      UFL = D1MACH(1)*1.0D+3
224      IF (AZ.LT.UFL) GO TO 230
225      IF (FNU.GT.FNUL) GO TO 90
226      IF (FN.LE.1.0D0) GO TO 70
227      IF (FN.GT.2.0D0) GO TO 60
228      IF (AZ.GT.TOL) GO TO 70
229      ARG = 0.5D0*AZ
230      ALN = -FN*DLOG(ARG)
231      IF (ALN.GT.ELIM) GO TO 230
232      GO TO 70
233   60 CONTINUE
234      CALL ZUOIK(ZNR, ZNI, FNU, KODE, 2, NN, CYR, CYI, NUF, TOL, ELIM,
235     * ALIM)
236      IF (NUF.LT.0) GO TO 230
237      NZ = NZ + NUF
238      NN = NN - NUF
239C-----------------------------------------------------------------------
240C     HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
241C     IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
242C-----------------------------------------------------------------------
243      IF (NN.EQ.0) GO TO 140
244   70 CONTINUE
245      IF ((ZNR.LT.0.0D0) .OR. (ZNR.EQ.0.0D0 .AND. ZNI.LT.0.0D0 .AND.
246     * M.EQ.2)) GO TO 80
247C-----------------------------------------------------------------------
248C     RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR.
249C     YN.GE.0. .OR. M=1)
250C-----------------------------------------------------------------------
251      CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NZ, TOL, ELIM, ALIM)
252      GO TO 110
253C-----------------------------------------------------------------------
254C     LEFT HALF PLANE COMPUTATION
255C-----------------------------------------------------------------------
256   80 CONTINUE
257      MR = -MM
258      CALL ZACON(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, RL, FNUL,
259     * TOL, ELIM, ALIM)
260      IF (NW.LT.0) GO TO 240
261      NZ=NW
262      GO TO 110
263   90 CONTINUE
264C-----------------------------------------------------------------------
265C     UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
266C-----------------------------------------------------------------------
267      MR = 0
268      IF ((ZNR.GE.0.0D0) .AND. (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0 .OR.
269     * M.NE.2)) GO TO 100
270      MR = -MM
271      IF (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0) GO TO 100
272      ZNR = -ZNR
273      ZNI = -ZNI
274  100 CONTINUE
275      CALL ZBUNK(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, TOL, ELIM,
276     * ALIM)
277      IF (NW.LT.0) GO TO 240
278      NZ = NZ + NW
279  110 CONTINUE
280C-----------------------------------------------------------------------
281C     H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT)
282C
283C     ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2
284C-----------------------------------------------------------------------
285      SGN = DSIGN(HPI,-FMM)
286C-----------------------------------------------------------------------
287C     CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
288C     WHEN FNU IS LARGE
289C-----------------------------------------------------------------------
290      INU = INT(SNGL(FNU))
291      INUH = INU/2
292      IR = INU - 2*INUH
293      ARG = (FNU-DBLE(FLOAT(INU-IR)))*SGN
294      RHPI = 1.0D0/SGN
295C     ZNI = RHPI*DCOS(ARG)
296C     ZNR = -RHPI*DSIN(ARG)
297      CSGNI = RHPI*DCOS(ARG)
298      CSGNR = -RHPI*DSIN(ARG)
299      IF (MOD(INUH,2).EQ.0) GO TO 120
300C     ZNR = -ZNR
301C     ZNI = -ZNI
302      CSGNR = -CSGNR
303      CSGNI = -CSGNI
304  120 CONTINUE
305      ZTI = -FMM
306      RTOL = 1.0D0/TOL
307      ASCLE = UFL*RTOL
308      DO 130 I=1,NN
309C       STR = CYR(I)*ZNR - CYI(I)*ZNI
310C       CYI(I) = CYR(I)*ZNI + CYI(I)*ZNR
311C       CYR(I) = STR
312C       STR = -ZNI*ZTI
313C       ZNI = ZNR*ZTI
314C       ZNR = STR
315        AA = CYR(I)
316        BB = CYI(I)
317        ATOL = 1.0D0
318        IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 135
319          AA = AA*RTOL
320          BB = BB*RTOL
321          ATOL = TOL
322  135 CONTINUE
323      STR = AA*CSGNR - BB*CSGNI
324      STI = AA*CSGNI + BB*CSGNR
325      CYR(I) = STR*ATOL
326      CYI(I) = STI*ATOL
327      STR = -CSGNI*ZTI
328      CSGNI = CSGNR*ZTI
329      CSGNR = STR
330  130 CONTINUE
331      RETURN
332  140 CONTINUE
333      IF (ZNR.LT.0.0D0) GO TO 230
334      RETURN
335  230 CONTINUE
336      NZ=0
337      IERR=2
338      RETURN
339  240 CONTINUE
340      IF(NW.EQ.(-1)) GO TO 230
341      NZ=0
342      IERR=5
343      RETURN
344  260 CONTINUE
345      NZ=0
346      IERR=4
347      RETURN
348      END
349