zbesy.f 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244
  1. SUBROUTINE ZBESY(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, CWRKR, CWRKI,
  2. * IERR)
  3. C***BEGIN PROLOGUE ZBESY
  4. C***DATE WRITTEN 830501 (YYMMDD)
  5. C***REVISION DATE 890801 (YYMMDD)
  6. C***CATEGORY NO. B5K
  7. C***KEYWORDS Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT,
  8. C BESSEL FUNCTION OF SECOND KIND
  9. C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
  10. C***PURPOSE TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT
  11. C***DESCRIPTION
  12. C
  13. C ***A DOUBLE PRECISION ROUTINE***
  14. C
  15. C ON KODE=1, CBESY COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
  16. C BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE
  17. C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
  18. C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESY RETURNS THE SCALED
  19. C FUNCTIONS
  20. C
  21. C CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z)
  22. C
  23. C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
  24. C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
  25. C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
  26. C (REF. 1).
  27. C
  28. C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION
  29. C ZR,ZI - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0),
  30. C -PI.LT.ARG(Z).LE.PI
  31. C FNU - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0D0
  32. C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
  33. C KODE= 1 RETURNS
  34. C CY(I)=Y(FNU+I-1,Z), I=1,...,N
  35. C = 2 RETURNS
  36. C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N
  37. C WHERE Y=AIMAG(Z)
  38. C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
  39. C CWRKR, - DOUBLE PRECISION WORK VECTORS OF DIMENSION AT
  40. C CWRKI AT LEAST N
  41. C
  42. C OUTPUT CYR,CYI ARE DOUBLE PRECISION
  43. C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
  44. C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
  45. C CY(I)=Y(FNU+I-1,Z) OR
  46. C CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N
  47. C DEPENDING ON KODE.
  48. C NZ - NZ=0 , A NORMAL RETURN
  49. C NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO
  50. C UNDERFLOW (GENERALLY ON KODE=2)
  51. C IERR - ERROR FLAG
  52. C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
  53. C IERR=1, INPUT ERROR - NO COMPUTATION
  54. C IERR=2, OVERFLOW - NO COMPUTATION, FNU IS
  55. C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
  56. C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
  57. C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
  58. C REDUCTION PRODUCE LESS THAN HALF OF MACHINE
  59. C ACCURACY
  60. C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
  61. C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
  62. C CANCE BY ARGUMENT REDUCTION
  63. C IERR=5, ERROR - NO COMPUTATION,
  64. C ALGORITHM TERMINATION CONDITION NOT MET
  65. C
  66. C***LONG DESCRIPTION
  67. C
  68. C THE COMPUTATION IS CARRIED OUT BY THE FORMULA
  69. C
  70. C Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I
  71. C
  72. C WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z)
  73. C AND H(2,FNU,Z) ARE CALCULATED IN CBESH.
  74. C
  75. C FOR NEGATIVE ORDERS,THE FORMULA
  76. C
  77. C Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
  78. C
  79. C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD
  80. C INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE
  81. C POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)*
  82. C SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS
  83. C NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A
  84. C LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM
  85. C CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS,
  86. C WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF
  87. C ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z).
  88. C
  89. C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
  90. C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
  91. C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
  92. C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
  93. C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
  94. C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
  95. C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
  96. C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
  97. C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
  98. C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
  99. C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
  100. C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
  101. C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
  102. C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
  103. C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
  104. C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
  105. C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
  106. C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
  107. C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
  108. C
  109. C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
  110. C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
  111. C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
  112. C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
  113. C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
  114. C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
  115. C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
  116. C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
  117. C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
  118. C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
  119. C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
  120. C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
  121. C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
  122. C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
  123. C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
  124. C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
  125. C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
  126. C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
  127. C OR -PI/2+P.
  128. C
  129. C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
  130. C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
  131. C COMMERCE, 1955.
  132. C
  133. C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
  134. C BY D. E. AMOS, SAND83-0083, MAY, 1983.
  135. C
  136. C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
  137. C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
  138. C
  139. C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
  140. C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
  141. C 1018, MAY, 1985
  142. C
  143. C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
  144. C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
  145. C MATH. SOFTWARE, 1986
  146. C
  147. C***ROUTINES CALLED ZBESH,I1MACH,D1MACH
  148. C***END PROLOGUE ZBESY
  149. C
  150. C COMPLEX CWRK,CY,C1,C2,EX,HCI,Z,ZU,ZV
  151. DOUBLE PRECISION CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2R,
  152. * ELIM, EXI, EXR, EY, FNU, HCII, STI, STR, TAY, ZI, ZR, DEXP,
  153. * D1MACH, ASCLE, RTOL, ATOL, AA, BB, TOL
  154. INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH
  155. DIMENSION CYR(N), CYI(N), CWRKR(N), CWRKI(N)
  156. C***FIRST EXECUTABLE STATEMENT ZBESY
  157. IERR = 0
  158. NZ=0
  159. IF (ZR.EQ.0.0D0 .AND. ZI.EQ.0.0D0) IERR=1
  160. IF (FNU.LT.0.0D0) IERR=1
  161. IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
  162. IF (N.LT.1) IERR=1
  163. IF (IERR.NE.0) RETURN
  164. HCII = 0.5D0
  165. CALL ZBESH(ZR, ZI, FNU, KODE, 1, N, CYR, CYI, NZ1, IERR)
  166. IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
  167. CALL ZBESH(ZR, ZI, FNU, KODE, 2, N, CWRKR, CWRKI, NZ2, IERR)
  168. IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
  169. NZ = MIN0(NZ1,NZ2)
  170. IF (KODE.EQ.2) GO TO 60
  171. DO 50 I=1,N
  172. STR = CWRKR(I) - CYR(I)
  173. STI = CWRKI(I) - CYI(I)
  174. CYR(I) = -STI*HCII
  175. CYI(I) = STR*HCII
  176. 50 CONTINUE
  177. RETURN
  178. 60 CONTINUE
  179. TOL = DMAX1(D1MACH(4),1.0D-18)
  180. K1 = I1MACH(15)
  181. K2 = I1MACH(16)
  182. K = MIN0(IABS(K1),IABS(K2))
  183. R1M5 = D1MACH(5)
  184. C-----------------------------------------------------------------------
  185. C ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
  186. C-----------------------------------------------------------------------
  187. ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
  188. EXR = DCOS(ZR)
  189. EXI = DSIN(ZR)
  190. EY = 0.0D0
  191. TAY = DABS(ZI+ZI)
  192. IF (TAY.LT.ELIM) EY = DEXP(-TAY)
  193. IF (ZI.LT.0.0D0) GO TO 90
  194. C1R = EXR*EY
  195. C1I = EXI*EY
  196. C2R = EXR
  197. C2I = -EXI
  198. 70 CONTINUE
  199. NZ = 0
  200. RTOL = 1.0D0/TOL
  201. ASCLE = D1MACH(1)*RTOL*1.0D+3
  202. DO 80 I=1,N
  203. C STR = C1R*CYR(I) - C1I*CYI(I)
  204. C STI = C1R*CYI(I) + C1I*CYR(I)
  205. C STR = -STR + C2R*CWRKR(I) - C2I*CWRKI(I)
  206. C STI = -STI + C2R*CWRKI(I) + C2I*CWRKR(I)
  207. C CYR(I) = -STI*HCII
  208. C CYI(I) = STR*HCII
  209. AA = CWRKR(I)
  210. BB = CWRKI(I)
  211. ATOL = 1.0D0
  212. IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 75
  213. AA = AA*RTOL
  214. BB = BB*RTOL
  215. ATOL = TOL
  216. 75 CONTINUE
  217. STR = (AA*C2R - BB*C2I)*ATOL
  218. STI = (AA*C2I + BB*C2R)*ATOL
  219. AA = CYR(I)
  220. BB = CYI(I)
  221. ATOL = 1.0D0
  222. IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 85
  223. AA = AA*RTOL
  224. BB = BB*RTOL
  225. ATOL = TOL
  226. 85 CONTINUE
  227. STR = STR - (AA*C1R - BB*C1I)*ATOL
  228. STI = STI - (AA*C1I + BB*C1R)*ATOL
  229. CYR(I) = -STI*HCII
  230. CYI(I) = STR*HCII
  231. IF (STR.EQ.0.0D0 .AND. STI.EQ.0.0D0 .AND. EY.EQ.0.0D0) NZ = NZ
  232. * + 1
  233. 80 CONTINUE
  234. RETURN
  235. 90 CONTINUE
  236. C1R = EXR
  237. C1I = EXI
  238. C2R = EXR*EY
  239. C2I = -EXI*EY
  240. GO TO 70
  241. 170 CONTINUE
  242. NZ = 0
  243. RETURN
  244. END