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,XZABS,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, XZABS, 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 = XZABS(ZR,ZI)
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   35 CONTINUE
224      UFL = D1MACH(1)*1.0D+3
225      IF (AZ.LT.UFL) GO TO 230
226      IF (FNU.GT.FNUL) GO TO 90
227      IF (FN.LE.1.0D0) GO TO 70
228      IF (FN.GT.2.0D0) GO TO 60
229      IF (AZ.GT.TOL) GO TO 70
230      ARG = 0.5D0*AZ
231      ALN = -FN*DLOG(ARG)
232      IF (ALN.GT.ELIM) GO TO 230
233      GO TO 70
234   60 CONTINUE
235      CALL ZUOIK(ZNR, ZNI, FNU, KODE, 2, NN, CYR, CYI, NUF, TOL, ELIM,
236     * ALIM)
237      IF (NUF.LT.0) GO TO 230
238      NZ = NZ + NUF
239      NN = NN - NUF
240C-----------------------------------------------------------------------
241C     HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
242C     IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
243C-----------------------------------------------------------------------
244      IF (NN.EQ.0) GO TO 140
245   70 CONTINUE
246      IF ((ZNR.LT.0.0D0) .OR. (ZNR.EQ.0.0D0 .AND. ZNI.LT.0.0D0 .AND.
247     * M.EQ.2)) GO TO 80
248C-----------------------------------------------------------------------
249C     RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR.
250C     YN.GE.0. .OR. M=1)
251C-----------------------------------------------------------------------
252      CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NZ, TOL, ELIM, ALIM)
253      GO TO 110
254C-----------------------------------------------------------------------
255C     LEFT HALF PLANE COMPUTATION
256C-----------------------------------------------------------------------
257   80 CONTINUE
258      MR = -MM
259      CALL ZACON(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, RL, FNUL,
260     * TOL, ELIM, ALIM)
261      IF (NW.LT.0) GO TO 240
262      NZ=NW
263      GO TO 110
264   90 CONTINUE
265C-----------------------------------------------------------------------
266C     UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
267C-----------------------------------------------------------------------
268      MR = 0
269      IF ((ZNR.GE.0.0D0) .AND. (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0 .OR.
270     * M.NE.2)) GO TO 100
271      MR = -MM
272      IF (ZNR.NE.0.0D0 .OR. ZNI.GE.0.0D0) GO TO 100
273      ZNR = -ZNR
274      ZNI = -ZNI
275  100 CONTINUE
276      CALL ZBUNK(ZNR, ZNI, FNU, KODE, MR, NN, CYR, CYI, NW, TOL, ELIM,
277     * ALIM)
278      IF (NW.LT.0) GO TO 240
279      NZ = NZ + NW
280  110 CONTINUE
281C-----------------------------------------------------------------------
282C     H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT)
283C
284C     ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2
285C-----------------------------------------------------------------------
286      SGN = DSIGN(HPI,-FMM)
287C-----------------------------------------------------------------------
288C     CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
289C     WHEN FNU IS LARGE
290C-----------------------------------------------------------------------
291      INU = INT(SNGL(FNU))
292      INUH = INU/2
293      IR = INU - 2*INUH
294      ARG = (FNU-DBLE(FLOAT(INU-IR)))*SGN
295      RHPI = 1.0D0/SGN
296C     ZNI = RHPI*DCOS(ARG)
297C     ZNR = -RHPI*DSIN(ARG)
298      CSGNI = RHPI*DCOS(ARG)
299      CSGNR = -RHPI*DSIN(ARG)
300      IF (MOD(INUH,2).EQ.0) GO TO 120
301C     ZNR = -ZNR
302C     ZNI = -ZNI
303      CSGNR = -CSGNR
304      CSGNI = -CSGNI
305  120 CONTINUE
306      ZTI = -FMM
307      RTOL = 1.0D0/TOL
308      ASCLE = UFL*RTOL
309      DO 130 I=1,NN
310C       STR = CYR(I)*ZNR - CYI(I)*ZNI
311C       CYI(I) = CYR(I)*ZNI + CYI(I)*ZNR
312C       CYR(I) = STR
313C       STR = -ZNI*ZTI
314C       ZNI = ZNR*ZTI
315C       ZNR = STR
316        AA = CYR(I)
317        BB = CYI(I)
318        ATOL = 1.0D0
319        IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 135
320          AA = AA*RTOL
321          BB = BB*RTOL
322          ATOL = TOL
323  135 CONTINUE
324      STR = AA*CSGNR - BB*CSGNI
325      STI = AA*CSGNI + BB*CSGNR
326      CYR(I) = STR*ATOL
327      CYI(I) = STI*ATOL
328      STR = -CSGNI*ZTI
329      CSGNI = CSGNR*ZTI
330      CSGNR = STR
331  130 CONTINUE
332      RETURN
333  140 CONTINUE
334      IF (ZNR.LT.0.0D0) GO TO 230
335      RETURN
336  230 CONTINUE
337      NZ=0
338      IERR=2
339      RETURN
340  240 CONTINUE
341      IF(NW.EQ.(-1)) GO TO 230
342      NZ=0
343      IERR=5
344      RETURN
345  260 CONTINUE
346      IERR=4
347      GO TO 35
348      END
349