*DECK CASYI SUBROUTINE CASYI (Z, FNU, KODE, N, Y, NZ, RL, TOL, ELIM, ALIM) C***BEGIN PROLOGUE CASYI C***SUBSIDIARY C***PURPOSE Subsidiary to CBESI and CBESK C***LIBRARY SLATEC C***TYPE ALL (CASYI-A, ZASYI-A) C***AUTHOR Amos, D. E., (SNL) C***DESCRIPTION C C CASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE ABS(Z) IN THE C REGION ABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN. C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1. C C***SEE ALSO CBESI, CBESK C***ROUTINES CALLED R1MACH C***REVISION HISTORY (YYMMDD) C 830501 DATE WRITTEN C 910415 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE CASYI COMPLEX AK1, CK, CONE, CS1, CS2, CZ, CZERO, DK, EZ, P1, RZ, S2, * Y, Z REAL AA, ACZ, AEZ, AK, ALIM, ARG, ARM, ATOL, AZ, BB, BK, DFNU, * DNU2, ELIM, FDN, FNU, PI, RL, RTPI, RTR1, S, SGN, SQK, TOL, X, * YY, R1MACH INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ DIMENSION Y(N) DATA PI, RTPI /3.14159265358979324E0 , 0.159154943091895336E0 / DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) / C***FIRST EXECUTABLE STATEMENT CASYI NZ = 0 AZ = ABS(Z) X = REAL(Z) ARM = 1.0E+3*R1MACH(1) RTR1 = SQRT(ARM) IL = MIN(2,N) DFNU = FNU + (N-IL) C----------------------------------------------------------------------- C OVERFLOW TEST C----------------------------------------------------------------------- AK1 = CMPLX(RTPI,0.0E0)/Z AK1 = CSQRT(AK1) CZ = Z IF (KODE.EQ.2) CZ = Z - CMPLX(X,0.0E0) ACZ = REAL(CZ) IF (ABS(ACZ).GT.ELIM) GO TO 80 DNU2 = DFNU + DFNU KODED = 1 IF ((ABS(ACZ).GT.ALIM) .AND. (N.GT.2)) GO TO 10 KODED = 0 AK1 = AK1*CEXP(CZ) 10 CONTINUE FDN = 0.0E0 IF (DNU2.GT.RTR1) FDN = DNU2*DNU2 EZ = Z*CMPLX(8.0E0,0.0E0) C----------------------------------------------------------------------- C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE C EXPANSION FOR THE IMAGINARY PART. C----------------------------------------------------------------------- AEZ = 8.0E0*AZ S = TOL/AEZ JL = RL+RL + 2 YY = AIMAG(Z) P1 = CZERO IF (YY.EQ.0.0E0) GO TO 20 C----------------------------------------------------------------------- C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF C SIGNIFICANCE WHEN FNU OR N IS LARGE C----------------------------------------------------------------------- INU = FNU ARG = (FNU-INU)*PI INU = INU + N - IL AK = -SIN(ARG) BK = COS(ARG) IF (YY.LT.0.0E0) BK = -BK P1 = CMPLX(AK,BK) IF (MOD(INU,2).EQ.1) P1 = -P1 20 CONTINUE DO 50 K=1,IL SQK = FDN - 1.0E0 ATOL = S*ABS(SQK) SGN = 1.0E0 CS1 = CONE CS2 = CONE CK = CONE AK = 0.0E0 AA = 1.0E0 BB = AEZ DK = EZ DO 30 J=1,JL CK = CK*CMPLX(SQK,0.0E0)/DK CS2 = CS2 + CK SGN = -SGN CS1 = CS1 + CK*CMPLX(SGN,0.0E0) DK = DK + EZ AA = AA*ABS(SQK)/BB BB = BB + AEZ AK = AK + 8.0E0 SQK = SQK - AK IF (AA.LE.ATOL) GO TO 40 30 CONTINUE GO TO 90 40 CONTINUE S2 = CS1 IF (X+X.LT.ELIM) S2 = S2 + P1*CS2*CEXP(-Z-Z) FDN = FDN + 8.0E0*DFNU + 4.0E0 P1 = -P1 M = N - IL + K Y(M) = S2*AK1 50 CONTINUE IF (N.LE.2) RETURN NN = N K = NN - 2 AK = K RZ = (CONE+CONE)/Z IB = 3 DO 60 I=IB,NN Y(K) = CMPLX(AK+FNU,0.0E0)*RZ*Y(K+1) + Y(K+2) AK = AK - 1.0E0 K = K - 1 60 CONTINUE IF (KODED.EQ.0) RETURN CK = CEXP(CZ) DO 70 I=1,NN Y(I) = Y(I)*CK 70 CONTINUE RETURN 80 CONTINUE NZ = -1 RETURN 90 CONTINUE NZ=-2 RETURN END