github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/internal/amos/amoslib/zmlri.f (about) 1 SUBROUTINE ZMLRI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL) 2 C***BEGIN PROLOGUE ZMLRI 3 C***REFER TO ZBESI,ZBESK 4 C 5 C ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE 6 C MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES. 7 C 8 C***ROUTINES CALLED DGAMLN,D1MACH,ZABS,ZEXP,ZLOG,ZMLT 9 C***END PROLOGUE ZMLRI 10 C COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z 11 DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI, 12 * CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, FNU, PTI, PTR, P1I, 13 * P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI, 14 * SUMR, TFNF, TOL, TST, YI, YR, ZEROI, ZEROR, ZI, ZR, DGAMLN, 15 * D1MACH, ZABS 16 INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ 17 DIMENSION YR(N), YI(N) 18 DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / 19 SCLE = D1MACH(1)/TOL 20 NZ=0 21 AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 22 IAZ = INT(SNGL(AZ)) 23 IFNU = INT(SNGL(FNU)) 24 INU = IFNU + N - 1 25 AT = DBLE(FLOAT(IAZ)) + 1.0D0 26 RAZ = 1.0D0/AZ 27 STR = ZR*RAZ 28 STI = -ZI*RAZ 29 CKR = STR*AT*RAZ 30 CKI = STI*AT*RAZ 31 RZR = (STR+STR)*RAZ 32 RZI = (STI+STI)*RAZ 33 P1R = ZEROR 34 P1I = ZEROI 35 P2R = CONER 36 P2I = CONEI 37 ACK = (AT+1.0D0)*RAZ 38 RHO = ACK + DSQRT(ACK*ACK-1.0D0) 39 RHO2 = RHO*RHO 40 TST = (RHO2+RHO2)/((RHO2-1.0D0)*(RHO-1.0D0)) 41 TST = TST/TOL 42 C----------------------------------------------------------------------- 43 C COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES 44 C----------------------------------------------------------------------- 45 AK = AT 46 DO 10 I=1,80 47 PTR = P2R 48 PTI = P2I 49 P2R = P1R - (CKR*PTR-CKI*PTI) 50 P2I = P1I - (CKI*PTR+CKR*PTI) 51 P1R = PTR 52 P1I = PTI 53 CKR = CKR + RZR 54 CKI = CKI + RZI 55 AP = ZABS(CMPLX(P2R,P2I,kind=KIND(1.0D0))) 56 IF (AP.GT.TST*AK*AK) THEN 57 GO TO 20 58 END IF 59 AK = AK + 1.0D0 60 10 CONTINUE 61 GO TO 110 62 20 CONTINUE 63 I = I + 1 64 K = 0 65 IF (INU.LT.IAZ) GO TO 40 66 C----------------------------------------------------------------------- 67 C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS 68 C----------------------------------------------------------------------- 69 P1R = ZEROR 70 P1I = ZEROI 71 P2R = CONER 72 P2I = CONEI 73 AT = DBLE(FLOAT(INU)) + 1.0D0 74 STR = ZR*RAZ 75 STI = -ZI*RAZ 76 CKR = STR*AT*RAZ 77 CKI = STI*AT*RAZ 78 ACK = AT*RAZ 79 TST = DSQRT(ACK/TOL) 80 ITIME = 1 81 DO 30 K=1,80 82 PTR = P2R 83 PTI = P2I 84 P2R = P1R - (CKR*PTR-CKI*PTI) 85 P2I = P1I - (CKR*PTI+CKI*PTR) 86 P1R = PTR 87 P1I = PTI 88 CKR = CKR + RZR 89 CKI = CKI + RZI 90 AP = ZABS(CMPLX(P2R,P2I,kind=KIND(1.0D0))) 91 IF (AP.LT.TST) GO TO 30 92 IF (ITIME.EQ.2) GO TO 40 93 ACK = ZABS(CMPLX(CKR,CKI,kind=KIND(1.0D0))) 94 FLAM = ACK + DSQRT(ACK*ACK-1.0D0) 95 FKAP = AP/ZABS(CMPLX(P1R,P1I,kind=KIND(1.0D0))) 96 RHO = DMIN1(FLAM,FKAP) 97 TST = TST*DSQRT(RHO/(RHO*RHO-1.0D0)) 98 ITIME = 2 99 30 CONTINUE 100 GO TO 110 101 40 CONTINUE 102 C----------------------------------------------------------------------- 103 C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION 104 C----------------------------------------------------------------------- 105 K = K + 1 106 KK = MAX0(I+IAZ,K+INU) 107 FKK = DBLE(FLOAT(KK)) 108 P1R = ZEROR 109 P1I = ZEROI 110 C----------------------------------------------------------------------- 111 C SCALE P2 AND SUM BY SCLE 112 C----------------------------------------------------------------------- 113 P2R = SCLE 114 P2I = ZEROI 115 FNF = FNU - DBLE(FLOAT(IFNU)) 116 TFNF = FNF + FNF 117 BK = DGAMLN(FKK+TFNF+1.0D0,IDUM) - DGAMLN(FKK+1.0D0,IDUM) - 118 * DGAMLN(TFNF+1.0D0,IDUM) 119 BK = DEXP(BK) 120 SUMR = ZEROR 121 SUMI = ZEROI 122 KM = KK - INU 123 DO 50 I=1,KM 124 PTR = P2R 125 PTI = P2I 126 P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) 127 P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI) 128 P1R = PTR 129 P1I = PTI 130 AK = 1.0D0 - TFNF/(FKK+TFNF) 131 ACK = BK*AK 132 SUMR = SUMR + (ACK+BK)*P1R 133 SUMI = SUMI + (ACK+BK)*P1I 134 BK = ACK 135 FKK = FKK - 1.0D0 136 50 CONTINUE 137 YR(N) = P2R 138 YI(N) = P2I 139 IF (N.EQ.1) GO TO 70 140 DO 60 I=2,N 141 PTR = P2R 142 PTI = P2I 143 P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) 144 P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI) 145 P1R = PTR 146 P1I = PTI 147 AK = 1.0D0 - TFNF/(FKK+TFNF) 148 ACK = BK*AK 149 SUMR = SUMR + (ACK+BK)*P1R 150 SUMI = SUMI + (ACK+BK)*P1I 151 BK = ACK 152 FKK = FKK - 1.0D0 153 M = N - I + 1 154 YR(M) = P2R 155 YI(M) = P2I 156 60 CONTINUE 157 70 CONTINUE 158 IF (IFNU.LE.0) GO TO 90 159 DO 80 I=1,IFNU 160 PTR = P2R 161 PTI = P2I 162 P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) 163 P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR) 164 P1R = PTR 165 P1I = PTI 166 AK = 1.0D0 - TFNF/(FKK+TFNF) 167 ACK = BK*AK 168 SUMR = SUMR + (ACK+BK)*P1R 169 SUMI = SUMI + (ACK+BK)*P1I 170 BK = ACK 171 FKK = FKK - 1.0D0 172 80 CONTINUE 173 90 CONTINUE 174 PTR = ZR 175 PTI = ZI 176 IF (KODE.EQ.2) PTR = ZEROR 177 CALL ZLOG(RZR, RZI, STR, STI, IDUM) 178 P1R = -FNF*STR + PTR 179 P1I = -FNF*STI + PTI 180 AP = DGAMLN(1.0D0+FNF,IDUM) 181 PTR = P1R - AP 182 PTI = P1I 183 C----------------------------------------------------------------------- 184 C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW 185 C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES 186 C----------------------------------------------------------------------- 187 P2R = P2R + SUMR 188 P2I = P2I + SUMI 189 AP = ZABS(CMPLX(P2R,P2I,kind=KIND(1.0D0))) 190 P1R = 1.0D0/AP 191 CALL ZEXP(PTR, PTI, STR, STI) 192 CKR = STR*P1R 193 CKI = STI*P1R 194 PTR = P2R*P1R 195 PTI = -P2I*P1R 196 CALL ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI) 197 DO 100 I=1,N 198 STR = YR(I)*CNORMR - YI(I)*CNORMI 199 YI(I) = YR(I)*CNORMI + YI(I)*CNORMR 200 YR(I) = STR 201 100 CONTINUE 202 RETURN 203 110 CONTINUE 204 NZ=-2 205 RETURN 206 END