1*DECK ZASYI
2      SUBROUTINE ZASYI (ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM,
3     +   ALIM)
4C***BEGIN PROLOGUE  ZASYI
5C***SUBSIDIARY
6C***PURPOSE  Subsidiary to ZBESI and ZBESK
7C***LIBRARY   SLATEC
8C***TYPE      ALL (CASYI-A, ZASYI-A)
9C***AUTHOR  Amos, D. E., (SNL)
10C***DESCRIPTION
11C
12C     ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
13C     MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE ABS(Z) IN THE
14C     REGION ABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
15C     NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
16C
17C***SEE ALSO  ZBESI, ZBESK
18C***ROUTINES CALLED  D1MACH, ZABS, ZDIV, ZEXP, ZMLT, ZSQRT
19C***REVISION HISTORY  (YYMMDD)
20C   830501  DATE WRITTEN
21C   910415  Prologue converted to Version 4.0 format.  (BAB)
22C   930122  Added ZEXP and ZSQRT to EXTERNAL statement.  (RWC)
23C***END PROLOGUE  ZASYI
24C     COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z
25      DOUBLE PRECISION AA, AEZ, AK, AK1I, AK1R, ALIM, ARG, ARM, ATOL,
26     * AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI,
27     * CZR, DFNU, DKI, DKR, DNU2, ELIM, EZI, EZR, FDN, FNU, PI, P1I,
28     * P1R, RAZ, RL, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I,
29     * S2R, TOL, TZI, TZR, YI, YR, ZEROI, ZEROR, ZI, ZR, D1MACH, ZABS
30      INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ
31      DIMENSION YR(N), YI(N)
32      EXTERNAL ZABS, ZEXP, ZSQRT
33      DATA PI, RTPI  /3.14159265358979324D0 , 0.159154943091895336D0 /
34      DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
35C***FIRST EXECUTABLE STATEMENT  ZASYI
36      NZ = 0
37      AZ = ZABS(ZR,ZI)
38      ARM = 1.0D+3*D1MACH(1)
39      RTR1 = SQRT(ARM)
40      IL = MIN(2,N)
41      DFNU = FNU + (N-IL)
42C-----------------------------------------------------------------------
43C     OVERFLOW TEST
44C-----------------------------------------------------------------------
45      RAZ = 1.0D0/AZ
46      STR = ZR*RAZ
47      STI = -ZI*RAZ
48      AK1R = RTPI*STR*RAZ
49      AK1I = RTPI*STI*RAZ
50      CALL ZSQRT(AK1R, AK1I, AK1R, AK1I)
51      CZR = ZR
52      CZI = ZI
53      IF (KODE.NE.2) GO TO 10
54      CZR = ZEROR
55      CZI = ZI
56   10 CONTINUE
57      IF (ABS(CZR).GT.ELIM) GO TO 100
58      DNU2 = DFNU + DFNU
59      KODED = 1
60      IF ((ABS(CZR).GT.ALIM) .AND. (N.GT.2)) GO TO 20
61      KODED = 0
62      CALL ZEXP(CZR, CZI, STR, STI)
63      CALL ZMLT(AK1R, AK1I, STR, STI, AK1R, AK1I)
64   20 CONTINUE
65      FDN = 0.0D0
66      IF (DNU2.GT.RTR1) FDN = DNU2*DNU2
67      EZR = ZR*8.0D0
68      EZI = ZI*8.0D0
69C-----------------------------------------------------------------------
70C     WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
71C     FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
72C     EXPANSION FOR THE IMAGINARY PART.
73C-----------------------------------------------------------------------
74      AEZ = 8.0D0*AZ
75      S = TOL/AEZ
76      JL = RL+RL + 2
77      P1R = ZEROR
78      P1I = ZEROI
79      IF (ZI.EQ.0.0D0) GO TO 30
80C-----------------------------------------------------------------------
81C     CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
82C     SIGNIFICANCE WHEN FNU OR N IS LARGE
83C-----------------------------------------------------------------------
84      INU = FNU
85      ARG = (FNU-INU)*PI
86      INU = INU + N - IL
87      AK = -SIN(ARG)
88      BK = COS(ARG)
89      IF (ZI.LT.0.0D0) BK = -BK
90      P1R = AK
91      P1I = BK
92      IF (MOD(INU,2).EQ.0) GO TO 30
93      P1R = -P1R
94      P1I = -P1I
95   30 CONTINUE
96      DO 70 K=1,IL
97        SQK = FDN - 1.0D0
98        ATOL = S*ABS(SQK)
99        SGN = 1.0D0
100        CS1R = CONER
101        CS1I = CONEI
102        CS2R = CONER
103        CS2I = CONEI
104        CKR = CONER
105        CKI = CONEI
106        AK = 0.0D0
107        AA = 1.0D0
108        BB = AEZ
109        DKR = EZR
110        DKI = EZI
111        DO 40 J=1,JL
112          CALL ZDIV(CKR, CKI, DKR, DKI, STR, STI)
113          CKR = STR*SQK
114          CKI = STI*SQK
115          CS2R = CS2R + CKR
116          CS2I = CS2I + CKI
117          SGN = -SGN
118          CS1R = CS1R + CKR*SGN
119          CS1I = CS1I + CKI*SGN
120          DKR = DKR + EZR
121          DKI = DKI + EZI
122          AA = AA*ABS(SQK)/BB
123          BB = BB + AEZ
124          AK = AK + 8.0D0
125          SQK = SQK - AK
126          IF (AA.LE.ATOL) GO TO 50
127   40   CONTINUE
128        GO TO 110
129   50   CONTINUE
130        S2R = CS1R
131        S2I = CS1I
132        IF (ZR+ZR.GE.ELIM) GO TO 60
133        TZR = ZR + ZR
134        TZI = ZI + ZI
135        CALL ZEXP(-TZR, -TZI, STR, STI)
136        CALL ZMLT(STR, STI, P1R, P1I, STR, STI)
137        CALL ZMLT(STR, STI, CS2R, CS2I, STR, STI)
138        S2R = S2R + STR
139        S2I = S2I + STI
140   60   CONTINUE
141        FDN = FDN + 8.0D0*DFNU + 4.0D0
142        P1R = -P1R
143        P1I = -P1I
144        M = N - IL + K
145        YR(M) = S2R*AK1R - S2I*AK1I
146        YI(M) = S2R*AK1I + S2I*AK1R
147   70 CONTINUE
148      IF (N.LE.2) RETURN
149      NN = N
150      K = NN - 2
151      AK = K
152      STR = ZR*RAZ
153      STI = -ZI*RAZ
154      RZR = (STR+STR)*RAZ
155      RZI = (STI+STI)*RAZ
156      IB = 3
157      DO 80 I=IB,NN
158        YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2)
159        YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2)
160        AK = AK - 1.0D0
161        K = K - 1
162   80 CONTINUE
163      IF (KODED.EQ.0) RETURN
164      CALL ZEXP(CZR, CZI, CKR, CKI)
165      DO 90 I=1,NN
166        STR = YR(I)*CKR - YI(I)*CKI
167        YI(I) = YR(I)*CKI + YI(I)*CKR
168        YR(I) = STR
169   90 CONTINUE
170      RETURN
171  100 CONTINUE
172      NZ = -1
173      RETURN
174  110 CONTINUE
175      NZ=-2
176      RETURN
177      END
178