1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
|
*DECK DASYIK
SUBROUTINE DASYIK (X, FNU, KODE, FLGIK, RA, ARG, IN, Y)
C***BEGIN PROLOGUE DASYIK
C***SUBSIDIARY
C***PURPOSE Subsidiary to DBESI and DBESK
C***LIBRARY SLATEC
C***TYPE DOUBLE PRECISION (ASYIK-S, DASYIK-D)
C***AUTHOR Amos, D. E., (SNLA)
C***DESCRIPTION
C
C DASYIK computes Bessel functions I and K
C for arguments X.GT.0.0 and orders FNU.GE.35
C on FLGIK = 1 and FLGIK = -1 respectively.
C
C INPUT
C
C X - Argument, X.GT.0.0D0
C FNU - Order of first Bessel function
C KODE - A parameter to indicate the scaling option
C KODE=1 returns Y(I)= I/SUB(FNU+I-1)/(X), I=1,IN
C or Y(I)= K/SUB(FNU+I-1)/(X), I=1,IN
C on FLGIK = 1.0D0 or FLGIK = -1.0D0
C KODE=2 returns Y(I)=EXP(-X)*I/SUB(FNU+I-1)/(X), I=1,IN
C or Y(I)=EXP( X)*K/SUB(FNU+I-1)/(X), I=1,IN
C on FLGIK = 1.0D0 or FLGIK = -1.0D0
C FLGIK - Selection parameter for I or K FUNCTION
C FLGIK = 1.0D0 gives the I function
C FLGIK = -1.0D0 gives the K function
C RA - SQRT(1.+Z*Z), Z=X/FNU
C ARG - Argument of the leading exponential
C IN - Number of functions desired, IN=1 or 2
C
C OUTPUT
C
C Y - A vector whose first IN components contain the sequence
C
C Abstract **** A double precision routine ****
C DASYIK implements the uniform asymptotic expansion of
C the I and K Bessel functions for FNU.GE.35 and real
C X.GT.0.0D0. The forms are identical except for a change
C in sign of some of the terms. This change in sign is
C accomplished by means of the FLAG FLGIK = 1 or -1.
C
C***SEE ALSO DBESI, DBESK
C***ROUTINES CALLED D1MACH
C***REVISION HISTORY (YYMMDD)
C 750101 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890911 Removed unnecessary intrinsics. (WRB)
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900328 Added TYPE section. (WRB)
C 910408 Updated the AUTHOR section. (WRB)
C***END PROLOGUE DASYIK
C
INTEGER IN, J, JN, K, KK, KODE, L
DOUBLE PRECISION AK,AP,ARG,C,COEF,CON,ETX,FLGIK,FN,FNU,GLN,RA,
1 S1, S2, T, TOL, T2, X, Y, Z
DOUBLE PRECISION D1MACH
DIMENSION Y(*), C(65), CON(2)
SAVE CON, C
DATA CON(1), CON(2) /
1 3.98942280401432678D-01, 1.25331413731550025D+00/
DATA C(1), C(2), C(3), C(4), C(5), C(6), C(7), C(8), C(9), C(10),
1 C(11), C(12), C(13), C(14), C(15), C(16), C(17), C(18),
2 C(19), C(20), C(21), C(22), C(23), C(24)/
3 -2.08333333333333D-01, 1.25000000000000D-01,
4 3.34201388888889D-01, -4.01041666666667D-01,
5 7.03125000000000D-02, -1.02581259645062D+00,
6 1.84646267361111D+00, -8.91210937500000D-01,
7 7.32421875000000D-02, 4.66958442342625D+00,
8 -1.12070026162230D+01, 8.78912353515625D+00,
9 -2.36408691406250D+00, 1.12152099609375D-01,
1 -2.82120725582002D+01, 8.46362176746007D+01,
2 -9.18182415432400D+01, 4.25349987453885D+01,
3 -7.36879435947963D+00, 2.27108001708984D-01,
4 2.12570130039217D+02, -7.65252468141182D+02,
5 1.05999045252800D+03, -6.99579627376133D+02/
DATA C(25), C(26), C(27), C(28), C(29), C(30), C(31), C(32),
1 C(33), C(34), C(35), C(36), C(37), C(38), C(39), C(40),
2 C(41), C(42), C(43), C(44), C(45), C(46), C(47), C(48)/
3 2.18190511744212D+02, -2.64914304869516D+01,
4 5.72501420974731D-01, -1.91945766231841D+03,
5 8.06172218173731D+03, -1.35865500064341D+04,
6 1.16553933368645D+04, -5.30564697861340D+03,
7 1.20090291321635D+03, -1.08090919788395D+02,
8 1.72772750258446D+00, 2.02042913309661D+04,
9 -9.69805983886375D+04, 1.92547001232532D+05,
1 -2.03400177280416D+05, 1.22200464983017D+05,
2 -4.11926549688976D+04, 7.10951430248936D+03,
3 -4.93915304773088D+02, 6.07404200127348D+00,
4 -2.42919187900551D+05, 1.31176361466298D+06,
5 -2.99801591853811D+06, 3.76327129765640D+06/
DATA C(49), C(50), C(51), C(52), C(53), C(54), C(55), C(56),
1 C(57), C(58), C(59), C(60), C(61), C(62), C(63), C(64),
2 C(65)/
3 -2.81356322658653D+06, 1.26836527332162D+06,
4 -3.31645172484564D+05, 4.52187689813627D+04,
5 -2.49983048181121D+03, 2.43805296995561D+01,
6 3.28446985307204D+06, -1.97068191184322D+07,
7 5.09526024926646D+07, -7.41051482115327D+07,
8 6.63445122747290D+07, -3.75671766607634D+07,
9 1.32887671664218D+07, -2.78561812808645D+06,
1 3.08186404612662D+05, -1.38860897537170D+04,
2 1.10017140269247D+02/
C***FIRST EXECUTABLE STATEMENT DASYIK
TOL = D1MACH(3)
TOL = MAX(TOL,1.0D-15)
FN = FNU
Z = (3.0D0-FLGIK)/2.0D0
KK = INT(Z)
DO 50 JN=1,IN
IF (JN.EQ.1) GO TO 10
FN = FN - FLGIK
Z = X/FN
RA = SQRT(1.0D0+Z*Z)
GLN = LOG((1.0D0+RA)/Z)
ETX = KODE - 1
T = RA*(1.0D0-ETX) + ETX/(Z+RA)
ARG = FN*(T-GLN)*FLGIK
10 COEF = EXP(ARG)
T = 1.0D0/RA
T2 = T*T
T = T/FN
T = SIGN(T,FLGIK)
S2 = 1.0D0
AP = 1.0D0
L = 0
DO 30 K=2,11
L = L + 1
S1 = C(L)
DO 20 J=2,K
L = L + 1
S1 = S1*T2 + C(L)
20 CONTINUE
AP = AP*T
AK = AP*S1
S2 = S2 + AK
IF (MAX(ABS(AK),ABS(AP)) .LT.TOL) GO TO 40
30 CONTINUE
40 CONTINUE
T = ABS(T)
Y(JN) = S2*COEF*SQRT(T)*CON(KK)
50 CONTINUE
RETURN
END
|