src/numerical/slatec/fortran/zbesi.f

   1 *DECK ZBESI
   2       SUBROUTINE ZBESI (ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR)
   3 C***BEGIN PROLOGUE  ZBESI
   4 C***PURPOSE  Compute a sequence of the Bessel functions I(a,z) for
   5 C            complex argument z and real nonnegative orders a=b,b+1,
   6 C            b+2,... where b>0.  A scaling option is available to
   7 C            help avoid overflow.
   8 C***LIBRARY   SLATEC
   9 C***CATEGORY  C10B4
  10 C***TYPE      COMPLEX (CBESI-C, ZBESI-C)
  11 C***KEYWORDS  BESSEL FUNCTIONS OF COMPLEX ARGUMENT, I BESSEL FUNCTIONS,
  12 C             MODIFIED BESSEL FUNCTIONS
  13 C***AUTHOR  Amos, D. E., (SNL)
  14 C***DESCRIPTION
  15 C
  16 C                    ***A DOUBLE PRECISION ROUTINE***
  17 C         On KODE=1, ZBESI computes an N-member sequence of complex
  18 C         Bessel functions CY(L)=I(FNU+L-1,Z) for real nonnegative
  19 C         orders FNU+L-1, L=1,...,N and complex Z in the cut plane
  20 C         -pi<arg(Z)<=pi where Z=ZR+i*ZI.  On KODE=2, CBESI returns
  21 C         the scaled functions
  22 C
  23 C            CY(L) = exp(-abs(X))*I(FNU+L-1,Z), L=1,...,N and X=Re(Z)
  24 C
  25 C         which removes the exponential growth in both the left and
  26 C         right half-planes as Z goes to infinity.
  27 C
  28 C         Input
  29 C           ZR     - DOUBLE PRECISION real part of argument Z
  30 C           ZI     - DOUBLE PRECISION imag part of argument Z
  31 C           FNU    - DOUBLE PRECISION initial order, FNU>=0
  32 C           KODE   - A parameter to indicate the scaling option
  33 C                    KODE=1  returns
  34 C                            CY(L)=I(FNU+L-1,Z), L=1,...,N
  35 C                        =2  returns
  36 C                            CY(L)=exp(-abs(X))*I(FNU+L-1,Z), L=1,...,N
  37 C                            where X=Re(Z)
  38 C           N      - Number of terms in the sequence, N>=1
  39 C
  40 C         Output
  41 C           CYR    - DOUBLE PRECISION real part of result vector
  42 C           CYI    - DOUBLE PRECISION imag part of result vector
  43 C           NZ     - Number of underflows set to zero
  44 C                    NZ=0    Normal return
  45 C                    NZ>0    CY(L)=0, L=N-NZ+1,...,N
  46 C           IERR   - Error flag
  47 C                    IERR=0  Normal return     - COMPUTATION COMPLETED
  48 C                    IERR=1  Input error       - NO COMPUTATION
  49 C                    IERR=2  Overflow          - NO COMPUTATION
  50 C                            (Re(Z) too large on KODE=1)
  51 C                    IERR=3  Precision warning - COMPUTATION COMPLETED
  52 C                            (Result has half precision or less
  53 C                            because abs(Z) or FNU+N-1 is large)
  54 C                    IERR=4  Precision error   - NO COMPUTATION
  55 C                            (Result has no precision because
  56 C                            abs(Z) or FNU+N-1 is too large)
  57 C                    IERR=5  Algorithmic error - NO COMPUTATION
  58 C                            (Termination condition not met)
  59 C
  60 C *Long Description:
  61 C
  62 C         The computation of I(a,z) is carried out by the power series
  63 C         for small abs(z), the asymptotic expansion for large abs(z),
  64 C         the Miller algorithm normalized by the Wronskian and a
  65 C         Neumann series for intermediate magnitudes of z, and the
  66 C         uniform asymptotic expansions for I(a,z) and J(a,z) for
  67 C         large orders a.  Backward recurrence is used to generate
  68 C         sequences or reduce orders when necessary.
  69 C
  70 C         The calculations above are done in the right half plane and
  71 C         continued into the left half plane by the formula
  72 C
  73 C            I(a,z*exp(t)) = exp(t*a)*I(a,z), Re(z)>0
  74 C                        t = i*pi or -i*pi
  75 C
  76 C         For negative orders, the formula
  77 C
  78 C            I(-a,z) = I(a,z) + (2/pi)*sin(pi*a)*K(a,z)
  79 C
  80 C         can be used.  However, for large orders close to integers the
  81 C         the function changes radically.  When a is a large positive
  82 C         integer, the magnitude of I(-a,z)=I(a,z) is a large
  83 C         negative power of ten. But when a is not an integer,
  84 C         K(a,z) dominates in magnitude with a large positive power of
  85 C         ten and the most that the second term can be reduced is by
  86 C         unit roundoff from the coefficient. Thus, wide changes can
  87 C         occur within unit roundoff of a large integer for a. Here,
  88 C         large means a>abs(z).
  89 C
  90 C         In most complex variable computation, one must evaluate ele-
  91 C         mentary functions.  When the magnitude of Z or FNU+N-1 is
  92 C         large, losses of significance by argument reduction occur.
  93 C         Consequently, if either one exceeds U1=SQRT(0.5/UR), then
  94 C         losses exceeding half precision are likely and an error flag
  95 C         IERR=3 is triggered where UR=MAX(D1MACH(4),1.0D-18) is double
  96 C         precision unit roundoff limited to 18 digits precision.  Also,
  97 C         if either is larger than U2=0.5/UR, then all significance is
  98 C         lost and IERR=4.  In order to use the INT function, arguments
  99 C         must be further restricted not to exceed the largest machine
 100 C         integer, U3=I1MACH(9).  Thus, the magnitude of Z and FNU+N-1
 101 C         is restricted by MIN(U2,U3).  In IEEE arithmetic, U1,U2, and
 102 C         U3 approximate 2.0E+3, 4.2E+6, 2.1E+9 in single precision
 103 C         and 4.7E+7, 2.3E+15 and 2.1E+9 in double precision.  This
 104 C         makes U2 limiting in single precision and U3 limiting in
 105 C         double precision.  This means that one can expect to retain,
 106 C         in the worst cases on IEEE machines, no digits in single pre-
 107 C         cision and only 6 digits in double precision.  Similar con-
 108 C         siderations hold for other machines.
 109 C
 110 C         The approximate relative error in the magnitude of a complex
 111 C         Bessel function can be expressed as P*10**S where P=MAX(UNIT
 112 C         ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre-
 113 C         sents the increase in error due to argument reduction in the
 114 C         elementary functions.  Here, S=MAX(1,ABS(LOG10(ABS(Z))),
 115 C         ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF
 116 C         ABS(Z),ABS(EXPONENT OF FNU)) ).  However, the phase angle may
 117 C         have only absolute accuracy.  This is most likely to occur
 118 C         when one component (in magnitude) is larger than the other by
 119 C         several orders of magnitude.  If one component is 10**K larger
 120 C         than the other, then one can expect only MAX(ABS(LOG10(P))-K,
 121 C         0) significant digits; or, stated another way, when K exceeds
 122 C         the exponent of P, no significant digits remain in the smaller
 123 C         component.  However, the phase angle retains absolute accuracy
 124 C         because, in complex arithmetic with precision P, the smaller
 125 C         component will not (as a rule) decrease below P times the
 126 C         magnitude of the larger component.  In these extreme cases,
 127 C         the principal phase angle is on the order of +P, -P, PI/2-P,
 128 C         or -PI/2+P.
 129 C
 130 C***REFERENCES  1. M. Abramowitz and I. A. Stegun, Handbook of Mathe-
 131 C                 matical Functions, National Bureau of Standards
 132 C                 Applied Mathematics Series 55, U. S. Department
 133 C                 of Commerce, Tenth Printing (1972) or later.
 134 C               2. D. E. Amos, Computation of Bessel Functions of
 135 C                 Complex Argument, Report SAND83-0086, Sandia National
 136 C                 Laboratories, Albuquerque, NM, May 1983.
 137 C               3. D. E. Amos, Computation of Bessel Functions of
 138 C                 Complex Argument and Large Order, Report SAND83-0643,
 139 C                 Sandia National Laboratories, Albuquerque, NM, May
 140 C                 1983.
 141 C               4. D. E. Amos, A Subroutine Package for Bessel Functions
 142 C                 of a Complex Argument and Nonnegative Order, Report
 143 C                 SAND85-1018, Sandia National Laboratory, Albuquerque,
 144 C                 NM, May 1985.
 145 C               5. D. E. Amos, A portable package for Bessel functions
 146 C                 of a complex argument and nonnegative order, ACM
 147 C                 Transactions on Mathematical Software, 12 (September
 148 C                 1986), pp. 265-273.
 149 C
 150 C***ROUTINES CALLED  D1MACH, I1MACH, ZABS, ZBINU
 151 C***REVISION HISTORY  (YYMMDD)
 152 C   830501  DATE WRITTEN
 153 C   890801  REVISION DATE from Version 3.2
 154 C   910415  Prologue converted to Version 4.0 format.  (BAB)
 155 C   920128  Category corrected.  (WRB)
 156 C   920811  Prologue revised.  (DWL)
 157 C***END PROLOGUE  ZBESI
 158 C     COMPLEX CONE,CSGN,CW,CY,CZERO,Z,ZN
 159       DOUBLE PRECISION AA, ALIM, ARG, CONEI, CONER, CSGNI, CSGNR, CYI,
 160      * CYR, DIG, ELIM, FNU, FNUL, PI, RL, R1M5, STR, TOL, ZI, ZNI, ZNR,
 161      * ZR, D1MACH, AZ, BB, FN, ZABS, ASCLE, RTOL, ATOL, STI
 162       INTEGER I, IERR, INU, K, KODE, K1,K2,N,NZ,NN, I1MACH
 163       DIMENSION CYR(N), CYI(N)
 164       EXTERNAL ZABS
 165       DATA PI /3.14159265358979324D0/
 166       DATA CONER, CONEI /1.0D0,0.0D0/
 167 C
 168 C***FIRST EXECUTABLE STATEMENT  ZBESI
 169       IERR = 0
 170       NZ=0
 171       IF (FNU.LT.0.0D0) IERR=1
 172       IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
 173       IF (N.LT.1) IERR=1
 174       IF (IERR.NE.0) RETURN
 175 C-----------------------------------------------------------------------
 176 C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
 177 C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
 178 C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
 179 C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
 180 C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
 181 C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
 182 C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
 183 C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
 184 C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU.
 185 C-----------------------------------------------------------------------
 186       TOL = MAX(D1MACH(4),1.0D-18)
 187       K1 = I1MACH(15)
 188       K2 = I1MACH(16)
 189       R1M5 = D1MACH(5)
 190       K = MIN(ABS(K1),ABS(K2))
 191       ELIM = 2.303D0*(K*R1M5-3.0D0)
 192       K1 = I1MACH(14) - 1
 193       AA = R1M5*K1
 194       DIG = MIN(AA,18.0D0)
 195       AA = AA*2.303D0
 196       ALIM = ELIM + MAX(-AA,-41.45D0)
 197       RL = 1.2D0*DIG + 3.0D0
 198       FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
 199 C-----------------------------------------------------------------------
 200 C     TEST FOR PROPER RANGE
 201 C-----------------------------------------------------------------------
 202       AZ = ZABS(ZR,ZI)
 203       FN = FNU+(N-1)
 204       AA = 0.5D0/TOL
 205       BB=I1MACH(9)*0.5D0
 206       AA = MIN(AA,BB)
 207       IF (AZ.GT.AA) GO TO 260
 208       IF (FN.GT.AA) GO TO 260
 209       AA = SQRT(AA)
 210       IF (AZ.GT.AA) IERR=3
 211       IF (FN.GT.AA) IERR=3
 212       ZNR = ZR
 213       ZNI = ZI
 214       CSGNR = CONER
 215       CSGNI = CONEI
 216       IF (ZR.GE.0.0D0) GO TO 40
 217       ZNR = -ZR
 218       ZNI = -ZI
 219 C-----------------------------------------------------------------------
 220 C     CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
 221 C     WHEN FNU IS LARGE
 222 C-----------------------------------------------------------------------
 223       INU = FNU
 224       ARG = (FNU-INU)*PI
 225       IF (ZI.LT.0.0D0) ARG = -ARG
 226       CSGNR = COS(ARG)
 227       CSGNI = SIN(ARG)
 228       IF (MOD(INU,2).EQ.0) GO TO 40
 229       CSGNR = -CSGNR
 230       CSGNI = -CSGNI
 231    40 CONTINUE
 232       CALL ZBINU(ZNR, ZNI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL,
 233      * ELIM, ALIM)
 234       IF (NZ.LT.0) GO TO 120
 235       IF (ZR.GE.0.0D0) RETURN
 236 C-----------------------------------------------------------------------
 237 C     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE
 238 C-----------------------------------------------------------------------
 239       NN = N - NZ
 240       IF (NN.EQ.0) RETURN
 241       RTOL = 1.0D0/TOL
 242       ASCLE = D1MACH(1)*RTOL*1.0D+3
 243       DO 50 I=1,NN
 244 C       STR = CYR(I)*CSGNR - CYI(I)*CSGNI
 245 C       CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR
 246 C       CYR(I) = STR
 247         AA = CYR(I)
 248         BB = CYI(I)
 249         ATOL = 1.0D0
 250         IF (MAX(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 55
 251           AA = AA*RTOL
 252           BB = BB*RTOL
 253           ATOL = TOL
 254    55   CONTINUE
 255         STR = AA*CSGNR - BB*CSGNI
 256         STI = AA*CSGNI + BB*CSGNR
 257         CYR(I) = STR*ATOL
 258         CYI(I) = STI*ATOL
 259         CSGNR = -CSGNR
 260         CSGNI = -CSGNI
 261    50 CONTINUE
 262       RETURN
 263   120 CONTINUE
 264       IF(NZ.EQ.(-2)) GO TO 130
 265       NZ = 0
 266       IERR=2
 267       RETURN
 268   130 CONTINUE
 269       NZ=0
 270       IERR=5
 271       RETURN
 272   260 CONTINUE
 273       NZ=0
 274       IERR=4
 275       RETURN
 276       END