1 SUBROUTINE CMLRI(Z, FNU, KODE, N, Y, NZ, TOL) 2C***BEGIN PROLOGUE CMLRI 3C***REFER TO CBESI,CBESK 4C 5C CMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE 6C MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES. 7C 8C***ROUTINES CALLED GAMLN,R1MACH 9C***END PROLOGUE CMLRI 10 COMPLEX CK, CNORM, CONE, CTWO, CZERO, PT, P1, P2, RZ, SUM, Y, Z 11 REAL ACK, AK, AP, AT, AZ, BK, FKAP, FKK, FLAM, FNF, FNU, RHO, 12 * RHO2, SCLE, TFNF, TOL, TST, X, GAMLN, R1MACH 13 INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N 14 DIMENSION Y(N) 15 DATA CZERO,CONE,CTWO /(0.0E0,0.0E0),(1.0E0,0.0E0),(2.0E0,0.0E0)/ 16 SCLE = 1.0E+3*R1MACH(1)/TOL 17 NZ=0 18 AZ = CABS(Z) 19 X = REAL(Z) 20 IAZ = INT(AZ) 21 IFNU = INT(FNU) 22 INU = IFNU + N - 1 23 AT = FLOAT(IAZ) + 1.0E0 24 CK = CMPLX(AT,0.0E0)/Z 25 RZ = CTWO/Z 26 P1 = CZERO 27 P2 = CONE 28 ACK = (AT+1.0E0)/AZ 29 RHO = ACK + SQRT(ACK*ACK-1.0E0) 30 RHO2 = RHO*RHO 31 TST = (RHO2+RHO2)/((RHO2-1.0E0)*(RHO-1.0E0)) 32 TST = TST/TOL 33C----------------------------------------------------------------------- 34C COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES 35C----------------------------------------------------------------------- 36 AK = AT 37 DO 10 I=1,80 38 PT = P2 39 P2 = P1 - CK*P2 40 P1 = PT 41 CK = CK + RZ 42 AP = CABS(P2) 43 IF (AP.GT.TST*AK*AK) GO TO 20 44 AK = AK + 1.0E0 45 10 CONTINUE 46 GO TO 110 47 20 CONTINUE 48 I = I + 1 49 K = 0 50 IF (INU.LT.IAZ) GO TO 40 51C----------------------------------------------------------------------- 52C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS 53C----------------------------------------------------------------------- 54 P1 = CZERO 55 P2 = CONE 56 AT = FLOAT(INU) + 1.0E0 57 CK = CMPLX(AT,0.0E0)/Z 58 ACK = AT/AZ 59 TST = SQRT(ACK/TOL) 60 ITIME = 1 61 DO 30 K=1,80 62 PT = P2 63 P2 = P1 - CK*P2 64 P1 = PT 65 CK = CK + RZ 66 AP = CABS(P2) 67 IF (AP.LT.TST) GO TO 30 68 IF (ITIME.EQ.2) GO TO 40 69 ACK = CABS(CK) 70 FLAM = ACK + SQRT(ACK*ACK-1.0E0) 71 FKAP = AP/CABS(P1) 72 RHO = AMIN1(FLAM,FKAP) 73 TST = TST*SQRT(RHO/(RHO*RHO-1.0E0)) 74 ITIME = 2 75 30 CONTINUE 76 GO TO 110 77 40 CONTINUE 78C----------------------------------------------------------------------- 79C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION 80C----------------------------------------------------------------------- 81 K = K + 1 82 KK = MAX0(I+IAZ,K+INU) 83 FKK = FLOAT(KK) 84 P1 = CZERO 85C----------------------------------------------------------------------- 86C SCALE P2 AND SUM BY SCLE 87C----------------------------------------------------------------------- 88 P2 = CMPLX(SCLE,0.0E0) 89 FNF = FNU - FLOAT(IFNU) 90 TFNF = FNF + FNF 91 BK = GAMLN(FKK+TFNF+1.0E0,IDUM) - GAMLN(FKK+1.0E0,IDUM) 92 * -GAMLN(TFNF+1.0E0,IDUM) 93 BK = EXP(BK) 94 SUM = CZERO 95 KM = KK - INU 96 DO 50 I=1,KM 97 PT = P2 98 P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2 99 P1 = PT 100 AK = 1.0E0 - TFNF/(FKK+TFNF) 101 ACK = BK*AK 102 SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1 103 BK = ACK 104 FKK = FKK - 1.0E0 105 50 CONTINUE 106 Y(N) = P2 107 IF (N.EQ.1) GO TO 70 108 DO 60 I=2,N 109 PT = P2 110 P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2 111 P1 = PT 112 AK = 1.0E0 - TFNF/(FKK+TFNF) 113 ACK = BK*AK 114 SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1 115 BK = ACK 116 FKK = FKK - 1.0E0 117 M = N - I + 1 118 Y(M) = P2 119 60 CONTINUE 120 70 CONTINUE 121 IF (IFNU.LE.0) GO TO 90 122 DO 80 I=1,IFNU 123 PT = P2 124 P2 = P1 + CMPLX(FKK+FNF,0.0E0)*RZ*P2 125 P1 = PT 126 AK = 1.0E0 - TFNF/(FKK+TFNF) 127 ACK = BK*AK 128 SUM = SUM + CMPLX(ACK+BK,0.0E0)*P1 129 BK = ACK 130 FKK = FKK - 1.0E0 131 80 CONTINUE 132 90 CONTINUE 133 PT = Z 134 IF (KODE.EQ.2) PT = PT - CMPLX(X,0.0E0) 135 P1 = -CMPLX(FNF,0.0E0)*CLOG(RZ) + PT 136 AP = GAMLN(1.0E0+FNF,IDUM) 137 PT = P1 - CMPLX(AP,0.0E0) 138C----------------------------------------------------------------------- 139C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW 140C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES 141C----------------------------------------------------------------------- 142 P2 = P2 + SUM 143 AP = CABS(P2) 144 P1 = CMPLX(1.0E0/AP,0.0E0) 145 CK = CEXP(PT)*P1 146 PT = CONJG(P2)*P1 147 CNORM = CK*PT 148 DO 100 I=1,N 149 Y(I) = Y(I)*CNORM 150 100 CONTINUE 151 RETURN 152 110 CONTINUE 153 NZ=-2 154 RETURN 155 END 156