2
SUBROUTINE DBESJ (X, ALPHA, N, Y, NZ)
3
C***BEGIN PROLOGUE DBESJ
4
C***PURPOSE Compute an N member sequence of J Bessel functions
5
C J/SUB(ALPHA+K-1)/(X), K=1,...,N for non-negative ALPHA
9
C***TYPE DOUBLE PRECISION (BESJ-S, DBESJ-D)
10
C***KEYWORDS J BESSEL FUNCTION, SPECIAL FUNCTIONS
11
C***AUTHOR Amos, D. E., (SNLA)
12
C Daniel, S. L., (SNLA)
13
C Weston, M. K., (SNLA)
16
C Abstract **** a double precision routine ****
17
C DBESJ computes an N member sequence of J Bessel functions
18
C J/sub(ALPHA+K-1)/(X), K=1,...,N for non-negative ALPHA and X.
19
C A combination of the power series, the asymptotic expansion
20
C for X to infinity and the uniform asymptotic expansion for
21
C NU to infinity are applied over subdivisions of the (NU,X)
22
C plane. For values of (NU,X) not covered by one of these
23
C formulae, the order is incremented or decremented by integer
24
C values into a region where one of the formulae apply. Backward
25
C recursion is applied to reduce orders by integer values except
26
C where the entire sequence lies in the oscillatory region. In
27
C this case forward recursion is stable and values from the
28
C asymptotic expansion for X to infinity start the recursion
29
C when it is efficient to do so. Leading terms of the series and
30
C uniform expansion are tested for underflow. 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 members are set to zero.
34
C Overflow cannot occur.
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.
40
C Description of Arguments
42
C Input X,ALPHA are double precision
44
C ALPHA - order of first member of the sequence,
46
C N - number of members in the sequence, N .GE. 1
48
C Output Y is double precision
49
C Y - a vector whose first N components contain
50
C values for J/sub(ALPHA+K-1)/(X), K=1,...,N
51
C NZ - number of components of Y set to zero due to
53
C NZ=0 , normal return, computation completed
54
C NZ .NE. 0, last NZ components of Y set to zero,
55
C Y(K)=0.0D0, K=N-NZ+1,...,N.
58
C Improper input arguments - a fatal error
59
C Underflow - a non-fatal error (NZ .NE. 0)
61
C***REFERENCES D. E. Amos, S. L. Daniel and M. K. Weston, CDC 6600
62
C subroutines IBESS and JBESS for Bessel functions
63
C I(NU,X) and J(NU,X), X .GE. 0, NU .GE. 0, ACM
64
C Transactions on Mathematical Software 3, (1977),
66
C F. W. J. Olver, Tables of Bessel Functions of Moderate
67
C or Large Orders, NPL Mathematical Tables 6, Her
68
C Majesty's Stationery Office, London, 1962.
69
C***ROUTINES CALLED D1MACH, DASYJY, DJAIRY, DLNGAM, I1MACH, XERMSG
70
C***REVISION HISTORY (YYMMDD)
72
C 890531 Changed all specific intrinsics to generic. (WRB)
73
C 890911 Removed unnecessary intrinsics. (WRB)
74
C 890911 REVISION DATE from Version 3.2
75
C 891214 Prologue converted to Version 4.0 format. (BAB)
76
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
77
C 900326 Removed duplicate information from DESCRIPTION section.
79
C 920501 Reformatted the REFERENCES section. (WRB)
80
C***END PROLOGUE DBESJ
82
INTEGER I,IALP,IDALP,IFLW,IN,INLIM,IS,I1,I2,K,KK,KM,KT,N,NN,
85
DOUBLE PRECISION AK,AKM,ALPHA,ANS,AP,ARG,COEF,DALPHA,DFN,DTM,
86
1 EARG,ELIM1,ETX,FIDAL,FLGJY,FN,FNF,FNI,FNP1,FNU,
87
2 FNULIM,GLN,PDF,PIDT,PP,RDEN,RELB,RTTP,RTWO,RTX,RZDEN,
88
3 S,SA,SB,SXO2,S1,S2,T,TA,TAU,TB,TEMP,TFN,TM,TOL,
89
4 TOLLN,TRX,TX,T1,T2,WK,X,XO2,XO2L,Y,SLIM,RTOL
90
SAVE RTWO, PDF, RTTP, PIDT, PP, INLIM, FNULIM
91
DOUBLE PRECISION D1MACH, DLNGAM
92
DIMENSION Y(*), TEMP(3), FNULIM(2), PP(4), WK(7)
93
DATA RTWO,PDF,RTTP,PIDT / 1.34839972492648D+00,
94
1 7.85398163397448D-01, 7.97884560802865D-01, 1.57079632679490D+00/
95
DATA PP(1), PP(2), PP(3), PP(4) / 8.72909153935547D+00,
96
1 2.65693932265030D-01, 1.24578576865586D-01, 7.70133747430388D-04/
98
DATA FNULIM(1), FNULIM(2) / 100.0D0, 60.0D0 /
99
C***FIRST EXECUTABLE STATEMENT DBESJ
103
C I1MACH(14) REPLACES I1MACH(11) IN A DOUBLE PRECISION CODE
104
C I1MACH(15) REPLACES I1MACH(12) IN A DOUBLE PRECISION CODE
106
TOL = MAX(TA,1.0D-15)
110
ELIM1 = -2.303D0*(I2*TB+3.0D0)
112
SLIM=D1MACH(1)*RTOL*1.0D+3
114
TOLLN = 2.303D0*TB*I1
115
TOLLN = MIN(TOLLN,34.5388D0)
120
30 IF (ALPHA) 710, 40, 50
131
IF (ALPHA.LT.0.0D0) GO TO 710
141
C DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X
142
C TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE
145
IF (SXO2.LE.(FNU+1.0D0)) GO TO 90
147
IF (X.GT.TA) GO TO 120
148
IF (X.GT.12.0D0) GO TO 110
150
NS = INT(SXO2-FNU) + 1
156
IF (X.LE.0.50D0) GO TO 330
163
IF (N-1+NS.GT.0) IS = 3
165
110 ANS = MAX(36.0D0-FNU,0.0D0)
171
IF (N-1+NS.GT.0) IS = 3
176
TA = TAU + FNULIM(KT)
177
IF (FNU.LE.TA) GO TO 480
181
C UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY
187
CALL DASYJY(DJAIRY,X,FN,FLGJY,I1,TEMP(IS),WK,IFLW)
188
IF(IFLW.NE.0) GO TO 380
189
GO TO (320, 450, 620), IS
190
310 TEMP(1) = TEMP(3)
196
IF(I1.EQ.2) GO TO 450
199
C SERIES FOR (X/2)**2.LE.NU+1
204
IF (ARG.LT.(-ELIM1)) GO TO 400
208
IF (X.LT.TOL) GO TO 360
217
IF (ABS(T).LT.TOL) GO TO 360
224
GO TO (370, 450, 610), IS
225
370 EARG = EARG*FN/XO2
232
C SET UNDERFLOW VALUE AND UPDATE PARAMETERS
233
C UNDERFLOW CAN ONLY OCCUR FOR NS=0 SINCE THE ORDER MUST BE LARGER
234
C THAN 36. THEREFORE, NS NEE NOT BE TESTED.
241
IF (NN-1) 440, 390, 130
251
IF (NN-1) 440, 410, 420
254
420 IF (SXO2.LE.FNP1) GO TO 430
256
430 ARG = ARG - XO2L + LOG(FNP1)
257
IF (ARG.LT.(-ELIM1)) GO TO 400
262
C BACKWARD RECURSION SECTION
265
IF(NS.NE.0) GO TO 451
267
IF (KT.EQ.2) GO TO 470
268
C BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
279
IF(ABS(TA).GT.SLIM) GO TO 455
286
IF(IN.EQ.0) GO TO 690
287
IF(NS.NE.0) GO TO 670
302
C ASYMPTOTIC EXPANSION FOR X TO INFINITY WITH FORWARD RECURSION IN
303
C OSCILLATORY REGION X.GT.MAX(20, NU), PROVIDED THE LAST MEMBER
304
C OF THE SEQUENCE IS ALSO IN THE REGION.
307
IN = INT(ALPHA-TAU+2.0D0)
308
IF (IN.LE.0) GO TO 490
309
IDALP = IALP - IN - 1
318
ARG = X - PIDT*DALPHA - PDF
327
IF (FIDAL.EQ.0.0D0 .AND. ABS(FNF).LT.TOL) GO TO 520
328
TM = 4.0D0*FNF*(FIDAL+FIDAL+FNF)
352
IF (ABS(T2).LE.RELB) GO TO 540
355
540 TEMP(IS) = COEF*(S1*SB-S2*SA)
356
IF(IS.EQ.2) GO TO 560
357
FIDAL = FIDAL + 1.0D0
365
C FORWARD RECURSION SECTION
367
560 IF (KT.EQ.2) GO TO 470
372
IF (IN.EQ.0) GO TO 580
374
C FORWARD RECUR TO INDEX ALPHA
382
IF (NN.EQ.1) GO TO 600
389
C FORWARD RECUR FROM INDEX ALPHA TO ALPHA+N-1
395
Y(I) = TM*Y(I-1) - Y(I-2)
402
C BACKWARD RECURSION WITH NORMALIZATION BY
403
C ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES.
406
C COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION
407
AKM = MAX(3.0D0-FN,0.0D0)
410
TA = (GLN+TFN-0.9189385332D0-0.0833333333D0/TFN)/(TFN+0.5D0)
412
TB = -(1.0D0-1.5D0/TFN)/TFN
413
AKM = TOLLN/(-TA+SQRT(TA*TA-TOLLN*TB)) + 1.5D0
417
C COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION
419
IF (WK(6).GT.30.0D0) GO TO 640
420
RDEN = (PP(4)*WK(6)+PP(3))*WK(6) + 1.0D0
421
RZDEN = PP(1) + PP(2)*WK(6)
423
IF (WK(1).LT.0.10D0) GO TO 630
426
630 TB=(1.259921049D0+(0.1679894730D0+0.0887944358D0*WK(1))*WK(1))
430
TA = 0.5D0*TOLLN/WK(4)
431
TA=((0.0493827160D0*TA-0.1111111111D0)*TA+0.6666666667D0)*TA*WK(6)
432
IF (WK(1).LT.0.10D0) GO TO 630
434
650 IN = INT(TA/TB+1.5D0)
435
IF (IN.GT.INLIM) GO TO 310
446
C BACKWARD RECUR UNINDEXED
456
IF (KK.NE.1) GO TO 690
461
IF(ABS(S).GT.SLIM) GO TO 685
469
IF (NS.NE.0) GO TO 670
483
C BACKWARD RECUR INDEXED
499
CALL XERMSG ('SLATEC', 'DBESJ', 'ORDER, ALPHA, LESS THAN ZERO.',
503
CALL XERMSG ('SLATEC', 'DBESJ', 'N LESS THAN ONE.', 2, 1)
506
CALL XERMSG ('SLATEC', 'DBESJ', 'X LESS THAN ZERO.', 2, 1)
added
removed
deps/openlibm/slatec/dbesj.f
