SUBROUTINE VK1 (M, X, F, WORK, IWORK, INFO) C***BEGIN PROLOGUE VK1 C***PURPOSE Computes the hyperbolic Bessel function of the third kind C of order one (K1) for a vector of real arguments C***LIBRARY VFNLIB C***CATEGORY C10B1 C***TYPE SINGLE PRECISION (VK1-S, DVK1-D) C***KEYWORDS BESSEL FUNCTION,THIRD KIND,HYPERBOLIC BESSEL FUNCTION, C MODIFIED BESSEL FUNCTION, 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 VK1 computes the modified (hyperbolic) Bessel function of the C third kind of order one (K1) 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 K1 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 3 No Error: Some X(i) so small K1 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 BESK1 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 WK1 C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C 930203 Minor modifications to prologue. C***END PROLOGUE VK1 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 VK1 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 ... WK1 DOES ALL THE WORK C CALL WK1(M,X,F,WORK(IWTC),WORK(IWXC),WORK(IWYC),WORK(IWZC), + IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO) C RETURN END SUBROUTINE WK1 (M, X, F, TCMP, XCMP, YCMP, ZCMP, INDX, B0, B1, + B2, INFO) C***BEGIN PROLOGUE WK1 C***SUBSIDIARY C***PURPOSE Computes the hyperbolic Bessel function of the third kind C of order one (K1) 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 WK1 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 LAK1, LAK12, LBI1, LBK1 PARAMETER ( LAK1=17, LAK12=14, LBI1=11, LBK1=11 ) C INTEGER I, IWCS, J, KEY, N, NTI1, NTK1, NTAK1, NTAK12, NTOT REAL AK1CS, AK12CS, BI1CS, BK1CS, EPMACH, EPS, R1MACH, + XMAX, XMIN, XSML C DIMENSION AK1CS(LAK1), AK12CS(LAK12), BI1CS(LBI1), BK1CS(LBK1) C SAVE AK1CS, AK12CS, BI1CS, BK1CS, NTAK1, NTAK12, NTI1, NTK1, + XMAX, XMIN, XSML C C---------------------------------------------------------------------- C C Series for BK1 on the interval 0. to 4.00000D+00 C with weighted error 7.02E-18 C log weighted error 17.15 C significant figures required 16.73 C decimal places required 17.67 C DATA BK1 CS( 1) / .0253002273 389477705E0 / DATA BK1 CS( 2) / -.3531559607 76544876E0 / DATA BK1 CS( 3) / -.1226111808 22657148E0 / DATA BK1 CS( 4) / -.0069757238 596398643E0 / DATA BK1 CS( 5) / -.0001730288 957513052E0 / DATA BK1 CS( 6) / -.0000024334 061415659E0 / DATA BK1 CS( 7) / -.0000000221 338763073E0 / DATA BK1 CS( 8) / -.0000000001 411488392E0 / DATA BK1 CS( 9) / -.0000000000 006666901E0 / DATA BK1 CS(10) / -.0000000000 000024274E0 / DATA BK1 CS(11) / -.0000000000 000000070E0 / C C---------------------------------------------------------------------- C C Series for BI1 ON THE INTERVAL 0. to 9.00000D+00 C WITH WEIGHTED ERROR 2.40E-17 C LOG WEIGHTED ERROR 16.62 C SIGNIFICANT FIGURES REQUIRED 16.23 C DECIMAL PLACES REQUIRED 17.14 C DATA BI1 CS( 1) / -.0019717132 61099859E0 / DATA BI1 CS( 2) / .4073488766 7546481E0 / DATA BI1 CS( 3) / .0348389942 99959456E0 / DATA BI1 CS( 4) / .0015453945 56300123E0 / DATA BI1 CS( 5) / .0000418885 21098377E0 / DATA BI1 CS( 6) / .0000007649 02676483E0 / DATA BI1 CS( 7) / .0000000100 42493924E0 / DATA BI1 CS( 8) / .0000000000 99322077E0 / DATA BI1 CS( 9) / .0000000000 00766380E0 / DATA BI1 CS(10) / .0000000000 00004741E0 / DATA BI1 CS(11) / .0000000000 00000024E0 / C C---------------------------------------------------------------------- C C Series for AK1 on the interval 1.25000D-01 to 5.00000D-01 C with weighted error 6.06E-17 C log weighted error 16.22 C significant figures required 15.41 C decimal places required 16.83 C DATA AK1 CS( 1) / .2744313406 973883E0 / DATA AK1 CS( 2) / .0757198995 3199368E0 / DATA AK1 CS( 3) / -.0014410515 5647540E0 / DATA AK1 CS( 4) / .0000665011 6955125E0 / DATA AK1 CS( 5) / -.0000043699 8470952E0 / DATA AK1 CS( 6) / .0000003540 2774997E0 / DATA AK1 CS( 7) / -.0000000331 1163779E0 / DATA AK1 CS( 8) / .0000000034 4597758E0 / DATA AK1 CS( 9) / -.0000000003 8989323E0 / DATA AK1 CS(10) / .0000000000 4720819E0 / DATA AK1 CS(11) / -.0000000000 0604783E0 / DATA AK1 CS(12) / .0000000000 0081284E0 / DATA AK1 CS(13) / -.0000000000 0011386E0 / DATA AK1 CS(14) / .0000000000 0001654E0 / DATA AK1 CS(15) / -.0000000000 0000248E0 / DATA AK1 CS(16) / .0000000000 0000038E0 / DATA AK1 CS(17) / -.0000000000 0000006E0 / C C---------------------------------------------------------------------- C C Series for AK12 on the interval 0. to 1.25000D-01 C with weighted error 2.58E-17 C log weighted error 16.59 C significant figures required 15.22 C decimal places required 17.16 C DATA AK12CS( 1) / .0637930834 3739001E0 / DATA AK12CS( 2) / .0283288781 3049721E0 / DATA AK12CS( 3) / -.0002475370 6739052E0 / DATA AK12CS( 4) / .0000057719 7245160E0 / DATA AK12CS( 5) / -.0000002068 9392195E0 / DATA AK12CS( 6) / .0000000097 3998344E0 / DATA AK12CS( 7) / -.0000000005 5853361E0 / DATA AK12CS( 8) / .0000000000 3732996E0 / DATA AK12CS( 9) / -.0000000000 0282505E0 / DATA AK12CS(10) / .0000000000 0023720E0 / DATA AK12CS(11) / -.0000000000 0002176E0 / DATA AK12CS(12) / .0000000000 0000215E0 / DATA AK12CS(13) / -.0000000000 0000022E0 / DATA AK12CS(14) / .0000000000 0000002E0 / C C---------------------------------------------------------------------- C DATA NTK1 / 0 / C C***FIRST EXECUTABLE STATEMENT WK1 C IF (M .LE. 0) GO TO 910 C IF (NTK1 .EQ. 0) THEN EPMACH = R1MACH(3) EPS = 0.10E0*EPMACH NTK1 = IWCS(BK1CS, LBK1, EPS) NTI1 = IWCS(BI1CS, LBI1, EPS) NTAK1 = IWCS(AK1CS, LAK1, EPS) NTAK12 = IWCS(AK12CS, LAK12, EPS) XMIN = EXP(MAX( LOG(R1MACH(1)), -LOG(R1MACH(2))) + 0.010E0) XSML = SQRT(4.0E0*EPMACH) XMAX = -LOG(R1MACH(1)) XMAX = XMAX - 0.50E0*XMAX*LOG(XMAX)/(XMAX + 0.50E0) ENDIF C NTOT = 0 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 .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 I1(X) ... RESULT IN ZCMP C DO 20 J=1,N TCMP(J) = XCMP(J)**2/4.50E0 - 1.0E0 20 CONTINUE CALL WCS(N,TCMP,BI1CS,NTI1,ZCMP,B0,B1,B2) DO 30 J=1,N ZCMP(J) = XCMP(J)*(0.8750E0 + ZCMP(J)) 30 CONTINUE C DO 40 J=1,N TCMP(J) = 0.50E0*XCMP(J)**2 - 1.0E0 40 CONTINUE CALL WCS(N,TCMP,BK1CS,NTK1,YCMP,B0,B1,B2) DO 50 J=1,N ZCMP(J) = LOG(0.50E0*XCMP(J))*ZCMP(J) + + (0.750E0 + YCMP(J))/XCMP(J) 50 CONTINUE CALL WSCTR(N,ZCMP,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 60 J=1,N TCMP(J) = (16.0E0/XCMP(J) - 5.0E0)/3.0E0 60 CONTINUE CALL WCS(N,TCMP,AK1CS,NTAK1,YCMP,B0,B1,B2) DO 70 J=1,N ZCMP(J) = EXP(-XCMP(J))*(1.250E0 + YCMP(J))/SQRT(XCMP(J)) 70 CONTINUE CALL WSCTR(N,ZCMP,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 80 J=1,N TCMP(J) = 16.0E0/XCMP(J) - 1.0E0 80 CONTINUE CALL WCS(N,TCMP,AK12CS,NTAK12,YCMP,B0,B1,B2) DO 90 J=1,N ZCMP(J) = EXP(-XCMP(J))*(1.250E0 + YCMP(J))/SQRT(XCMP(J)) 90 CONTINUE CALL WSCTR(N,ZCMP,INDX,F) ENDIF C C ----- C EXITS C ----- C C ... NORMAL C INFO = 0 IF (NTOT .NE. M) INFO = -1 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 K1 OVERFLOWS C 930 CONTINUE INFO = 3 INDX(1) = KEY GO TO 999 C 999 CONTINUE RETURN END