1 SUBROUTINE ZRATI(ZR, ZI, FNU, N, CYR, CYI, TOL) 2C***BEGIN PROLOGUE ZRATI 3C***REFER TO ZBESI,ZBESK,ZBESH 4C 5C ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD 6C RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD 7C RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B, 8C MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973, 9C BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER, 10C BY D. J. SOOKNE. 11C 12C***ROUTINES CALLED AZABS,ZDIV 13C***END PROLOGUE ZRATI 14C 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, AZABS 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 = AZABS(ZR,ZI) 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 = AZABS(P2R,P2I) 46 AP1 = AZABS(P1R,P1I) 47C----------------------------------------------------------------------- 48C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU 49C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT 50C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR 51C PREMATURELY. 52C----------------------------------------------------------------------- 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 = AZABS(P2R,P2I) 74 IF (AP1.LE.TEST) GO TO 10 75 IF (ITIME.EQ.2) GO TO 20 76 AK = AZABS(T1R,T1I)*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 = AZABS(PTR,PTI) 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 133