src/numerical/slatec/fortran/zbesj.f

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