      SUBROUTINE VY0 (M, X, F, WORK, IWORK, INFO)
C***BEGIN PROLOGUE  VY0
C***PURPOSE  Computes the Bessel function of the second kind
C            of order zero (Y0) for a vector of real arguments
C***LIBRARY   VFNLIB
C***CATEGORY  C10A1
C***TYPE      SINGLE PRECISION (VY0-S, DVY0-D)
C***KEYWORDS  BESSEL FUNCTION,SECOND KIND, ORDER ZERO, VECTORIZED
C***AUTHOR  SAUNDERS, B. V., (NIST)
C             COMPUTING AND APPLIED MATHEMATICS LABORATORY
C             NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C             GAITHERSBURG, MD 20899
C           BOISVERT, R. F., (NIST)
C             COMPUTING AND APPLIED MATHEMATICS LABORATORY
C             NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C             GAITHERSBURG, MD 20899
C***DESCRIPTION
C
C   VY0 computes the  Bessel function of the second kind of order
C   zero (Y0) for real arguments using uniform approximation by
C   Chebyshev polynomials.
C
C
C   P A R A M E T E R S
C
C   M      (Input) Integer (M .GT. 0)
C          The number of arguments at which the function is to be
C          evaluated.
C
C   X      (Input) Real array of length M
C          The arguments at which the function is to be evaluated are
C          stored in X(1) to X(M) in any order.
C
C   F      (Output) Real array of length M
C          F(i) contains the value of the function at X(i), i=1,..,M.
C
C   WORK   (Work) Real vector of length 7*M
C
C   IWORK  (Work) Integer vector of length M
C
C   INFO   (Output) Integer
C          Indicates status of computed result. The following table
C          lists possible values and their meanings.  If OK=Yes then
C          all F(i) have been set, otherwise none have been set.
C
C          INFO  OK            Description
C          ------------------------------------------------------------
C            0   Yes  Successfull execution.
C            1   No   Error: M .LE. 0
C            2   No   Error: X(i) zero or negative for some i
C            3   No   Error: Some X(i) so big that no precision possible
C                     in computing Y0. The index of the first offending
C                     argument is returned in IWORK(1).
C
C
C *********************************************************************
C   This routine is a modification of the function BESY0 developed by
C   W. Fullerton of LANL.
C *********************************************************************
C
C***REFERENCES  Ronald F. Boisvert and Bonita V. Saunders, Portable
C               Vectorized Software for Bessel Function Evaluation,
C               ACM Transactions on Mathematical Software 18 (1992),
C               pp 456-469.
C***ROUTINES CALLED  WY0
C***REVISION HISTORY  (YYMMDD)
C   910226  DATE WRITTEN
C   930203  Minor modifications to prologue.
C***END PROLOGUE  VY0
C
C     ----------
C     PARAMETERS
C     ----------
C
      INTEGER INFO, IWORK, M
      REAL    F, X, WORK
C
      DIMENSION X(M), F(M), WORK(7*M), IWORK(M)
C
C     ---------------
C     LOCAL VARIABLES
C     ---------------
C
      INTEGER IWB0, IWB1, IWB2, IWTC, IWXC, IWYC, IWZC, JIN
C
C
C***FIRST EXECUTABLE STATEMENT  VY0
C
C     ... PARTITION WORK ARRAYS
C
      IWTC   = 1
      IWXC  = IWTC + M
      IWYC  = IWXC + M
      IWZC  = IWYC + M
      IWB0  = IWZC + M
      IWB1  = IWB0 + M
      IWB2  = IWB1 + M
C     Total = IB2  + M
C
      JIN   = 1
C     Total = JIN  + M
C
C     ... WY0 DOES ALL THE WORK
C
      CALL WY0(M,X,F,WORK(IWTC),WORK(IWXC),WORK(IWYC),WORK(IWZC),
     +         IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO)
C
      RETURN
      END

      SUBROUTINE WY0 (M, X, F, TCMP, XCMP, YCMP, ZCMP, INDX, B0, B1,
     +   B2, INFO)
C***BEGIN PROLOGUE  WY0
C***SUBSIDIARY
C***PURPOSE  Computes the Bessel function of the second kind
C            of order zero (Y0) for a vector of arguments
C***LIBRARY   VFNLIB
C***AUTHOR  SAUNDERS, B. V., (NIST)
C             COMPUTING AND APPLIED MATHEMATICS LABORATORY
C             NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C             GAITHERSBURG, MD 20899
C           BOISVERT, R. F., (NIST)
C             COMPUTING AND APPLIED MATHEMATICS LABORATORY
C             NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C             GAITHERSBURG, MD 20899
C***ROUTINES CALLED  IWCS, R1MACH, WNGT, WGTHR, WGTLE, WGT, WLE,
C                    WSCTR, WCS
C***REVISION HISTORY  (YYMMDD)
C   910226  DATE WRITTEN
C***END PROLOGUE  WY0
C
C  ----------
C  PARAMETERS
C  ----------
C
      INTEGER INFO, INDX, M
      REAL    B0, B1, B2, F, X, TCMP, XCMP, YCMP, ZCMP
C
      DIMENSION B0(M), B1(M), B2(M), F(M), INDX(M), X(M), TCMP(M),
     +          XCMP(M), YCMP(M), ZCMP(M)
C
C  ---------------
C  LOCAL VARIABLES
C  ---------------
C
      INTEGER LBJ0, LBM0, LBM02, LBTH0, LBT02, LBY0
      PARAMETER (LBJ0=13, LBM0=21, LBM02=2, LBTH0=24, LBT02=2,
     +           LBY0=13)
C
      INTEGER IWCS, J, JH, KEY, N, NA, NB, NTJ0, NTM0, NTM02, NTTH0,
     +        NTT02, NTY0
      REAL    BJ0CS, BM0CS, BM02CS, BTH0CS, BT02CS, BY0CS, C1,
     +         C2, EPMACH, EPS, R1MACH, PI4, TWODPI, XSML, XMAX
C
      DIMENSION BY0CS(LBY0), BM0CS(LBM0), BM02CS(LBM02), BTH0CS(LBTH0),
     +          BT02CS(LBT02), BJ0CS(LBJ0)
C
      SAVE BJ0CS, BM0CS, BTH0CS, BY0CS, NTJ0, NTM0, NTTH0, NTY0,
     +     PI4, TWODPI, XSML, XMAX
C
C----------------------------------------------------------------------
C
C Series for BY0        on the interval  0.          to  1.60000D+01
C                                        with weighted error   1.20E-17
C                                         log weighted error  16.92
C                               significant figures required  16.15
C                                    decimal places required  17.48
C
      DATA BY0 CS( 1) /   -.0112778393 92865573E0 /
      DATA BY0 CS( 2) /   -.1283452375 6042035E0 /
      DATA BY0 CS( 3) /   -.1043788479 9794249E0 /
      DATA BY0 CS( 4) /    .0236627491 83969695E0 /
      DATA BY0 CS( 5) /   -.0020903916 47700486E0 /
      DATA BY0 CS( 6) /    .0001039754 53939057E0 /
      DATA BY0 CS( 7) /   -.0000033697 47162423E0 /
      DATA BY0 CS( 8) /    .0000000772 93842676E0 /
      DATA BY0 CS( 9) /   -.0000000013 24976772E0 /
      DATA BY0 CS(10) /    .0000000000 17648232E0 /
      DATA BY0 CS(11) /   -.0000000000 00188105E0 /
      DATA BY0 CS(12) /    .0000000000 00001641E0 /
      DATA BY0 CS(13) /   -.0000000000 00000011E0 /
C
C----------------------------------------------------------------------
C
C Series for BM0        on the interval  0.          to  6.25000D-02
C                                        with weighted error   4.98E-17
C                                         log weighted error  16.30
C                               significant figures required  14.97
C                                    decimal places required  16.96
C
      DATA BM0 CS( 1) /    .0928496163 7381644E0 /
      DATA BM0 CS( 2) /   -.0014298770 7403484E0 /
      DATA BM0 CS( 3) /    .0000283057 9271257E0 /
      DATA BM0 CS( 4) /   -.0000014330 0611424E0 /
      DATA BM0 CS( 5) /    .0000001202 8628046E0 /
      DATA BM0 CS( 6) /   -.0000000139 7113013E0 /
      DATA BM0 CS( 7) /    .0000000020 4076188E0 /
      DATA BM0 CS( 8) /   -.0000000003 5399669E0 /
      DATA BM0 CS( 9) /    .0000000000 7024759E0 /
      DATA BM0 CS(10) /   -.0000000000 1554107E0 /
      DATA BM0 CS(11) /    .0000000000 0376226E0 /
      DATA BM0 CS(12) /   -.0000000000 0098282E0 /
      DATA BM0 CS(13) /    .0000000000 0027408E0 /
      DATA BM0 CS(14) /   -.0000000000 0008091E0 /
      DATA BM0 CS(15) /    .0000000000 0002511E0 /
      DATA BM0 CS(16) /   -.0000000000 0000814E0 /
      DATA BM0 CS(17) /    .0000000000 0000275E0 /
      DATA BM0 CS(18) /   -.0000000000 0000096E0 /
      DATA BM0 CS(19) /    .0000000000 0000034E0 /
      DATA BM0 CS(20) /   -.0000000000 0000012E0 /
      DATA BM0 CS(21) /    .0000000000 0000004E0 /
C
C----------------------------------------------------------------------
C
C Series for BM02       Used in double precision version only
C
      DATA BM02CS( 1) /   1.0E0 /
      DATA BM02CS( 2) /   0.0E0 /
C
C----------------------------------------------------------------------
C
C Series for BTH0       on the interval  0.          to  6.25000D-02
C                                        with weighted error   3.67E-17
C                                         log weighted error  16.44
C                               significant figures required  15.53
C                                    decimal places required  17.13
C
C
      DATA BTH0CS( 1) /   -.2463916377 4300119E0 /
      DATA BTH0CS( 2) /    .0017370983 07508963E0 /
      DATA BTH0CS( 3) /   -.0000621836 33402968E0 /
      DATA BTH0CS( 4) /    .0000043680 50165742E0 /
      DATA BTH0CS( 5) /   -.0000004560 93019869E0 /
      DATA BTH0CS( 6) /    .0000000621 97400101E0 /
      DATA BTH0CS( 7) /   -.0000000103 00442889E0 /
      DATA BTH0CS( 8) /    .0000000019 79526776E0 /
      DATA BTH0CS( 9) /   -.0000000004 28198396E0 /
      DATA BTH0CS(10) /    .0000000001 02035840E0 /
      DATA BTH0CS(11) /   -.0000000000 26363898E0 /
      DATA BTH0CS(12) /    .0000000000 07297935E0 /
      DATA BTH0CS(13) /   -.0000000000 02144188E0 /
      DATA BTH0CS(14) /    .0000000000 00663693E0 /
      DATA BTH0CS(15) /   -.0000000000 00215126E0 /
      DATA BTH0CS(16) /    .0000000000 00072659E0 /
      DATA BTH0CS(17) /   -.0000000000 00025465E0 /
      DATA BTH0CS(18) /    .0000000000 00009229E0 /
      DATA BTH0CS(19) /   -.0000000000 00003448E0 /
      DATA BTH0CS(20) /    .0000000000 00001325E0 /
      DATA BTH0CS(21) /   -.0000000000 00000522E0 /
      DATA BTH0CS(22) /    .0000000000 00000210E0 /
      DATA BTH0CS(23) /   -.0000000000 00000087E0 /
      DATA BTH0CS(24) /    .0000000000 00000036E0 /
C
C----------------------------------------------------------------------
C
C Series for BT02       Used in double precision version only
C
      DATA BT02CS( 1) /   1.0E0 /
      DATA BT02CS( 2) /   0.0E0 /
C
C----------------------------------------------------------------------
C
C Series for BJ0        on the interval  0.          to  1.60000D+01
C                                        with weighted error   7.47E-18
C                                         log weighted error  17.13
C                               significant figures required  16.98
C                                    decimal places required  17.68
C
      DATA BJ0 CS( 1) /    .1002541619 68939137E0 /
      DATA BJ0 CS( 2) /   -.6652230077 64405132E0 /
      DATA BJ0 CS( 3) /    .2489837034 98281314E0 /
      DATA BJ0 CS( 4) /   -.0332527231 700357697E0 /
      DATA BJ0 CS( 5) /    .0023114179 304694015E0 /
      DATA BJ0 CS( 6) /   -.0000991127 741995080E0 /
      DATA BJ0 CS( 7) /    .0000028916 708643998E0 /
      DATA BJ0 CS( 8) /   -.0000000612 108586630E0 /
      DATA BJ0 CS( 9) /    .0000000009 838650793E0 /
      DATA BJ0 CS(10) /   -.0000000000 124235515E0 /
      DATA BJ0 CS(11) /    .0000000000 001265433E0 /
      DATA BJ0 CS(12) /   -.0000000000 000010619E0 /
      DATA BJ0 CS(13) /    .0000000000 000000074E0 /
C
C----------------------------------------------------------------------
C
      DATA NTY0   / 0 /
C
C***FIRST EXECUTABLE STATEMENT  WY0
C
      IF (M .LE. 0)  GOTO 910
C
      IF (NTY0 .EQ. 0) THEN
         PI4 = ATAN(1.0E0)
         TWODPI = 0.50E0/PI4
         EPMACH = R1MACH(3)
         EPS = 0.10E0*EPMACH
         NTY0 = IWCS(BY0CS, LBY0, EPS)
         NTM0 = IWCS(BM0CS, LBM0, EPS)
         NTM02 = IWCS(BM02CS, LBM02, EPS)
         NTTH0 = IWCS(BTH0CS, LBTH0, EPS)
         NTT02 = IWCS(BT02CS, LBT02, EPS)
         NTJ0 = IWCS(BJ0CS, LBJ0, EPS)
         XSML = SQRT (4.0E0*EPMACH)
         XMAX = 1.0E0/R1MACH(4)
      ENDIF
C
      CALL WNLE(M,X,0.0E0,KEY)
      IF (KEY .NE. 0) GO TO 920
C
      CALL WNGT(M,X,XMAX,KEY)
      IF (KEY .NE. 0) GO TO 930
C
C  ----------------
C  CASE  X .LE. 4.0
C  ----------------
C
      CALL WLE(M,X,4.0E0,N,INDX)
      IF (N .GT. 0) THEN
         CALL WGTHR(N,X,INDX,XCMP)
C
C        ... COMPUTE J0(X) ... RESULT IN ZCMP
C
         DO 10 J=1,N
            TCMP(J) = 0.125E0*XCMP(J)**2 - 1.0E0
   10    CONTINUE
         CALL WCS(N,TCMP,BJ0CS,NTJ0,ZCMP,B0,B1,B2)
C
         CALL WCS(N,TCMP,BY0CS,NTY0,YCMP,B0,B1,B2)
         DO 30 J=1,N
            YCMP(J) = TWODPI*LOG(0.50E0*XCMP(J))*ZCMP(J)
     +                 + 0.375E0 + YCMP(J)
   30    CONTINUE
         CALL WSCTR(N,YCMP,INDX,F)
      ENDIF
C
C  ----------------
C  CASE  X .GT. 4.0
C  ----------------
C
      CALL WGT(M,X,4.0E0,N,INDX)
      IF (N .GT. 0) THEN
         CALL WGTHR(N,X,INDX,XCMP)
         DO 50 J=1,N
            TCMP(J) = 32.0E0/XCMP(J)**2 - 1.0E0
  50     CONTINUE
         CALL WCS(N,TCMP,BTH0CS,NTTH0,ZCMP,B0,B1,B2)
         CALL WCS(N,TCMP,BM0CS,NTM0,YCMP,B0,B1,B2)
         DO 60 J=1,N
            YCMP(J) = (0.750E0 + YCMP(J)) / SQRT(XCMP(J))
            ZCMP(J) = (XCMP(J) - PI4) + ZCMP(J) / XCMP(J)
            YCMP(J) = YCMP(J) * SIN(ZCMP(J))
  60     CONTINUE
         CALL WSCTR(N,YCMP,INDX,F)
      ENDIF
C
C  -----
C  EXITS
C  -----
C
C     ... NORMAL
C
      INFO = 0
      GO TO 999
C
C     ... M .LE. 0
C
  910 CONTINUE
      INFO = 1
      GO TO 999
C
C     ... X I ZERO OR NEGATIVE
C
  920 CONTINUE
      INFO = 2
      INDX(1) = KEY
      GO TO 999
C
C     ... NO PRECISION BECAUSE X BIG
C
  930 CONTINUE
      INFO = 3
      INDX(1) = KEY
      GO TO 999
C
  999 CONTINUE
      RETURN
      END