1*DECK ZRATI 2 SUBROUTINE ZRATI (ZR, ZI, FNU, N, CYR, CYI, TOL) 3C***BEGIN PROLOGUE ZRATI 4C***SUBSIDIARY 5C***PURPOSE Subsidiary to ZBESH, ZBESI and ZBESK 6C***LIBRARY SLATEC 7C***TYPE ALL (CRATI-A, ZRATI-A) 8C***AUTHOR Amos, D. E., (SNL) 9C***DESCRIPTION 10C 11C ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD 12C RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD 13C RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B, 14C MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973, 15C BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER, 16C BY D. J. SOOKNE. 17C 18C***SEE ALSO ZBESH, ZBESI, ZBESK 19C***ROUTINES CALLED ZABS, ZDIV 20C***REVISION HISTORY (YYMMDD) 21C 830501 DATE WRITTEN 22C 910415 Prologue converted to Version 4.0 format. (BAB) 23C***END PROLOGUE ZRATI 24 DOUBLE PRECISION AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR, 25 * CONEI, CONER, CYI, CYR, CZEROI, CZEROR, DFNU, FDNU, FLAM, FNU, 26 * FNUP, PTI, PTR, P1I, P1R, P2I, P2R, RAK, RAP1, RHO, RT2, RZI, 27 * RZR, TEST, TEST1, TOL, TTI, TTR, T1I, T1R, ZI, ZR, ZABS 28 INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N 29 DIMENSION CYR(N), CYI(N) 30 EXTERNAL ZABS 31 DATA CZEROR,CZEROI,CONER,CONEI,RT2/ 32 1 0.0D0, 0.0D0, 1.0D0, 0.0D0, 1.41421356237309505D0 / 33C***FIRST EXECUTABLE STATEMENT ZRATI 34 AZ = ZABS(ZR,ZI) 35 INU = FNU 36 IDNU = INU + N - 1 37 MAGZ = AZ 38 AMAGZ = MAGZ+1 39 FDNU = IDNU 40 FNUP = MAX(AMAGZ,FDNU) 41 ID = IDNU - MAGZ - 1 42 ITIME = 1 43 K = 1 44 PTR = 1.0D0/AZ 45 RZR = PTR*(ZR+ZR)*PTR 46 RZI = -PTR*(ZI+ZI)*PTR 47 T1R = RZR*FNUP 48 T1I = RZI*FNUP 49 P2R = -T1R 50 P2I = -T1I 51 P1R = CONER 52 P1I = CONEI 53 T1R = T1R + RZR 54 T1I = T1I + RZI 55 IF (ID.GT.0) ID = 0 56 AP2 = ZABS(P2R,P2I) 57 AP1 = ZABS(P1R,P1I) 58C----------------------------------------------------------------------- 59C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU 60C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT 61C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR 62C PREMATURELY. 63C----------------------------------------------------------------------- 64 ARG = (AP2+AP2)/(AP1*TOL) 65 TEST1 = SQRT(ARG) 66 TEST = TEST1 67 RAP1 = 1.0D0/AP1 68 P1R = P1R*RAP1 69 P1I = P1I*RAP1 70 P2R = P2R*RAP1 71 P2I = P2I*RAP1 72 AP2 = AP2*RAP1 73 10 CONTINUE 74 K = K + 1 75 AP1 = AP2 76 PTR = P2R 77 PTI = P2I 78 P2R = P1R - (T1R*PTR-T1I*PTI) 79 P2I = P1I - (T1R*PTI+T1I*PTR) 80 P1R = PTR 81 P1I = PTI 82 T1R = T1R + RZR 83 T1I = T1I + RZI 84 AP2 = ZABS(P2R,P2I) 85 IF (AP1.LE.TEST) GO TO 10 86 IF (ITIME.EQ.2) GO TO 20 87 AK = ZABS(T1R,T1I)*0.5D0 88 FLAM = AK + SQRT(AK*AK-1.0D0) 89 RHO = MIN(AP2/AP1,FLAM) 90 TEST = TEST1*SQRT(RHO/(RHO*RHO-1.0D0)) 91 ITIME = 2 92 GO TO 10 93 20 CONTINUE 94 KK = K + 1 - ID 95 AK = KK 96 T1R = AK 97 T1I = CZEROI 98 DFNU = FNU + (N-1) 99 P1R = 1.0D0/AP2 100 P1I = CZEROI 101 P2R = CZEROR 102 P2I = CZEROI 103 DO 30 I=1,KK 104 PTR = P1R 105 PTI = P1I 106 RAP1 = DFNU + T1R 107 TTR = RZR*RAP1 108 TTI = RZI*RAP1 109 P1R = (PTR*TTR-PTI*TTI) + P2R 110 P1I = (PTR*TTI+PTI*TTR) + P2I 111 P2R = PTR 112 P2I = PTI 113 T1R = T1R - CONER 114 30 CONTINUE 115 IF (P1R.NE.CZEROR .OR. P1I.NE.CZEROI) GO TO 40 116 P1R = TOL 117 P1I = TOL 118 40 CONTINUE 119 CALL ZDIV(P2R, P2I, P1R, P1I, CYR(N), CYI(N)) 120 IF (N.EQ.1) RETURN 121 K = N - 1 122 AK = K 123 T1R = AK 124 T1I = CZEROI 125 CDFNUR = FNU*RZR 126 CDFNUI = FNU*RZI 127 DO 60 I=2,N 128 PTR = CDFNUR + (T1R*RZR-T1I*RZI) + CYR(K+1) 129 PTI = CDFNUI + (T1R*RZI+T1I*RZR) + CYI(K+1) 130 AK = ZABS(PTR,PTI) 131 IF (AK.NE.CZEROR) GO TO 50 132 PTR = TOL 133 PTI = TOL 134 AK = TOL*RT2 135 50 CONTINUE 136 RAK = CONER/AK 137 CYR(K) = RAK*PTR*RAK 138 CYI(K) = -RAK*PTI*RAK 139 T1R = T1R - CONER 140 K = K - 1 141 60 CONTINUE 142 RETURN 143 END 144