SUBROUTINE VY1 (M, X, F, WORK, IWORK, INFO) C***BEGIN PROLOGUE VY1 C***PURPOSE Computes the Bessel function of the second kind C of order one (Y1) for a vector of real arguments C***LIBRARY VFNLIB C***CATEGORY C10A1 C***TYPE SINGLE PRECISION (VY1-S, DVY1-D) C***KEYWORDS BESSEL FUNCTION,SECOND KIND, ORDER ONE, 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 VY1 computes the Bessel function of the second kind of order C one (Y1) 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 small Y1 overflows. C The index of the first offending argument is C returned in IWORK(1). C C C ********************************************************************* C This routine is a modification of the function BESY1 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 WY1 C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C 930203 Minor modifications to prologue. C***END PROLOGUE VY1 C C ---------- C PARAMETERS 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 VY1 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 ... WY1 DOES ALL THE WORK C CALL WY1(M,X,F,WORK(IWTC),WORK(IWXC),WORK(IWYC),WORK(IWZC), + IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO) C RETURN END SUBROUTINE WY1 (M, X, F, TCMP, XCMP, YCMP, ZCMP, INDX, B0, B1, + B2, INFO) C***BEGIN PROLOGUE WY1 C***SUBSIDIARY C***PURPOSE Computes the Bessel function of the second kind C of order one (Y1) 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, WNLE C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WY1 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 LBJ1, LBM1, LBM12, LBTH1, LBT12, LBY1 PARAMETER (LBJ1=12, LBM1=21, LBM12=2, LBTH1=24, LBT12=2, + LBY1=14) C INTEGER IWCS, J, JH, KEY, N, NA, NB, NTJ1, NTM1, NTM12, NTTH1, + NTT12, NTY1 REAL BJ1CS, BM1CS, BM12CS, BTH1CS, BT12CS, BY1CS, + C1, C2, EPMACH, EPS, PI4, R1MACH, TPI4, TWODPI, XMIN, + XSML, XMAX C DIMENSION BJ1CS(LBJ1), BM1CS(LBM1), BM12CS(LBM12), BTH1CS(LBTH1), + BT12CS(LBT12), BY1CS(LBY1) C SAVE BJ1CS, BM1CS, BTH1CS, BY1CS, NTJ1, NTM1, NTM12, NTTH1, NTT12, + NTY1, PI4, TPI4, TWODPI, XMIN, XSML, XMAX C C---------------------------------------------------------------------- C C C Series for BY1 on the interval 0. to 1.60000D+01 C with weighted error 1.87E-18 C log weighted error 17.73 C significant figures required 17.83 C decimal places required 18.30 C DATA BY1 CS( 1) / .0320804710 0611908629E0 / DATA BY1 CS( 2) / 1.2627078974 33500450E0 / DATA BY1 CS( 3) / .0064999618 9992317500E0 / DATA BY1 CS( 4) / -.0893616452 8860504117E0 / DATA BY1 CS( 5) / .0132508812 2175709545E0 / DATA BY1 CS( 6) / -.0008979059 1196483523E0 / DATA BY1 CS( 7) / .0000364736 1487958306E0 / DATA BY1 CS( 8) / -.0000010013 7438166600E0 / DATA BY1 CS( 9) / .0000000199 4539657390E0 / DATA BY1 CS(10) / -.0000000003 0230656018E0 / DATA BY1 CS(11) / .0000000000 0360987815E0 / DATA BY1 CS(12) / -.0000000000 0003487488E0 / DATA BY1 CS(13) / .0000000000 0000027838E0 / DATA BY1 CS(14) / -.0000000000 0000000186E0 / C C---------------------------------------------------------------------- C C Series for BM1 on the interval 0. to 6.25000D-02 C with weighted error 5.61E-17 C log weighted error 16.25 C significant figures required 14.97 C decimal places required 16.91 C DATA BM1 CS( 1) / .1047362510 931285E0 / DATA BM1 CS( 2) / .0044244389 3702345E0 / DATA BM1 CS( 3) / -.0000566163 9504035E0 / DATA BM1 CS( 4) / .0000023134 9417339E0 / DATA BM1 CS( 5) / -.0000001737 7182007E0 / DATA BM1 CS( 6) / .0000000189 3209930E0 / DATA BM1 CS( 7) / -.0000000026 5416023E0 / DATA BM1 CS( 8) / .0000000004 4740209E0 / DATA BM1 CS( 9) / -.0000000000 8691795E0 / DATA BM1 CS(10) / .0000000000 1891492E0 / DATA BM1 CS(11) / -.0000000000 0451884E0 / DATA BM1 CS(12) / .0000000000 0116765E0 / DATA BM1 CS(13) / -.0000000000 0032265E0 / DATA BM1 CS(14) / .0000000000 0009450E0 / DATA BM1 CS(15) / -.0000000000 0002913E0 / DATA BM1 CS(16) / .0000000000 0000939E0 / DATA BM1 CS(17) / -.0000000000 0000315E0 / DATA BM1 CS(18) / .0000000000 0000109E0 / DATA BM1 CS(19) / -.0000000000 0000039E0 / DATA BM1 CS(20) / .0000000000 0000014E0 / DATA BM1 CS(21) / -.0000000000 0000005E0 / C C---------------------------------------------------------------------- C C Series for BM12 Used in double precision version only C DATA BM12CS( 1) / 1.0E0 / DATA BM12CS( 2) / 0.0E0 / C C---------------------------------------------------------------------- C C Series for BTH1 on the interval 0. to 6.25000D-02 C with weighted error 4.10E-17 C log weighted error 16.39 C significant figures required 15.96 C decimal places required 17.08 C DATA BTH1CS( 1) / .7406014102 6313850E0 / DATA BTH1CS( 2) / -.0045717556 59637690E0 / DATA BTH1CS( 3) / .0001198185 10964326E0 / DATA BTH1CS( 4) / -.0000069645 61891648E0 / DATA BTH1CS( 5) / .0000006554 95621447E0 / DATA BTH1CS( 6) / -.0000000840 66228945E0 / DATA BTH1CS( 7) / .0000000133 76886564E0 / DATA BTH1CS( 8) / -.0000000024 99565654E0 / DATA BTH1CS( 9) / .0000000005 29495100E0 / DATA BTH1CS(10) / -.0000000001 24135944E0 / DATA BTH1CS(11) / .0000000000 31656485E0 / DATA BTH1CS(12) / -.0000000000 08668640E0 / DATA BTH1CS(13) / .0000000000 02523758E0 / DATA BTH1CS(14) / -.0000000000 00775085E0 / DATA BTH1CS(15) / .0000000000 00249527E0 / DATA BTH1CS(16) / -.0000000000 00083773E0 / DATA BTH1CS(17) / .0000000000 00029205E0 / DATA BTH1CS(18) / -.0000000000 00010534E0 / DATA BTH1CS(19) / .0000000000 00003919E0 / DATA BTH1CS(20) / -.0000000000 00001500E0 / DATA BTH1CS(21) / .0000000000 00000589E0 / DATA BTH1CS(22) / -.0000000000 00000237E0 / DATA BTH1CS(23) / .0000000000 00000097E0 / DATA BTH1CS(24) / -.0000000000 00000040E0 / C C---------------------------------------------------------------------- C C Series for BT12 Used in double precision version only C DATA BT12CS( 1) / 1.0E0 / DATA BT12CS( 2) / 0.0E0 / C C---------------------------------------------------------------------- C C Series for BJ1 on the interval 0. to 1.60000D+01 C with weighted error 4.48E-17 C log weighted error 16.35 C significant figures required 15.77 C decimal places required 16.89 C DATA BJ1 CS( 1) / -.1172614151 3332787E0 / DATA BJ1 CS( 2) / -.2536152183 0790640E0 / DATA BJ1 CS( 3) / .0501270809 84469569E0 / DATA BJ1 CS( 4) / -.0046315148 09625081E0 / DATA BJ1 CS( 5) / .0002479962 29415914E0 / DATA BJ1 CS( 6) / -.0000086789 48686278E0 / DATA BJ1 CS( 7) / .0000002142 93917143E0 / DATA BJ1 CS( 8) / -.0000000039 36093079E0 / DATA BJ1 CS( 9) / .0000000000 55911823E0 / DATA BJ1 CS(10) / -.0000000000 00632761E0 / DATA BJ1 CS(11) / .0000000000 00005840E0 / DATA BJ1 CS(12) / -.0000000000 00000044E0 / C C---------------------------------------------------------------------- C DATA NTY1 / 0 / C C***FIRST EXECUTABLE STATEMENT WY1 C IF (M .LE. 0) GO TO 910 C IF (NTY1 .EQ. 0) THEN EPMACH = R1MACH(3) EPS = 0.10E0*EPMACH NTY1 = IWCS(BY1CS, LBY1, EPS) NTM1 = IWCS(BM1CS, LBM1, EPS) NTM12 = IWCS(BM12CS, LBM12, EPS) NTTH1 = IWCS(BTH1CS, LBTH1, EPS) NTT12 = IWCS(BT12CS, LBT12, EPS) NTJ1 = IWCS(BJ1CS, LBJ1, EPS) XMIN = EXP(MAX( LOG(R1MACH(1)), -LOG(R1MACH(2))) + 0.010E0) XSML = SQRT(4.0E0*EPMACH) XMAX = 1.0/R1MACH(4) PI4 = ATAN(1.0E0) TPI4 = 3.0E0*PI4 TWODPI = 0.50E0/PI4 ENDIF C CALL WNLE(M,X,0.0E0,KEY) IF (KEY .NE. 0) GO TO 920 C CALL WNLE(M,X,XMIN,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 J1(X) ... RESULT IN ZCMP C DO 20 J=1,N TCMP(J) = 0.125E0*XCMP(J)**2 - 1.0E0 20 CONTINUE CALL WCS(N,TCMP,BJ1CS,NTJ1,ZCMP,B0,B1,B2) DO 30 J=1,N ZCMP(J) = XCMP(J)*(0.250E0 + ZCMP(J)) 30 CONTINUE C CALL WCS(N,TCMP,BY1CS,NTY1,YCMP,B0,B1,B2) DO 40 J=1,N ZCMP(J) = TWODPI*LOG(0.50E0*XCMP(J))*ZCMP(J) + + (0.50E0 + YCMP(J))/XCMP(J) 40 CONTINUE CALL WSCTR(N,ZCMP,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,BTH1CS,NTTH1,ZCMP,B0,B1,B2) CALL WCS(N,TCMP,BM1CS,NTM1,YCMP,B0,B1,B2) DO 60 J=1,N YCMP(J) = (0.750E0 + YCMP(J)) / SQRT(XCMP(J)) ZCMP(J) = (XCMP(J) - TPI4) + 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 IS ZERO OR NEGATIVE C 920 CONTINUE INFO = 2 INDX(1) = KEY GO TO 999 C C ... X SO SMALL THAT Y1 OVERFLOWS C 930 CONTINUE INFO = 3 INDX(1) = KEY GO TO 999 C 999 CONTINUE RETURN END