zmlri.f 6.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204
  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(COMPLEX(ZR,ZI))
  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(COMPLEX(P2R,P2I))
  56. IF (AP.GT.TST*AK*AK) GO TO 20
  57. AK = AK + 1.0D0
  58. 10 CONTINUE
  59. GO TO 110
  60. 20 CONTINUE
  61. I = I + 1
  62. K = 0
  63. IF (INU.LT.IAZ) GO TO 40
  64. C-----------------------------------------------------------------------
  65. C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
  66. C-----------------------------------------------------------------------
  67. P1R = ZEROR
  68. P1I = ZEROI
  69. P2R = CONER
  70. P2I = CONEI
  71. AT = DBLE(FLOAT(INU)) + 1.0D0
  72. STR = ZR*RAZ
  73. STI = -ZI*RAZ
  74. CKR = STR*AT*RAZ
  75. CKI = STI*AT*RAZ
  76. ACK = AT*RAZ
  77. TST = DSQRT(ACK/TOL)
  78. ITIME = 1
  79. DO 30 K=1,80
  80. PTR = P2R
  81. PTI = P2I
  82. P2R = P1R - (CKR*PTR-CKI*PTI)
  83. P2I = P1I - (CKR*PTI+CKI*PTR)
  84. P1R = PTR
  85. P1I = PTI
  86. CKR = CKR + RZR
  87. CKI = CKI + RZI
  88. AP = ZABS(COMPLEX(P2R,P2I))
  89. IF (AP.LT.TST) GO TO 30
  90. IF (ITIME.EQ.2) GO TO 40
  91. ACK = ZABS(COMPLEX(CKR,CKI))
  92. FLAM = ACK + DSQRT(ACK*ACK-1.0D0)
  93. FKAP = AP/ZABS(COMPLEX(P1R,P1I))
  94. RHO = DMIN1(FLAM,FKAP)
  95. TST = TST*DSQRT(RHO/(RHO*RHO-1.0D0))
  96. ITIME = 2
  97. 30 CONTINUE
  98. GO TO 110
  99. 40 CONTINUE
  100. C-----------------------------------------------------------------------
  101. C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
  102. C-----------------------------------------------------------------------
  103. K = K + 1
  104. KK = MAX0(I+IAZ,K+INU)
  105. FKK = DBLE(FLOAT(KK))
  106. P1R = ZEROR
  107. P1I = ZEROI
  108. C-----------------------------------------------------------------------
  109. C SCALE P2 AND SUM BY SCLE
  110. C-----------------------------------------------------------------------
  111. P2R = SCLE
  112. P2I = ZEROI
  113. FNF = FNU - DBLE(FLOAT(IFNU))
  114. TFNF = FNF + FNF
  115. BK = DGAMLN(FKK+TFNF+1.0D0,IDUM) - DGAMLN(FKK+1.0D0,IDUM) -
  116. * DGAMLN(TFNF+1.0D0,IDUM)
  117. BK = DEXP(BK)
  118. SUMR = ZEROR
  119. SUMI = ZEROI
  120. KM = KK - INU
  121. DO 50 I=1,KM
  122. PTR = P2R
  123. PTI = P2I
  124. P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
  125. P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
  126. P1R = PTR
  127. P1I = PTI
  128. AK = 1.0D0 - TFNF/(FKK+TFNF)
  129. ACK = BK*AK
  130. SUMR = SUMR + (ACK+BK)*P1R
  131. SUMI = SUMI + (ACK+BK)*P1I
  132. BK = ACK
  133. FKK = FKK - 1.0D0
  134. 50 CONTINUE
  135. YR(N) = P2R
  136. YI(N) = P2I
  137. IF (N.EQ.1) GO TO 70
  138. DO 60 I=2,N
  139. PTR = P2R
  140. PTI = P2I
  141. P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
  142. P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI)
  143. P1R = PTR
  144. P1I = PTI
  145. AK = 1.0D0 - TFNF/(FKK+TFNF)
  146. ACK = BK*AK
  147. SUMR = SUMR + (ACK+BK)*P1R
  148. SUMI = SUMI + (ACK+BK)*P1I
  149. BK = ACK
  150. FKK = FKK - 1.0D0
  151. M = N - I + 1
  152. YR(M) = P2R
  153. YI(M) = P2I
  154. 60 CONTINUE
  155. 70 CONTINUE
  156. IF (IFNU.LE.0) GO TO 90
  157. DO 80 I=1,IFNU
  158. PTR = P2R
  159. PTI = P2I
  160. P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI)
  161. P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR)
  162. P1R = PTR
  163. P1I = PTI
  164. AK = 1.0D0 - TFNF/(FKK+TFNF)
  165. ACK = BK*AK
  166. SUMR = SUMR + (ACK+BK)*P1R
  167. SUMI = SUMI + (ACK+BK)*P1I
  168. BK = ACK
  169. FKK = FKK - 1.0D0
  170. 80 CONTINUE
  171. 90 CONTINUE
  172. PTR = ZR
  173. PTI = ZI
  174. IF (KODE.EQ.2) PTR = ZEROR
  175. CALL ZLOG(RZR, RZI, STR, STI, IDUM)
  176. P1R = -FNF*STR + PTR
  177. P1I = -FNF*STI + PTI
  178. AP = DGAMLN(1.0D0+FNF,IDUM)
  179. PTR = P1R - AP
  180. PTI = P1I
  181. C-----------------------------------------------------------------------
  182. C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
  183. C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
  184. C-----------------------------------------------------------------------
  185. P2R = P2R + SUMR
  186. P2I = P2I + SUMI
  187. AP = ZABS(COMPLEX(P2R,P2I))
  188. P1R = 1.0D0/AP
  189. CALL ZEXP(PTR, PTI, STR, STI)
  190. CKR = STR*P1R
  191. CKI = STI*P1R
  192. PTR = P2R*P1R
  193. PTI = -P2I*P1R
  194. CALL ZMLT(CKR, CKI, PTR, PTI, CNORMR, CNORMI)
  195. DO 100 I=1,N
  196. STR = YR(I)*CNORMR - YI(I)*CNORMI
  197. YI(I) = YR(I)*CNORMI + YI(I)*CNORMR
  198. YR(I) = STR
  199. 100 CONTINUE
  200. RETURN
  201. 110 CONTINUE
  202. NZ=-2
  203. RETURN
  204. END