github.com/gopherd/gonum@v0.0.4/mathext/internal/amos/amoslib/zrati.f (about) 1 SUBROUTINE ZRATI(ZR, ZI, FNU, N, CYR, CYI, TOL) 2 C***BEGIN PROLOGUE ZRATI 3 C***REFER TO ZBESI,ZBESK,ZBESH 4 C 5 C ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD 6 C RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD 7 C RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B, 8 C MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973, 9 C BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER, 10 C BY D. J. SOOKNE. 11 C 12 C***ROUTINES CALLED ZABS,ZDIV 13 C***END PROLOGUE ZRATI 14 C COMPLEX Z,CY(1),CONE,CZERO,P1,P2,T1,RZ,PT,CDFNU 15 DOUBLE PRECISION AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR, 16 * CONEI, CONER, CYI, CYR, CZEROI, CZEROR, DFNU, FDNU, FLAM, FNU, 17 * FNUP, PTI, PTR, P1I, P1R, P2I, P2R, RAK, RAP1, RHO, RT2, RZI, 18 * RZR, TEST, TEST1, TOL, TTI, TTR, T1I, T1R, ZI, ZR, ZABS 19 INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N 20 DIMENSION CYR(N), CYI(N) 21 DATA CZEROR,CZEROI,CONER,CONEI,RT2/ 22 1 0.0D0, 0.0D0, 1.0D0, 0.0D0, 1.41421356237309505D0 / 23 AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 24 INU = INT(SNGL(FNU)) 25 IDNU = INU + N - 1 26 MAGZ = INT(SNGL(AZ)) 27 AMAGZ = DBLE(FLOAT(MAGZ+1)) 28 FDNU = DBLE(FLOAT(IDNU)) 29 FNUP = DMAX1(AMAGZ,FDNU) 30 ID = IDNU - MAGZ - 1 31 ITIME = 1 32 K = 1 33 PTR = 1.0D0/AZ 34 RZR = PTR*(ZR+ZR)*PTR 35 RZI = -PTR*(ZI+ZI)*PTR 36 T1R = RZR*FNUP 37 T1I = RZI*FNUP 38 P2R = -T1R 39 P2I = -T1I 40 P1R = CONER 41 P1I = CONEI 42 T1R = T1R + RZR 43 T1I = T1I + RZI 44 IF (ID.GT.0) ID = 0 45 AP2 = ZABS(CMPLX(P2R,P2I,kind=KIND(1.0D0))) 46 AP1 = ZABS(CMPLX(P1R,P1I,kind=KIND(1.0D0))) 47 C----------------------------------------------------------------------- 48 C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU 49 C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT 50 C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR 51 C PREMATURELY. 52 C----------------------------------------------------------------------- 53 ARG = (AP2+AP2)/(AP1*TOL) 54 TEST1 = DSQRT(ARG) 55 TEST = TEST1 56 RAP1 = 1.0D0/AP1 57 P1R = P1R*RAP1 58 P1I = P1I*RAP1 59 P2R = P2R*RAP1 60 P2I = P2I*RAP1 61 AP2 = AP2*RAP1 62 10 CONTINUE 63 K = K + 1 64 AP1 = AP2 65 PTR = P2R 66 PTI = P2I 67 P2R = P1R - (T1R*PTR-T1I*PTI) 68 P2I = P1I - (T1R*PTI+T1I*PTR) 69 P1R = PTR 70 P1I = PTI 71 T1R = T1R + RZR 72 T1I = T1I + RZI 73 AP2 = ZABS(CMPLX(P2R,P2I,kind=KIND(1.0D0))) 74 IF (AP1.LE.TEST) GO TO 10 75 IF (ITIME.EQ.2) GO TO 20 76 AK = ZABS(CMPLX(T1R,T1I,kind=KIND(1.0D0))*0.5D0) 77 FLAM = AK + DSQRT(AK*AK-1.0D0) 78 RHO = DMIN1(AP2/AP1,FLAM) 79 TEST = TEST1*DSQRT(RHO/(RHO*RHO-1.0D0)) 80 ITIME = 2 81 GO TO 10 82 20 CONTINUE 83 KK = K + 1 - ID 84 AK = DBLE(FLOAT(KK)) 85 T1R = AK 86 T1I = CZEROI 87 DFNU = FNU + DBLE(FLOAT(N-1)) 88 P1R = 1.0D0/AP2 89 P1I = CZEROI 90 P2R = CZEROR 91 P2I = CZEROI 92 DO 30 I=1,KK 93 PTR = P1R 94 PTI = P1I 95 RAP1 = DFNU + T1R 96 TTR = RZR*RAP1 97 TTI = RZI*RAP1 98 P1R = (PTR*TTR-PTI*TTI) + P2R 99 P1I = (PTR*TTI+PTI*TTR) + P2I 100 P2R = PTR 101 P2I = PTI 102 T1R = T1R - CONER 103 30 CONTINUE 104 IF (P1R.NE.CZEROR .OR. P1I.NE.CZEROI) GO TO 40 105 P1R = TOL 106 P1I = TOL 107 40 CONTINUE 108 CALL ZDIV(P2R, P2I, P1R, P1I, CYR(N), CYI(N)) 109 IF (N.EQ.1) RETURN 110 K = N - 1 111 AK = DBLE(FLOAT(K)) 112 T1R = AK 113 T1I = CZEROI 114 CDFNUR = FNU*RZR 115 CDFNUI = FNU*RZI 116 DO 60 I=2,N 117 PTR = CDFNUR + (T1R*RZR-T1I*RZI) + CYR(K+1) 118 PTI = CDFNUI + (T1R*RZI+T1I*RZR) + CYI(K+1) 119 AK = ZABS(CMPLX(PTR,PTI,kind=KIND(1.0D0))) 120 IF (AK.NE.CZEROR) GO TO 50 121 PTR = TOL 122 PTI = TOL 123 AK = TOL*RT2 124 50 CONTINUE 125 RAK = CONER/AK 126 CYR(K) = RAK*PTR*RAK 127 CYI(K) = -RAK*PTI*RAK 128 T1R = T1R - CONER 129 K = K - 1 130 60 CONTINUE 131 RETURN 132 END