zasyi.f 4.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165
  1. SUBROUTINE ZASYI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM,
  2. * ALIM)
  3. C***BEGIN PROLOGUE ZASYI
  4. C***REFER TO ZBESI,ZBESK
  5. C
  6. C ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
  7. C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
  8. C REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
  9. C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
  10. C
  11. C***ROUTINES CALLED D1MACH,ZABS,ZDIV,ZEXP,ZMLT,ZSQRT
  12. C***END PROLOGUE ZASYI
  13. C COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z
  14. DOUBLE PRECISION AA, AEZ, AK, AK1I, AK1R, ALIM, ARG, ARM, ATOL,
  15. * AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI,
  16. * CZR, DFNU, DKI, DKR, DNU2, ELIM, EZI, EZR, FDN, FNU, PI, P1I,
  17. * P1R, RAZ, RL, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I,
  18. * S2R, TOL, TZI, TZR, YI, YR, ZEROI, ZEROR, ZI, ZR, D1MACH, ZABS
  19. INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ
  20. DIMENSION YR(N), YI(N)
  21. DATA PI, RTPI /3.14159265358979324D0 , 0.159154943091895336D0 /
  22. DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
  23. C
  24. NZ = 0
  25. AZ = ZABS(COMPLEX(ZR,ZI))
  26. ARM = 1.0D+3*D1MACH(1)
  27. RTR1 = DSQRT(ARM)
  28. IL = MIN0(2,N)
  29. DFNU = FNU + DBLE(FLOAT(N-IL))
  30. C-----------------------------------------------------------------------
  31. C OVERFLOW TEST
  32. C-----------------------------------------------------------------------
  33. RAZ = 1.0D0/AZ
  34. STR = ZR*RAZ
  35. STI = -ZI*RAZ
  36. AK1R = RTPI*STR*RAZ
  37. AK1I = RTPI*STI*RAZ
  38. CALL ZSQRT(AK1R, AK1I, AK1R, AK1I)
  39. CZR = ZR
  40. CZI = ZI
  41. IF (KODE.NE.2) GO TO 10
  42. CZR = ZEROR
  43. CZI = ZI
  44. 10 CONTINUE
  45. IF (DABS(CZR).GT.ELIM) GO TO 100
  46. DNU2 = DFNU + DFNU
  47. KODED = 1
  48. IF ((DABS(CZR).GT.ALIM) .AND. (N.GT.2)) GO TO 20
  49. KODED = 0
  50. CALL ZEXP(CZR, CZI, STR, STI)
  51. CALL ZMLT(AK1R, AK1I, STR, STI, AK1R, AK1I)
  52. 20 CONTINUE
  53. FDN = 0.0D0
  54. IF (DNU2.GT.RTR1) FDN = DNU2*DNU2
  55. EZR = ZR*8.0D0
  56. EZI = ZI*8.0D0
  57. C-----------------------------------------------------------------------
  58. C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
  59. C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
  60. C EXPANSION FOR THE IMAGINARY PART.
  61. C-----------------------------------------------------------------------
  62. AEZ = 8.0D0*AZ
  63. S = TOL/AEZ
  64. JL = INT(SNGL(RL+RL)) + 2
  65. P1R = ZEROR
  66. P1I = ZEROI
  67. IF (ZI.EQ.0.0D0) GO TO 30
  68. C-----------------------------------------------------------------------
  69. C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
  70. C SIGNIFICANCE WHEN FNU OR N IS LARGE
  71. C-----------------------------------------------------------------------
  72. INU = INT(SNGL(FNU))
  73. ARG = (FNU-DBLE(FLOAT(INU)))*PI
  74. INU = INU + N - IL
  75. AK = -DSIN(ARG)
  76. BK = DCOS(ARG)
  77. IF (ZI.LT.0.0D0) BK = -BK
  78. P1R = AK
  79. P1I = BK
  80. IF (MOD(INU,2).EQ.0) GO TO 30
  81. P1R = -P1R
  82. P1I = -P1I
  83. 30 CONTINUE
  84. DO 70 K=1,IL
  85. SQK = FDN - 1.0D0
  86. ATOL = S*DABS(SQK)
  87. SGN = 1.0D0
  88. CS1R = CONER
  89. CS1I = CONEI
  90. CS2R = CONER
  91. CS2I = CONEI
  92. CKR = CONER
  93. CKI = CONEI
  94. AK = 0.0D0
  95. AA = 1.0D0
  96. BB = AEZ
  97. DKR = EZR
  98. DKI = EZI
  99. DO 40 J=1,JL
  100. CALL ZDIV(CKR, CKI, DKR, DKI, STR, STI)
  101. CKR = STR*SQK
  102. CKI = STI*SQK
  103. CS2R = CS2R + CKR
  104. CS2I = CS2I + CKI
  105. SGN = -SGN
  106. CS1R = CS1R + CKR*SGN
  107. CS1I = CS1I + CKI*SGN
  108. DKR = DKR + EZR
  109. DKI = DKI + EZI
  110. AA = AA*DABS(SQK)/BB
  111. BB = BB + AEZ
  112. AK = AK + 8.0D0
  113. SQK = SQK - AK
  114. IF (AA.LE.ATOL) GO TO 50
  115. 40 CONTINUE
  116. GO TO 110
  117. 50 CONTINUE
  118. S2R = CS1R
  119. S2I = CS1I
  120. IF (ZR+ZR.GE.ELIM) GO TO 60
  121. TZR = ZR + ZR
  122. TZI = ZI + ZI
  123. CALL ZEXP(-TZR, -TZI, STR, STI)
  124. CALL ZMLT(STR, STI, P1R, P1I, STR, STI)
  125. CALL ZMLT(STR, STI, CS2R, CS2I, STR, STI)
  126. S2R = S2R + STR
  127. S2I = S2I + STI
  128. 60 CONTINUE
  129. FDN = FDN + 8.0D0*DFNU + 4.0D0
  130. P1R = -P1R
  131. P1I = -P1I
  132. M = N - IL + K
  133. YR(M) = S2R*AK1R - S2I*AK1I
  134. YI(M) = S2R*AK1I + S2I*AK1R
  135. 70 CONTINUE
  136. IF (N.LE.2) RETURN
  137. NN = N
  138. K = NN - 2
  139. AK = DBLE(FLOAT(K))
  140. STR = ZR*RAZ
  141. STI = -ZI*RAZ
  142. RZR = (STR+STR)*RAZ
  143. RZI = (STI+STI)*RAZ
  144. IB = 3
  145. DO 80 I=IB,NN
  146. YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2)
  147. YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2)
  148. AK = AK - 1.0D0
  149. K = K - 1
  150. 80 CONTINUE
  151. IF (KODED.EQ.0) RETURN
  152. CALL ZEXP(CZR, CZI, CKR, CKI)
  153. DO 90 I=1,NN
  154. STR = YR(I)*CKR - YI(I)*CKI
  155. YI(I) = YR(I)*CKI + YI(I)*CKR
  156. YR(I) = STR
  157. 90 CONTINUE
  158. RETURN
  159. 100 CONTINUE
  160. NZ = -1
  161. RETURN
  162. 110 CONTINUE
  163. NZ=-2
  164. RETURN
  165. END