amos/zbesh.f

   1       SUBROUTINE ZBESH(ZR, ZI, FNU, KODE, M, N, CYR, CYI, NZ, IERR)
   2 C***BEGIN PROLOGUE  ZBESH
   3 C***DATE WRITTEN   830501   (YYMMDD)
   4 C***REVISION DATE  890801   (YYMMDD)
   5 C***CATEGORY NO.  B5K
   6 C***KEYWORDS  H-BESSEL FUNCTIONS,BESSEL FUNCTIONS OF COMPLEX ARGUMENT,
   7 C             BESSEL FUNCTIONS OF THIRD KIND,HANKEL FUNCTIONS
   8 C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
   9 C***PURPOSE  TO COMPUTE THE H-BESSEL FUNCTIONS OF A COMPLEX ARGUMENT
  10 C***DESCRIPTION
  11 C
  12 C                      ***A DOUBLE PRECISION ROUTINE***
  13 C         ON KODE=1, ZBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
  14 C         HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1
  15 C         OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX
  16 C         Z.NE.CMPLX(0.0,0.0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI.
  17 C         ON KODE=2, ZBESH RETURNS THE SCALED HANKEL FUNCTIONS
  18 C
  19 C         CY(I)=EXP(-MM*Z*I)*H(M,FNU+J-1,Z)       MM=3-2*M,   I**2=-1.
  20 C
  21 C         WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER AND
  22 C         LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN THE
  23 C         NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1).
  24 C
  25 C         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
  26 C           ZR,ZI  - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0),
  27 C                    -PT.LT.ARG(Z).LE.PI
  28 C           FNU    - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0D0
  29 C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
  30 C                    KODE= 1  RETURNS
  31 C                             CY(J)=H(M,FNU+J-1,Z),   J=1,...,N
  32 C                        = 2  RETURNS
  33 C                             CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M))
  34 C                                  J=1,...,N  ,  I**2=-1
  35 C           M      - KIND OF HANKEL FUNCTION, M=1 OR 2
  36 C           N      - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1
  37 C
  38 C         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
  39 C           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
  40 C                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
  41 C                    CY(J)=H(M,FNU+J-1,Z)  OR
  42 C                    CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M))  J=1,...,N
  43 C                    DEPENDING ON KODE, I**2=-1.
  44 C           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
  45 C                    NZ= 0   , NORMAL RETURN
  46 C                    NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE
  47 C                              TO UNDERFLOW, CY(J)=CMPLX(0.0D0,0.0D0)
  48 C                              J=1,...,NZ WHEN Y.GT.0.0 AND M=1 OR
  49 C                              Y.LT.0.0 AND M=2. FOR THE COMPLMENTARY
  50 C                              HALF PLANES, NZ STATES ONLY THE NUMBER
  51 C                              OF UNDERFLOWS.
  52 C           IERR   - ERROR FLAG
  53 C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
  54 C                    IERR=1, INPUT ERROR   - NO COMPUTATION
  55 C                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU TOO
  56 C                            LARGE OR CABS(Z) TOO SMALL OR BOTH
  57 C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
  58 C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
  59 C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
  60 C                            ACCURACY
  61 C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
  62 C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
  63 C                            CANCE BY ARGUMENT REDUCTION
  64 C                    IERR=5, ERROR              - NO COMPUTATION,
  65 C                            ALGORITHM TERMINATION CONDITION NOT MET
  66 C
  67 C***LONG DESCRIPTION
  68 C
  69 C         THE COMPUTATION IS CARRIED OUT BY THE RELATION
  70 C
  71 C         H(M,FNU,Z)=(1/MP)*EXP(-MP*FNU)*K(FNU,Z*EXP(-MP))
  72 C             MP=MM*HPI*I,  MM=3-2*M,  HPI=PI/2,  I**2=-1
  73 C
  74 C         FOR M=1 OR 2 WHERE THE K BESSEL FUNCTION IS COMPUTED FOR THE
  75 C         RIGHT HALF PLANE RE(Z).GE.0.0. THE K FUNCTION IS CONTINUED
  76 C         TO THE LEFT HALF PLANE BY THE RELATION
  77 C
  78 C         K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z)
  79 C         MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1
  80 C
  81 C         WHERE I(FNU,Z) IS THE I BESSEL FUNCTION.
  82 C
  83 C         EXPONENTIAL DECAY OF H(M,FNU,Z) OCCURS IN THE UPPER HALF Z
  84 C         PLANE FOR M=1 AND THE LOWER HALF Z PLANE FOR M=2.  EXPONENTIAL
  85 C         GROWTH OCCURS IN THE COMPLEMENTARY HALF PLANES.  SCALING
  86 C         BY EXP(-MM*Z*I) REMOVES THE EXPONENTIAL BEHAVIOR IN THE
  87 C         WHOLE Z PLANE FOR Z TO INFINITY.
  88 C
  89 C         FOR NEGATIVE ORDERS,THE FORMULAE
  90 C
  91 C               H(1,-FNU,Z) = H(1,FNU,Z)*CEXP( PI*FNU*I)
  92 C               H(2,-FNU,Z) = H(2,FNU,Z)*CEXP(-PI*FNU*I)
  93 C                         I**2=-1
  94 C
  95 C         CAN BE USED.
  96 C
  97 C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
  98 C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
  99 C         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
 100 C         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
 101 C         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
 102 C         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
 103 C         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
 104 C         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
 105 C         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
 106 C         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
 107 C         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
 108 C         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
 109 C         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
 110 C         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
 111 C         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
 112 C         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
 113 C         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
 114 C         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
 115 C         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
 116 C
 117 C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
 118 C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
 119 C         ROUNDOFF,1.0D-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
 120 C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
 121 C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
 122 C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
 123 C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
 124 C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
 125 C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
 126 C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
 127 C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
 128 C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
 129 C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
 130 C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
 131 C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
 132 C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
 133 C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
 134 C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
 135 C         OR -PI/2+P.
 136 C
 137 C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
 138 C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
 139 C                 COMMERCE, 1955.
 140 C
 141 C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
 142 C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
 143 C
 144 C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
 145 C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
 146 C
 147 C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
 148 C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
 149 C                 1018, MAY, 1985
 150 C
 151 C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
 152 C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
 153 C                 MATH. SOFTWARE, 1986
 154 C
 155 C***ROUTINES CALLED  ZACON,ZBKNU,ZBUNK,ZUOIK,AZABS,I1MACH,D1MACH
 156 C***END PROLOGUE  ZBESH
 157 C
 158 C     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, AZABS, 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)
 166 C
 167       DATA HPI /1.57079632679489662D0/
 168 C
 169 C***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
 179 C-----------------------------------------------------------------------
 180 C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
 181 C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
 182 C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
 183 C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
 184 C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
 185 C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
 186 C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
 187 C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
 188 C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
 189 C-----------------------------------------------------------------------
 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
 208 C-----------------------------------------------------------------------
 209 C     TEST FOR PROPER RANGE
 210 C-----------------------------------------------------------------------
 211       AZ = AZABS(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
 220 C-----------------------------------------------------------------------
 221 C     OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
 222 C-----------------------------------------------------------------------
 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
 239 C-----------------------------------------------------------------------
 240 C     HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
 241 C     IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
 242 C-----------------------------------------------------------------------
 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
 247 C-----------------------------------------------------------------------
 248 C     RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR.
 249 C     YN.GE.0. .OR. M=1)
 250 C-----------------------------------------------------------------------
 251       CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NZ, TOL, ELIM, ALIM)
 252       GO TO 110
 253 C-----------------------------------------------------------------------
 254 C     LEFT HALF PLANE COMPUTATION
 255 C-----------------------------------------------------------------------
 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
 264 C-----------------------------------------------------------------------
 265 C     UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
 266 C-----------------------------------------------------------------------
 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
 280 C-----------------------------------------------------------------------
 281 C     H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT)
 282 C
 283 C     ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2
 284 C-----------------------------------------------------------------------
 285       SGN = DSIGN(HPI,-FMM)
 286 C-----------------------------------------------------------------------
 287 C     CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
 288 C     WHEN FNU IS LARGE
 289 C-----------------------------------------------------------------------
 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
 295 C     ZNI = RHPI*DCOS(ARG)
 296 C     ZNR = -RHPI*DSIN(ARG)
 297       CSGNI = RHPI*DCOS(ARG)
 298       CSGNR = -RHPI*DSIN(ARG)
 299       IF (MOD(INUH,2).EQ.0) GO TO 120
 300 C     ZNR = -ZNR
 301 C     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
 309 C       STR = CYR(I)*ZNR - CYI(I)*ZNI
 310 C       CYI(I) = CYR(I)*ZNI + CYI(I)*ZNR
 311 C       CYR(I) = STR
 312 C       STR = -ZNI*ZTI
 313 C       ZNI = ZNR*ZTI
 314 C       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