SUBROUTINE VK0 (M, X, F, WORK, IWORK, INFO) C***BEGIN PROLOGUE VK0 C***PURPOSE Computes the hyperbolic Bessel function of the third kind C of order zero (K0) for a vector of real arguments C***LIBRARY VFNLIB C***CATEGORY C10B1 C***TYPE SINGLE PRECISION (VK0-S, DVK0-D) C***KEYWORDS BESSEL FUNCTION,THIRD KIND,HYPERBOLIC BESSEL FUNCTION, C MODIFIED BESSEL FUNCTION, 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 VK0 computes the modified (hyperbolic) Bessel function of the C third kind of order zero (K0) for real arguments using uniform C approximation by Chebyshev polynomials. C C C P A R A M E T E R S C C M (Input) Integer 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 -1 Yes Warning: Some X(i) so big K0 underflows. C The corresponding F(i) are set to zero. C 0 Yes Successfull execution. C 1 No Error: M .LE. 0 C 2 No Error: Some X(i) is zero or negative. 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 BESK0 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 WK0 C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C 930203 Minor modifications to prologue. C***END PROLOGUE VK0 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***FIRST EXECUTABLE STATEMENT VK0 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 ... WK0 DOES ALL THE WORK C CALL WK0(M,X,F,WORK(IWTC),WORK(IWXC),WORK(IWYC),WORK(IWZC), + IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO) C RETURN END SUBROUTINE WK0 (M, X, F, TCMP, XCMP, YCMP, ZCMP, INDX, B0, B1, + B2, INFO) C***BEGIN PROLOGUE WK0 C***SUBSIDIARY C***PURPOSE Computes the hyperbolic Bessel function of the third kind C of order zero (K0) 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, WNLE, WGTHR, WGTLE, WGT, WLE, C WSCTR, WCS C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WK0 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 LAK0, LAK02, LBI0, LBK0 PARAMETER ( LAK0=17, LAK02=14, LBI0=12, LBK0=11 ) C INTEGER I, IWCS, J, KEY, N, NTI0, NTK0, NTAK0, NTAK02, NTOT REAL AK0CS, AK02CS, BI0CS, BK0CS, EPMACH, EPS, R1MACH, + XMAX, XSML C DIMENSION AK0CS(LAK0), AK02CS(LAK02), BI0CS(LBI0), BK0CS(LBK0) C SAVE BK0CS, AK0CS, AK02CS, BI0CS, N, NTAK0, NTAK02, NTI0, NTK0, + XMAX, XSML C C---------------------------------------------------------------------- C C Series for BK0 on the interval 0. to 4.00000D+00 C with weighted error 3.57E-19 C log weighted error 18.45 C significant figures required 17.99 C decimal places required 18.97 C DATA BK0 CS( 1) / -.0353273932 3390276872E0 / DATA BK0 CS( 2) / .3442898999 246284869E0 / DATA BK0 CS( 3) / .0359799365 1536150163E0 / DATA BK0 CS( 4) / .0012646154 1144692592E0 / DATA BK0 CS( 5) / .0000228621 2103119451E0 / DATA BK0 CS( 6) / .0000002534 7910790261E0 / DATA BK0 CS( 7) / .0000000019 0451637722E0 / DATA BK0 CS( 8) / .0000000000 1034969525E0 / DATA BK0 CS( 9) / .0000000000 0004259816E0 / DATA BK0 CS(10) / .0000000000 0000013744E0 / DATA BK0 CS(11) / .0000000000 0000000035E0 / C C---------------------------------------------------------------------- C C Series for AK0 on the interval 1.25000D-01 to 5.00000D-01 C with weighted error 5.34E-17 C log weighted error 16.27 C significant figures required 14.92 C decimal places required 16.89 C DATA AK0 CS( 1) / -.0764394790 3327941E0 / DATA AK0 CS( 2) / -.0223565260 5699819E0 / DATA AK0 CS( 3) / .0007734181 1546938E0 / DATA AK0 CS( 4) / -.0000428100 6688886E0 / DATA AK0 CS( 5) / .0000030817 0017386E0 / DATA AK0 CS( 6) / -.0000002639 3672220E0 / DATA AK0 CS( 7) / .0000000256 3713036E0 / DATA AK0 CS( 8) / -.0000000027 4270554E0 / DATA AK0 CS( 9) / .0000000003 1694296E0 / DATA AK0 CS(10) / -.0000000000 3902353E0 / DATA AK0 CS(11) / .0000000000 0506804E0 / DATA AK0 CS(12) / -.0000000000 0068895E0 / DATA AK0 CS(13) / .0000000000 0009744E0 / DATA AK0 CS(14) / -.0000000000 0001427E0 / DATA AK0 CS(15) / .0000000000 0000215E0 / DATA AK0 CS(16) / -.0000000000 0000033E0 / DATA AK0 CS(17) / .0000000000 0000005E0 / C C---------------------------------------------------------------------- C C Series for AK02 on the interval 0. to 1.25000D-01 C with weighted error 2.34E-17 C log weighted error 16.63 C significant figures required 14.67 C decimal places required 17.20 C DATA AK02CS( 1) / -.0120186982 6307592E0 / DATA AK02CS( 2) / -.0091748526 9102569E0 / DATA AK02CS( 3) / .0001444550 9317750E0 / DATA AK02CS( 4) / -.0000040136 1417543E0 / DATA AK02CS( 5) / .0000001567 8318108E0 / DATA AK02CS( 6) / -.0000000077 7011043E0 / DATA AK02CS( 7) / .0000000004 6111825E0 / DATA AK02CS( 8) / -.0000000000 3158592E0 / DATA AK02CS( 9) / .0000000000 0243501E0 / DATA AK02CS(10) / -.0000000000 0020743E0 / DATA AK02CS(11) / .0000000000 0001925E0 / DATA AK02CS(12) / -.0000000000 0000192E0 / DATA AK02CS(13) / .0000000000 0000020E0 / DATA AK02CS(14) / -.0000000000 0000002E0 / C C---------------------------------------------------------------------- C C Series for BI0 on the interval 0. to 9.00000D+00 C with weighted error 2.46E-18 C log weighted error 17.61 C significant figures required 17.90 C decimal places required 18.15 C DATA BI0 CS( 1) / -.0766054725 2839144951E0 / DATA BI0 CS( 2) / 1.9273379539 93808270E0 / DATA BI0 CS( 3) / .2282644586 920301339E0 / DATA BI0 CS( 4) / .0130489146 6707290428E0 / DATA BI0 CS( 5) / .0004344270 9008164874E0 / DATA BI0 CS( 6) / .0000094226 5768600193E0 / DATA BI0 CS( 7) / .0000001434 0062895106E0 / DATA BI0 CS( 8) / .0000000016 1384906966E0 / DATA BI0 CS( 9) / .0000000000 1396650044E0 / DATA BI0 CS(10) / .0000000000 0009579451E0 / DATA BI0 CS(11) / .0000000000 0000053339E0 / DATA BI0 CS(12) / .0000000000 0000000245E0 / C C---------------------------------------------------------------------- C DATA NTK0 / 0 / C C***FIRST EXECUTABLE STATEMENT WK0 C IF (M .LE. 0) GO TO 910 C IF (NTK0 .EQ. 0) THEN EPMACH = R1MACH(3) EPS = 0.10E0*EPMACH NTK0 = IWCS(BK0CS, LBK0, EPS) NTAK0 = IWCS(AK0CS, LAK0, EPS) NTAK02 = IWCS(AK02CS, LAK02, EPS) NTI0 = IWCS(BI0CS, LBI0, EPS) XSML = SQRT (4.0E0*EPMACH) XMAX = -LOG(R1MACH(1)) XMAX = XMAX - 0.50E0*XMAX*LOG(XMAX)/(XMAX+0.50E0) - 0.010E0 ENDIF C NTOT = 0 C CALL WNLE(M,X,0.0E0,KEY) IF (KEY .NE. 0) GO TO 920 C C ------------------ C CASE X .GT. XMAX C ------------------ C C NOTE -- K0 UNDERFLOWS FOR X .GT. XMAX C DO 5 I=1,M F(I) = 0.0E0 5 CONTINUE C C ---------------- C CASE X .LE. 2.0 C ---------------- C CALL WLE(M,X,2.0E0,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL WGTHR(N,X,INDX,XCMP) C C ... COMPUTE I0(X) ... RESULT IN ZCMP C DO 10 J=1,N TCMP(J) = XCMP(J)**2/4.50E0 - 1.0E0 10 CONTINUE CALL WCS(N,TCMP,BI0CS,NTI0,ZCMP,B0,B1,B2) DO 15 J=1,N ZCMP(J) = 2.750E0 + ZCMP(J) 15 CONTINUE C DO 20 J=1,N TCMP(J) = 0.50E0*XCMP(J)**2 - 1.0E0 20 CONTINUE CALL WCS(N,TCMP,BK0CS,NTK0,YCMP,B0,B1,B2) DO 30 J=1,N YCMP(J) = -LOG(0.50E0*XCMP(J))*ZCMP(J) - 0.250E0 + YCMP(J) 30 CONTINUE CALL WSCTR(N,YCMP,INDX,F) ENDIF C C ------------------------- C CASE 2.0 .LT. X .LE. 8.0 C ------------------------- C CALL WGTLE(M,X,2.0E0,8.0E0,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL WGTHR(N,X,INDX,XCMP) DO 40 J=1,N TCMP(J) = (16.0E0/XCMP(J) - 5.0E0)/3.0E0 40 CONTINUE CALL WCS(N,TCMP,AK0CS,NTAK0,ZCMP,B0,B1,B2) DO 50 J=1,N YCMP(J) = EXP(-XCMP(J))*(1.250E0 + ZCMP(J))/SQRT(XCMP(J)) 50 CONTINUE CALL WSCTR(N,YCMP,INDX,F) ENDIF C C -------------------------- C CASE 8.0 .LT. X .LE. XMAX C -------------------------- C CALL WGTLE(M,X,8.0E0,XMAX,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL WGTHR(N,X,INDX,XCMP) DO 60 J=1,N TCMP(J) = 16.0E0/XCMP(J) - 1.0E0 60 CONTINUE CALL WCS(N,TCMP,AK02CS,NTAK02,ZCMP,B0,B1,B2) DO 70 J=1,N YCMP(J) = EXP(-XCMP(J))*(1.250E0 + ZCMP(J))/SQRT(XCMP(J)) 70 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 IF (NTOT .NE. M) INFO = -1 GO TO 999 C C ... X IS ZERO OR NEGATIVE C 920 CONTINUE INFO = 2 INDX(1) = KEY GO TO 999 C 999 CONTINUE RETURN END