src/numerical/slatec/fortran/dbesy0.f

   1 *DECK DBESY0
   2       DOUBLE PRECISION FUNCTION DBESY0 (X)
   3 C***BEGIN PROLOGUE  DBESY0
   4 C***PURPOSE  Compute the Bessel function of the second kind of order
   5 C            zero.
   6 C***LIBRARY   SLATEC (FNLIB)
   7 C***CATEGORY  C10A1
   8 C***TYPE      DOUBLE PRECISION (BESY0-S, DBESY0-D)
   9 C***KEYWORDS  BESSEL FUNCTION, FNLIB, ORDER ZERO, SECOND KIND,
  10 C             SPECIAL FUNCTIONS
  11 C***AUTHOR  Fullerton, W., (LANL)
  12 C***DESCRIPTION
  13 C
  14 C DBESY0(X) calculates the double precision Bessel function of the
  15 C second kind of order zero for double precision argument X.
  16 C
  17 C Series for BY0        on the interval  0.          to  1.60000E+01
  18 C                                        with weighted error   8.14E-32
  19 C                                         log weighted error  31.09
  20 C                               significant figures required  30.31
  21 C                                    decimal places required  31.73
  22 C
  23 C***REFERENCES  (NONE)
  24 C***ROUTINES CALLED  D1MACH, D9B0MP, DBESJ0, DCSEVL, INITDS, XERMSG
  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   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
  31 C***END PROLOGUE  DBESY0
  32       DOUBLE PRECISION X, BY0CS(19), AMPL, THETA, TWODPI, XSML,
  33      1  Y, D1MACH, DCSEVL, DBESJ0
  34       LOGICAL FIRST
  35       SAVE BY0CS, TWODPI, NTY0, XSML, FIRST
  36       DATA BY0CS(  1) / -.1127783939 2865573217 9398054602 8 D-1      /
  37       DATA BY0CS(  2) / -.1283452375 6042034604 8088453183 8 D+0      /
  38       DATA BY0CS(  3) / -.1043788479 9794249365 8176227661 8 D+0      /
  39       DATA BY0CS(  4) / +.2366274918 3969695409 2415926461 3 D-1      /
  40       DATA BY0CS(  5) / -.2090391647 7004862391 9622395034 2 D-2      /
  41       DATA BY0CS(  6) / +.1039754539 3905725209 9924657638 1 D-3      /
  42       DATA BY0CS(  7) / -.3369747162 4239720967 1877534503 7 D-5      /
  43       DATA BY0CS(  8) / +.7729384267 6706671585 2136721637 1 D-7      /
  44       DATA BY0CS(  9) / -.1324976772 6642595914 4347606896 4 D-8      /
  45       DATA BY0CS( 10) / +.1764823261 5404527921 0038936315 8 D-10     /
  46       DATA BY0CS( 11) / -.1881055071 5801962006 0282301206 9 D-12     /
  47       DATA BY0CS( 12) / +.1641865485 3661495027 9223718574 9 D-14     /
  48       DATA BY0CS( 13) / -.1195659438 6046060857 4599100672 0 D-16     /
  49       DATA BY0CS( 14) / +.7377296297 4401858424 9411242666 6 D-19     /
  50       DATA BY0CS( 15) / -.3906843476 7104373307 4090666666 6 D-21     /
  51       DATA BY0CS( 16) / +.1795503664 4361579498 2912000000 0 D-23     /
  52       DATA BY0CS( 17) / -.7229627125 4480104789 3333333333 3 D-26     /
  53       DATA BY0CS( 18) / +.2571727931 6351685973 3333333333 3 D-28     /
  54       DATA BY0CS( 19) / -.8141268814 1636949333 3333333333 3 D-31     /
  55       DATA TWODPI / 0.6366197723 6758134307 5535053490 057 D0 /
  56       DATA FIRST /.TRUE./
  57 C***FIRST EXECUTABLE STATEMENT  DBESY0
  58       IF (FIRST) THEN
  59          NTY0 = INITDS (BY0CS, 19, 0.1*REAL(D1MACH(3)))
  60          XSML = SQRT(4.0D0*D1MACH(3))
  61       ENDIF
  62       FIRST = .FALSE.
  63 C
  64       IF (X .LE. 0.D0) CALL XERMSG ('SLATEC', 'DBESY0',
  65      +   'X IS ZERO OR NEGATIVE', 1, 2)
  66       IF (X.GT.4.0D0) GO TO 20
  67 C
  68       Y = 0.D0
  69       IF (X.GT.XSML) Y = X*X
  70       DBESY0 = TWODPI*LOG(0.5D0*X)*DBESJ0(X) + .375D0 + DCSEVL (
  71      1  .125D0*Y-1.D0, BY0CS, NTY0)
  72       RETURN
  73 C
  74  20   CALL D9B0MP (X, AMPL, THETA)
  75       DBESY0 = AMPL * SIN(THETA)
  76       RETURN
  77 C
  78       END