dbesj0.f 3.1 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  1. *DECK DBESJ0
  2. DOUBLE PRECISION FUNCTION DBESJ0 (X)
  3. C***BEGIN PROLOGUE DBESJ0
  4. C***PURPOSE Compute the Bessel function of the first kind of order
  5. C zero.
  6. C***LIBRARY SLATEC (FNLIB)
  7. C***CATEGORY C10A1
  8. C***TYPE DOUBLE PRECISION (BESJ0-S, DBESJ0-D)
  9. C***KEYWORDS BESSEL FUNCTION, FIRST KIND, FNLIB, ORDER ZERO,
  10. C SPECIAL FUNCTIONS
  11. C***AUTHOR Fullerton, W., (LANL)
  12. C***DESCRIPTION
  13. C
  14. C DBESJ0(X) calculates the double precision Bessel function of
  15. C the first kind of order zero for double precision argument X.
  16. C
  17. C Series for BJ0 on the interval 0. to 1.60000E+01
  18. C with weighted error 4.39E-32
  19. C log weighted error 31.36
  20. C significant figures required 31.21
  21. C decimal places required 32.00
  22. C
  23. C***REFERENCES (NONE)
  24. C***ROUTINES CALLED D1MACH, D9B0MP, DCSEVL, INITDS
  25. C***REVISION HISTORY (YYMMDD)
  26. C 770701 DATE WRITTEN
  27. C 890531 Changed all specific intrinsics to generic. (WRB)
  28. C 890531 REVISION DATE from Version 3.2
  29. C 891214 Prologue converted to Version 4.0 format. (BAB)
  30. C***END PROLOGUE DBESJ0
  31. DOUBLE PRECISION X, BJ0CS(19), AMPL, THETA, XSML, Y, D1MACH,
  32. 1 DCSEVL
  33. LOGICAL FIRST
  34. SAVE BJ0CS, NTJ0, XSML, FIRST
  35. DATA BJ0CS( 1) / +.1002541619 6893913701 0731272640 74 D+0 /
  36. DATA BJ0CS( 2) / -.6652230077 6440513177 6787578311 24 D+0 /
  37. DATA BJ0CS( 3) / +.2489837034 9828131370 4604687266 80 D+0 /
  38. DATA BJ0CS( 4) / -.3325272317 0035769653 8843415038 54 D-1 /
  39. DATA BJ0CS( 5) / +.2311417930 4694015462 9049241177 29 D-2 /
  40. DATA BJ0CS( 6) / -.9911277419 9508092339 0485193365 49 D-4 /
  41. DATA BJ0CS( 7) / +.2891670864 3998808884 7339037470 78 D-5 /
  42. DATA BJ0CS( 8) / -.6121085866 3032635057 8184074815 16 D-7 /
  43. DATA BJ0CS( 9) / +.9838650793 8567841324 7687486364 15 D-9 /
  44. DATA BJ0CS( 10) / -.1242355159 7301765145 5158970068 36 D-10 /
  45. DATA BJ0CS( 11) / +.1265433630 2559045797 9158272103 63 D-12 /
  46. DATA BJ0CS( 12) / -.1061945649 5287244546 9148175129 59 D-14 /
  47. DATA BJ0CS( 13) / +.7470621075 8024567437 0989155840 00 D-17 /
  48. DATA BJ0CS( 14) / -.4469703227 4412780547 6270079999 99 D-19 /
  49. DATA BJ0CS( 15) / +.2302428158 4337436200 5230933333 33 D-21 /
  50. DATA BJ0CS( 16) / -.1031914479 4166698148 5226666666 66 D-23 /
  51. DATA BJ0CS( 17) / +.4060817827 4873322700 8000000000 00 D-26 /
  52. DATA BJ0CS( 18) / -.1414383600 5240913919 9999999999 99 D-28 /
  53. DATA BJ0CS( 19) / +.4391090549 6698880000 0000000000 00 D-31 /
  54. DATA FIRST /.TRUE./
  55. C***FIRST EXECUTABLE STATEMENT DBESJ0
  56. IF (FIRST) THEN
  57. NTJ0 = INITDS (BJ0CS, 19, 0.1*REAL(D1MACH(3)))
  58. XSML = SQRT(8.0D0*D1MACH(3))
  59. ENDIF
  60. FIRST = .FALSE.
  61. C
  62. Y = ABS(X)
  63. IF (Y.GT.4.0D0) GO TO 20
  64. C
  65. DBESJ0 = 1.0D0
  66. IF (Y.GT.XSML) DBESJ0 = DCSEVL (.125D0*Y*Y-1.D0, BJ0CS, NTJ0)
  67. RETURN
  68. C
  69. 20 CALL D9B0MP (Y, AMPL, THETA)
  70. DBESJ0 = AMPL * COS(THETA)
  71. C
  72. RETURN
  73. END