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