zrati.f 3.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132
  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(COMPLEX(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 = ZABS(COMPLEX(P2R,P2I))
  46. AP1 = ZABS(COMPLEX(P1R,P1I))
  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(COMPLEX(P2R,P2I))
  74. IF (AP1.LE.TEST) GO TO 10
  75. IF (ITIME.EQ.2) GO TO 20
  76. AK = ZABS(COMPLEX(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 = ZABS(COMPLEX(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