src/numerical/slatec/fortran/dbesi.f

   1 *DECK DBESI
   2       SUBROUTINE DBESI (X, ALPHA, KODE, N, Y, NZ)
   3 C***BEGIN PROLOGUE  DBESI
   4 C***PURPOSE  Compute an N member sequence of I Bessel functions
   5 C            I/SUB(ALPHA+K-1)/(X), K=1,...,N or scaled Bessel functions
   6 C            EXP(-X)*I/SUB(ALPHA+K-1)/(X), K=1,...,N for nonnegative
   7 C            ALPHA and X.
   8 C***LIBRARY   SLATEC
   9 C***CATEGORY  C10B3
  10 C***TYPE      DOUBLE PRECISION (BESI-S, DBESI-D)
  11 C***KEYWORDS  I BESSEL FUNCTION, SPECIAL FUNCTIONS
  12 C***AUTHOR  Amos, D. E., (SNLA)
  13 C           Daniel, S. L., (SNLA)
  14 C***DESCRIPTION
  15 C
  16 C     Abstract  **** a double precision routine ****
  17 C         DBESI computes an N member sequence of I Bessel functions
  18 C         I/sub(ALPHA+K-1)/(X), K=1,...,N or scaled Bessel functions
  19 C         EXP(-X)*I/sub(ALPHA+K-1)/(X), K=1,...,N for nonnegative ALPHA
  20 C         and X.  A combination of the power series, the asymptotic
  21 C         expansion for X to infinity, and the uniform asymptotic
  22 C         expansion for NU to infinity are applied over subdivisions of
  23 C         the (NU,X) plane.  For values not covered by one of these
  24 C         formulae, the order is incremented by an integer so that one
  25 C         of these formulae apply.  Backward recursion is used to reduce
  26 C         orders by integer values.  The asymptotic expansion for X to
  27 C         infinity is used only when the entire sequence (specifically
  28 C         the last member) lies within the region covered by the
  29 C         expansion.  Leading terms of these expansions are used to test
  30 C         for over or underflow where appropriate.  If a sequence is
  31 C         requested and the last member would underflow, the result is
  32 C         set to zero and the next lower order tried, etc., until a
  33 C         member comes on scale or all are set to zero.  An overflow
  34 C         cannot occur with scaling.
  35 C
  36 C         The maximum number of significant digits obtainable
  37 C         is the smaller of 14 and the number of digits carried in
  38 C         double precision arithmetic.
  39 C
  40 C     Description of Arguments
  41 C
  42 C         Input      X,ALPHA are double precision
  43 C           X      - X .GE. 0.0D0
  44 C           ALPHA  - order of first member of the sequence,
  45 C                    ALPHA .GE. 0.0D0
  46 C           KODE   - a parameter to indicate the scaling option
  47 C                    KODE=1 returns
  48 C                           Y(K)=        I/sub(ALPHA+K-1)/(X),
  49 C                                K=1,...,N
  50 C                    KODE=2 returns
  51 C                           Y(K)=EXP(-X)*I/sub(ALPHA+K-1)/(X),
  52 C                                K=1,...,N
  53 C           N      - number of members in the sequence, N .GE. 1
  54 C
  55 C         Output     Y is double precision
  56 C           Y      - a vector whose first N components contain
  57 C                    values for I/sub(ALPHA+K-1)/(X) or scaled
  58 C                    values for EXP(-X)*I/sub(ALPHA+K-1)/(X),
  59 C                    K=1,...,N depending on KODE
  60 C           NZ     - number of components of Y set to zero due to
  61 C                    underflow,
  62 C                    NZ=0   , normal return, computation completed
  63 C                    NZ .NE. 0, last NZ components of Y set to zero,
  64 C                             Y(K)=0.0D0, K=N-NZ+1,...,N.
  65 C
  66 C     Error Conditions
  67 C         Improper input arguments - a fatal error
  68 C         Overflow with KODE=1 - a fatal error
  69 C         Underflow - a non-fatal error(NZ .NE. 0)
  70 C
  71 C***REFERENCES  D. E. Amos, S. L. Daniel and M. K. Weston, CDC 6600
  72 C                 subroutines IBESS and JBESS for Bessel functions
  73 C                 I(NU,X) and J(NU,X), X .GE. 0, NU .GE. 0, ACM
  74 C                 Transactions on Mathematical Software 3, (1977),
  75 C                 pp. 76-92.
  76 C               F. W. J. Olver, Tables of Bessel Functions of Moderate
  77 C                 or Large Orders, NPL Mathematical Tables 6, Her
  78 C                 Majesty's Stationery Office, London, 1962.
  79 C***ROUTINES CALLED  D1MACH, DASYIK, DLNGAM, I1MACH, XERMSG
  80 C***REVISION HISTORY  (YYMMDD)
  81 C   750101  DATE WRITTEN
  82 C   890531  Changed all specific intrinsics to generic.  (WRB)
  83 C   890911  Removed unnecessary intrinsics.  (WRB)
  84 C   890911  REVISION DATE from Version 3.2
  85 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  86 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
  87 C   900326  Removed duplicate information from DESCRIPTION section.
  88 C           (WRB)
  89 C   920501  Reformatted the REFERENCES section.  (WRB)
  90 C***END PROLOGUE  DBESI
  91 C
  92       INTEGER I, IALP, IN, INLIM, IS, I1, K, KK, KM, KODE, KT,
  93      1 N, NN, NS, NZ
  94       INTEGER I1MACH
  95       DOUBLE PRECISION AIN,AK,AKM,ALPHA,ANS,AP,ARG,ATOL,TOLLN,DFN,
  96      1 DTM, DX, EARG, ELIM, ETX, FLGIK,FN, FNF, FNI,FNP1,FNU,GLN,RA,
  97      2 RTTPI, S, SX, SXO2, S1, S2, T, TA, TB, TEMP, TFN, TM, TOL,
  98      3 TRX, T2, X, XO2, XO2L, Y, Z
  99       DOUBLE PRECISION D1MACH, DLNGAM
 100       DIMENSION Y(*), TEMP(3)
 101       SAVE RTTPI, INLIM
 102       DATA RTTPI           / 3.98942280401433D-01/
 103       DATA INLIM           /          80         /
 104 C***FIRST EXECUTABLE STATEMENT  DBESI
 105       NZ = 0
 106       KT = 1
 107 C     I1MACH(15) REPLACES I1MACH(12) IN A DOUBLE PRECISION CODE
 108 C     I1MACH(14) REPLACES I1MACH(11) IN A DOUBLE PRECISION CODE
 109       RA = D1MACH(3)
 110       TOL = MAX(RA,1.0D-15)
 111       I1 = -I1MACH(15)
 112       GLN = D1MACH(5)
 113       ELIM = 2.303D0*(I1*GLN-3.0D0)
 114 C     TOLLN = -LN(TOL)
 115       I1 = I1MACH(14)+1
 116       TOLLN = 2.303D0*GLN*I1
 117       TOLLN = MIN(TOLLN,34.5388D0)
 118       IF (N-1) 590, 10, 20
 119    10 KT = 2
 120    20 NN = N
 121       IF (KODE.LT.1 .OR. KODE.GT.2) GO TO 570
 122       IF (X) 600, 30, 80
 123    30 IF (ALPHA) 580, 40, 50
 124    40 Y(1) = 1.0D0
 125       IF (N.EQ.1) RETURN
 126       I1 = 2
 127       GO TO 60
 128    50 I1 = 1
 129    60 DO 70 I=I1,N
 130         Y(I) = 0.0D0
 131    70 CONTINUE
 132       RETURN
 133    80 CONTINUE
 134       IF (ALPHA.LT.0.0D0) GO TO 580
 135 C
 136       IALP = INT(ALPHA)
 137       FNI = IALP + N - 1
 138       FNF = ALPHA - IALP
 139       DFN = FNI + FNF
 140       FNU = DFN
 141       IN = 0
 142       XO2 = X*0.5D0
 143       SXO2 = XO2*XO2
 144       ETX = KODE - 1
 145       SX = ETX*X
 146 C
 147 C     DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X
 148 C     TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE
 149 C     APPLIED.
 150 C
 151       IF (SXO2.LE.(FNU+1.0D0)) GO TO 90
 152       IF (X.LE.12.0D0) GO TO 110
 153       FN = 0.55D0*FNU*FNU
 154       FN = MAX(17.0D0,FN)
 155       IF (X.GE.FN) GO TO 430
 156       ANS = MAX(36.0D0-FNU,0.0D0)
 157       NS = INT(ANS)
 158       FNI = FNI + NS
 159       DFN = FNI + FNF
 160       FN = DFN
 161       IS = KT
 162       KM = N - 1 + NS
 163       IF (KM.GT.0) IS = 3
 164       GO TO 120
 165    90 FN = FNU
 166       FNP1 = FN + 1.0D0
 167       XO2L = LOG(XO2)
 168       IS = KT
 169       IF (X.LE.0.5D0) GO TO 230
 170       NS = 0
 171   100 FNI = FNI + NS
 172       DFN = FNI + FNF
 173       FN = DFN
 174       FNP1 = FN + 1.0D0
 175       IS = KT
 176       IF (N-1+NS.GT.0) IS = 3
 177       GO TO 230
 178   110 XO2L = LOG(XO2)
 179       NS = INT(SXO2-FNU)
 180       GO TO 100
 181   120 CONTINUE
 182 C
 183 C     OVERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION
 184 C
 185       IF (KODE.EQ.2) GO TO 130
 186       IF (ALPHA.LT.1.0D0) GO TO 150
 187       Z = X/ALPHA
 188       RA = SQRT(1.0D0+Z*Z)
 189       GLN = LOG((1.0D0+RA)/Z)
 190       T = RA*(1.0D0-ETX) + ETX/(Z+RA)
 191       ARG = ALPHA*(T-GLN)
 192       IF (ARG.GT.ELIM) GO TO 610
 193       IF (KM.EQ.0) GO TO 140
 194   130 CONTINUE
 195 C
 196 C     UNDERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION
 197 C
 198       Z = X/FN
 199       RA = SQRT(1.0D0+Z*Z)
 200       GLN = LOG((1.0D0+RA)/Z)
 201       T = RA*(1.0D0-ETX) + ETX/(Z+RA)
 202       ARG = FN*(T-GLN)
 203   140 IF (ARG.LT.(-ELIM)) GO TO 280
 204       GO TO 190
 205   150 IF (X.GT.ELIM) GO TO 610
 206       GO TO 130
 207 C
 208 C     UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY
 209 C
 210   160 IF (KM.NE.0) GO TO 170
 211       Y(1) = TEMP(3)
 212       RETURN
 213   170 TEMP(1) = TEMP(3)
 214       IN = NS
 215       KT = 1
 216       I1 = 0
 217   180 CONTINUE
 218       IS = 2
 219       FNI = FNI - 1.0D0
 220       DFN = FNI + FNF
 221       FN = DFN
 222       IF(I1.EQ.2) GO TO 350
 223       Z = X/FN
 224       RA = SQRT(1.0D0+Z*Z)
 225       GLN = LOG((1.0D0+RA)/Z)
 226       T = RA*(1.0D0-ETX) + ETX/(Z+RA)
 227       ARG = FN*(T-GLN)
 228   190 CONTINUE
 229       I1 = ABS(3-IS)
 230       I1 = MAX(I1,1)
 231       FLGIK = 1.0D0
 232       CALL DASYIK(X,FN,KODE,FLGIK,RA,ARG,I1,TEMP(IS))
 233       GO TO (180, 350, 510), IS
 234 C
 235 C     SERIES FOR (X/2)**2.LE.NU+1
 236 C
 237   230 CONTINUE
 238       GLN = DLNGAM(FNP1)
 239       ARG = FN*XO2L - GLN - SX
 240       IF (ARG.LT.(-ELIM)) GO TO 300
 241       EARG = EXP(ARG)
 242   240 CONTINUE
 243       S = 1.0D0
 244       IF (X.LT.TOL) GO TO 260
 245       AK = 3.0D0
 246       T2 = 1.0D0
 247       T = 1.0D0
 248       S1 = FN
 249       DO 250 K=1,17
 250         S2 = T2 + S1
 251         T = T*SXO2/S2
 252         S = S + T
 253         IF (ABS(T).LT.TOL) GO TO 260
 254         T2 = T2 + AK
 255         AK = AK + 2.0D0
 256         S1 = S1 + FN
 257   250 CONTINUE
 258   260 CONTINUE
 259       TEMP(IS) = S*EARG
 260       GO TO (270, 350, 500), IS
 261   270 EARG = EARG*FN/XO2
 262       FNI = FNI - 1.0D0
 263       DFN = FNI + FNF
 264       FN = DFN
 265       IS = 2
 266       GO TO 240
 267 C
 268 C     SET UNDERFLOW VALUE AND UPDATE PARAMETERS
 269 C
 270   280 Y(NN) = 0.0D0
 271       NN = NN - 1
 272       FNI = FNI - 1.0D0
 273       DFN = FNI + FNF
 274       FN = DFN
 275       IF (NN-1) 340, 290, 130
 276   290 KT = 2
 277       IS = 2
 278       GO TO 130
 279   300 Y(NN) = 0.0D0
 280       NN = NN - 1
 281       FNP1 = FN
 282       FNI = FNI - 1.0D0
 283       DFN = FNI + FNF
 284       FN = DFN
 285       IF (NN-1) 340, 310, 320
 286   310 KT = 2
 287       IS = 2
 288   320 IF (SXO2.LE.FNP1) GO TO 330
 289       GO TO 130
 290   330 ARG = ARG - XO2L + LOG(FNP1)
 291       IF (ARG.LT.(-ELIM)) GO TO 300
 292       GO TO 230
 293   340 NZ = N - NN
 294       RETURN
 295 C
 296 C     BACKWARD RECURSION SECTION
 297 C
 298   350 CONTINUE
 299       NZ = N - NN
 300   360 CONTINUE
 301       IF(KT.EQ.2) GO TO 420
 302       S1 = TEMP(1)
 303       S2 = TEMP(2)
 304       TRX = 2.0D0/X
 305       DTM = FNI
 306       TM = (DTM+FNF)*TRX
 307       IF (IN.EQ.0) GO TO 390
 308 C     BACKWARD RECUR TO INDEX ALPHA+NN-1
 309       DO 380 I=1,IN
 310         S = S2
 311         S2 = TM*S2 + S1
 312         S1 = S
 313         DTM = DTM - 1.0D0
 314         TM = (DTM+FNF)*TRX
 315   380 CONTINUE
 316       Y(NN) = S1
 317       IF (NN.EQ.1) RETURN
 318       Y(NN-1) = S2
 319       IF (NN.EQ.2) RETURN
 320       GO TO 400
 321   390 CONTINUE
 322 C     BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
 323       Y(NN) = S1
 324       Y(NN-1) = S2
 325       IF (NN.EQ.2) RETURN
 326   400 K = NN + 1
 327       DO 410 I=3,NN
 328         K = K - 1
 329         Y(K-2) = TM*Y(K-1) + Y(K)
 330         DTM = DTM - 1.0D0
 331         TM = (DTM+FNF)*TRX
 332   410 CONTINUE
 333       RETURN
 334   420 Y(1) = TEMP(2)
 335       RETURN
 336 C
 337 C     ASYMPTOTIC EXPANSION FOR X TO INFINITY
 338 C
 339   430 CONTINUE
 340       EARG = RTTPI/SQRT(X)
 341       IF (KODE.EQ.2) GO TO 440
 342       IF (X.GT.ELIM) GO TO 610
 343       EARG = EARG*EXP(X)
 344   440 ETX = 8.0D0*X
 345       IS = KT
 346       IN = 0
 347       FN = FNU
 348   450 DX = FNI + FNI
 349       TM = 0.0D0
 350       IF (FNI.EQ.0.0D0 .AND. ABS(FNF).LT.TOL) GO TO 460
 351       TM = 4.0D0*FNF*(FNI+FNI+FNF)
 352   460 CONTINUE
 353       DTM = DX*DX
 354       S1 = ETX
 355       TRX = DTM - 1.0D0
 356       DX = -(TRX+TM)/ETX
 357       T = DX
 358       S = 1.0D0 + DX
 359       ATOL = TOL*ABS(S)
 360       S2 = 1.0D0
 361       AK = 8.0D0
 362       DO 470 K=1,25
 363         S1 = S1 + ETX
 364         S2 = S2 + AK
 365         DX = DTM - S2
 366         AP = DX + TM
 367         T = -T*AP/S1
 368         S = S + T
 369         IF (ABS(T).LE.ATOL) GO TO 480
 370         AK = AK + 8.0D0
 371   470 CONTINUE
 372   480 TEMP(IS) = S*EARG
 373       IF(IS.EQ.2) GO TO 360
 374       IS = 2
 375       FNI = FNI - 1.0D0
 376       DFN = FNI + FNF
 377       FN = DFN
 378       GO TO 450
 379 C
 380 C     BACKWARD RECURSION WITH NORMALIZATION BY
 381 C     ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES.
 382 C
 383   500 CONTINUE
 384 C     COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION
 385       AKM = MAX(3.0D0-FN,0.0D0)
 386       KM = INT(AKM)
 387       TFN = FN + KM
 388       TA = (GLN+TFN-0.9189385332D0-0.0833333333D0/TFN)/(TFN+0.5D0)
 389       TA = XO2L - TA
 390       TB = -(1.0D0-1.0D0/TFN)/TFN
 391       AIN = TOLLN/(-TA+SQRT(TA*TA-TOLLN*TB)) + 1.5D0
 392       IN = INT(AIN)
 393       IN = IN + KM
 394       GO TO 520
 395   510 CONTINUE
 396 C     COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION
 397       T = 1.0D0/(FN*RA)
 398       AIN = TOLLN/(GLN+SQRT(GLN*GLN+T*TOLLN)) + 1.5D0
 399       IN = INT(AIN)
 400       IF (IN.GT.INLIM) GO TO 160
 401   520 CONTINUE
 402       TRX = 2.0D0/X
 403       DTM = FNI + IN
 404       TM = (DTM+FNF)*TRX
 405       TA = 0.0D0
 406       TB = TOL
 407       KK = 1
 408   530 CONTINUE
 409 C
 410 C     BACKWARD RECUR UNINDEXED
 411 C
 412       DO 540 I=1,IN
 413         S = TB
 414         TB = TM*TB + TA
 415         TA = S
 416         DTM = DTM - 1.0D0
 417         TM = (DTM+FNF)*TRX
 418   540 CONTINUE
 419 C     NORMALIZATION
 420       IF (KK.NE.1) GO TO 550
 421       TA = (TA/TB)*TEMP(3)
 422       TB = TEMP(3)
 423       KK = 2
 424       IN = NS
 425       IF (NS.NE.0) GO TO 530
 426   550 Y(NN) = TB
 427       NZ = N - NN
 428       IF (NN.EQ.1) RETURN
 429       TB = TM*TB + TA
 430       K = NN - 1
 431       Y(K) = TB
 432       IF (NN.EQ.2) RETURN
 433       DTM = DTM - 1.0D0
 434       TM = (DTM+FNF)*TRX
 435       KM = K - 1
 436 C
 437 C     BACKWARD RECUR INDEXED
 438 C
 439       DO 560 I=1,KM
 440         Y(K-1) = TM*Y(K) + Y(K+1)
 441         DTM = DTM - 1.0D0
 442         TM = (DTM+FNF)*TRX
 443         K = K - 1
 444   560 CONTINUE
 445       RETURN
 446 C
 447 C
 448 C
 449   570 CONTINUE
 450       CALL XERMSG ('SLATEC', 'DBESI',
 451      +   'SCALING OPTION, KODE, NOT 1 OR 2.', 2, 1)
 452       RETURN
 453   580 CONTINUE
 454       CALL XERMSG ('SLATEC', 'DBESI', 'ORDER, ALPHA, LESS THAN ZERO.',
 455      +   2, 1)
 456       RETURN
 457   590 CONTINUE
 458       CALL XERMSG ('SLATEC', 'DBESI', 'N LESS THAN ONE.', 2, 1)
 459       RETURN
 460   600 CONTINUE
 461       CALL XERMSG ('SLATEC', 'DBESI', 'X LESS THAN ZERO.', 2, 1)
 462       RETURN
 463   610 CONTINUE
 464       CALL XERMSG ('SLATEC', 'DBESI',
 465      +   'OVERFLOW, X TOO LARGE FOR KODE = 1.', 6, 1)
 466       RETURN
 467       END