1      SUBROUTINE CASYI(Z, FNU, KODE, N, Y, NZ, RL, TOL, ELIM, ALIM)
2C***BEGIN PROLOGUE  CASYI
3C***REFER TO  CBESI,CBESK
4C
5C     CASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
6C     MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
7C     REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
8C     NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
9C
10C***ROUTINES CALLED  R1MACH
11C***END PROLOGUE  CASYI
12      COMPLEX AK1, CK, CONE, CS1, CS2, CZ, CZERO, DK, EZ, P1, RZ, S2,
13     * Y, Z
14      REAL AA, ACZ, AEZ, AK, ALIM, ARG, ARM, ATOL, AZ, BB, BK, DFNU,
15     * DNU2, ELIM, FDN, FNU, PI, RL, RTPI, RTR1, S, SGN, SQK, TOL, X,
16     * YY, R1MACH
17      INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ
18      DIMENSION Y(N)
19      DATA PI, RTPI  /3.14159265358979324E0 , 0.159154943091895336E0 /
20      DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) /
21C
22      NZ = 0
23      AZ = CABS(Z)
24      X = REAL(Z)
25      ARM = 1.0E+3*R1MACH(1)
26      RTR1 = SQRT(ARM)
27      IL = MIN0(2,N)
28      DFNU = FNU + FLOAT(N-IL)
29C-----------------------------------------------------------------------
30C     OVERFLOW TEST
31C-----------------------------------------------------------------------
32      AK1 = CMPLX(RTPI,0.0E0)/Z
33      AK1 = CSQRT(AK1)
34      CZ = Z
35      IF (KODE.EQ.2) CZ = Z - CMPLX(X,0.0E0)
36      ACZ = REAL(CZ)
37      IF (ABS(ACZ).GT.ELIM) GO TO 80
38      DNU2 = DFNU + DFNU
39      KODED = 1
40      IF ((ABS(ACZ).GT.ALIM) .AND. (N.GT.2)) GO TO 10
41      KODED = 0
42      AK1 = AK1*CEXP(CZ)
43   10 CONTINUE
44      FDN = 0.0E0
45      IF (DNU2.GT.RTR1) FDN = DNU2*DNU2
46      EZ = Z*CMPLX(8.0E0,0.0E0)
47C-----------------------------------------------------------------------
48C     WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
49C     FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
50C     EXPANSION FOR THE IMAGINARY PART.
51C-----------------------------------------------------------------------
52      AEZ = 8.0E0*AZ
53      S = TOL/AEZ
54      JL = INT(RL+RL) + 2
55      YY = AIMAG(Z)
56      P1 = CZERO
57      IF (YY.EQ.0.0E0) GO TO 20
58C-----------------------------------------------------------------------
59C     CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
60C     SIGNIFICANCE WHEN FNU OR N IS LARGE
61C-----------------------------------------------------------------------
62      INU = INT(FNU)
63      ARG = (FNU-FLOAT(INU))*PI
64      INU = INU + N - IL
65      AK = -SIN(ARG)
66      BK = COS(ARG)
67      IF (YY.LT.0.0E0) BK = -BK
68      P1 = CMPLX(AK,BK)
69      IF (MOD(INU,2).EQ.1) P1 = -P1
70   20 CONTINUE
71      DO 50 K=1,IL
72        SQK = FDN - 1.0E0
73        ATOL = S*ABS(SQK)
74        SGN = 1.0E0
75        CS1 = CONE
76        CS2 = CONE
77        CK = CONE
78        AK = 0.0E0
79        AA = 1.0E0
80        BB = AEZ
81        DK = EZ
82        DO 30 J=1,JL
83          CK = CK*CMPLX(SQK,0.0E0)/DK
84          CS2 = CS2 + CK
85          SGN = -SGN
86          CS1 = CS1 + CK*CMPLX(SGN,0.0E0)
87          DK = DK + EZ
88          AA = AA*ABS(SQK)/BB
89          BB = BB + AEZ
90          AK = AK + 8.0E0
91          SQK = SQK - AK
92          IF (AA.LE.ATOL) GO TO 40
93   30   CONTINUE
94        GO TO 90
95   40   CONTINUE
96        S2 = CS1
97        IF (X+X.LT.ELIM) S2 = S2 + P1*CS2*CEXP(-Z-Z)
98        FDN = FDN + 8.0E0*DFNU + 4.0E0
99        P1 = -P1
100        M = N - IL + K
101        Y(M) = S2*AK1
102   50 CONTINUE
103      IF (N.LE.2) RETURN
104      NN = N
105      K = NN - 2
106      AK = FLOAT(K)
107      RZ = (CONE+CONE)/Z
108      IB = 3
109      DO 60 I=IB,NN
110        Y(K) = CMPLX(AK+FNU,0.0E0)*RZ*Y(K+1) + Y(K+2)
111        AK = AK - 1.0E0
112        K = K - 1
113   60 CONTINUE
114      IF (KODED.EQ.0) RETURN
115      CK = CEXP(CZ)
116      DO 70 I=1,NN
117        Y(I) = Y(I)*CK
118   70 CONTINUE
119      RETURN
120   80 CONTINUE
121      NZ = -1
122      RETURN
123   90 CONTINUE
124      NZ=-2
125      RETURN
126      END
127