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