src/numerical/slatec/fortran/dasyik.f

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