SUBROUTINE VI0 (M, X, F, WORK, IWORK, INFO) C***BEGIN PROLOGUE VI0 C***PURPOSE Computes the hyperbolic Bessel function of the first kind C of order zero (I0) for a vector of real arguments C***LIBRARY VFNLIB C***CATEGORY C10B1 C***TYPE SINGLE PRECISION (VI0-S, DVI0-D) C***KEYWORDS BESSEL FUNCTION,FIRST 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 VI0 computes the modified (hyperbolic) Bessel function of the C first kind of order zero (I0) 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 (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: Some abs(X(i)) so big I0 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 BESI0 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 WI0 C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C 930203 Minor modifications to prologue. C***END PROLOGUE VI0 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, IWY, IWTC, IWYC, IWZC, JIN C C***FIRST EXECUTABLE STATEMENT VI0 C C ... PARTITION WORK ARRAYS C IWY = 1 IWTC = IWY + M IWYC = IWTC + M IWZC = IWYC + M IWB0 = IWZC + M IWB1 = IWB0 + M IWB2 = IWB1 + M C Total = IWB2 + M C JIN = 1 C Total = JIN + M C C ... WI0 DOES ALL THE WORK C CALL WI0(M,X,F,WORK(IWY),WORK(IWTC),WORK(IWYC),WORK(IWZC), + IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO) C RETURN END SUBROUTINE WI0 (M, X, F, Y, TCMP, YCMP, ZCMP, INDX, B0, B1, B2, + INFO) C***BEGIN PROLOGUE WI0 C***SUBSIDIARY C***PURPOSE Computes the hyperbolic Bessel function of the first kind C of order zero (I0) 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 WI0 C C ---------- C PARAMETERS C ---------- C INTEGER INFO, INDX, M REAL B0, B1, B2, F, X, Y, TCMP, YCMP, ZCMP C DIMENSION B0(M), B1(M), B2(M), F(M), INDX(M), X(M), Y(M), + TCMP(M), YCMP(M), ZCMP(M) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER LAI0, LAI02, LBI0 PARAMETER ( LAI0=21, LAI02=22, LBI0=12 ) C INTEGER I, IWCS, J, KEY, N, NTI0, NTAI0, NTAI02 REAL AI0CS, AI02CS, BI0CS, EPMACH, EPS, R1MACH, XSML, + XMAX C DIMENSION AI0CS(LAI0), AI02CS(LAI02), BI0CS(LBI0) C SAVE AI0CS, AI02CS, BI0CS, N, NTAI0, NTAI02, NTI0, XSML, XMAX 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 C Series for AI0 on the interval 1.25000D-01 to 3.33333D-01 C with weighted error 7.87E-17 C log weighted error 16.10 C significant figures required 14.69 C decimal places required 16.76 C DATA AI0 CS( 1) / .0757599449 4023796E0 / DATA AI0 CS( 2) / .0075913808 1082334E0 / DATA AI0 CS( 3) / .0004153131 3389237E0 / DATA AI0 CS( 4) / .0000107007 6463439E0 / DATA AI0 CS( 5) / -.0000079011 7997921E0 / DATA AI0 CS( 6) / -.0000007826 1435014E0 / DATA AI0 CS( 7) / .0000002783 8499429E0 / DATA AI0 CS( 8) / .0000000082 5247260E0 / DATA AI0 CS( 9) / -.0000000120 4463945E0 / DATA AI0 CS(10) / .0000000015 5964859E0 / DATA AI0 CS(11) / .0000000002 2925563E0 / DATA AI0 CS(12) / -.0000000001 1916228E0 / DATA AI0 CS(13) / .0000000000 1757854E0 / DATA AI0 CS(14) / .0000000000 0112822E0 / DATA AI0 CS(15) / -.0000000000 0114684E0 / DATA AI0 CS(16) / .0000000000 0027155E0 / DATA AI0 CS(17) / -.0000000000 0002415E0 / DATA AI0 CS(18) / -.0000000000 0000608E0 / DATA AI0 CS(19) / .0000000000 0000314E0 / DATA AI0 CS(20) / -.0000000000 0000071E0 / DATA AI0 CS(21) / .0000000000 0000007E0 / C C---------------------------------------------------------------------- C C Series for AI02 on the interval 0. to 1.25000D-01 C with weighted error 3.79E-17 C log weighted error 16.42 C significant figures required 14.86 C decimal places required 17.09 C DATA AI02CS( 1) / .0544904110 1410882E0 / DATA AI02CS( 2) / .0033691164 7825569E0 / DATA AI02CS( 3) / .0000688975 8346918E0 / DATA AI02CS( 4) / .0000028913 7052082E0 / DATA AI02CS( 5) / .0000002048 9185893E0 / DATA AI02CS( 6) / .0000000226 6668991E0 / DATA AI02CS( 7) / .0000000033 9623203E0 / DATA AI02CS( 8) / .0000000004 9406022E0 / DATA AI02CS( 9) / .0000000000 1188914E0 / DATA AI02CS(10) / -.0000000000 3149915E0 / DATA AI02CS(11) / -.0000000000 1321580E0 / DATA AI02CS(12) / -.0000000000 0179419E0 / DATA AI02CS(13) / .0000000000 0071801E0 / DATA AI02CS(14) / .0000000000 0038529E0 / DATA AI02CS(15) / .0000000000 0001539E0 / DATA AI02CS(16) / -.0000000000 0004151E0 / DATA AI02CS(17) / -.0000000000 0000954E0 / DATA AI02CS(18) / .0000000000 0000382E0 / DATA AI02CS(19) / .0000000000 0000176E0 / DATA AI02CS(20) / -.0000000000 0000034E0 / DATA AI02CS(21) / -.0000000000 0000027E0 / DATA AI02CS(22) / .0000000000 0000003E0 / C C---------------------------------------------------------------------- C DATA NTI0 / 0 / C C***FIRST EXECUTABLE STATEMENT WI0 C IF (M .LE. 0) GO TO 910 C IF (NTI0.EQ.0) THEN EPMACH = R1MACH(3) EPS = 0.10E0*EPMACH NTI0 = IWCS(BI0CS, LBI0, EPS) NTAI0 = IWCS(AI0CS, LAI0, EPS) NTAI02 = IWCS(AI02CS, LAI02, EPS) XSML = SQRT(4.0E0*EPMACH) XMAX = LOG(R1MACH(2)) ENDIF C DO 10 I=1,M Y(I) = ABS(X(I)) 10 CONTINUE C CALL WNGT(M,Y,XMAX,KEY) IF (KEY .NE. 0) GO TO 920 C C ---------------- C CASE Y .LE. XSML C ---------------- C DO 15 I=1,M F(I) = 1.0E0 15 CONTINUE C C -------------------------- C CASE XSML .LT. Y .LE. 3.0 C -------------------------- C CALL WGTLE(M,Y,XSML,3.0E0,N,INDX) IF (N .GT. 0) THEN CALL WGTHR(N,Y,INDX,YCMP) DO 20 J=1,N TCMP(J) = YCMP(J)**2/4.50E0 - 1.0E0 20 CONTINUE CALL WCS(N,TCMP,BI0CS,NTI0,ZCMP,B0,B1,B2) DO 30 J=1,N ZCMP(J) = 2.750E0 + ZCMP(J) 30 CONTINUE CALL WSCTR(N,ZCMP,INDX,F) ENDIF C C ------------------------- C CASE 3.0 .LT. Y .LE. 8.0 C ------------------------- C CALL WGTLE(M,Y,3.0E0,8.0E0,N,INDX) IF (N .GT. 0) THEN CALL WGTHR(N,Y,INDX,YCMP) DO 50 J=1,N TCMP(J) = (48.0E0/YCMP(J) - 11.0E0)/5.0E0 50 CONTINUE CALL WCS(N,TCMP,AI0CS,NTAI0,ZCMP,B0,B1,B2) DO 60 J=1,N ZCMP(J) = EXP(YCMP(J))*(0.3750E0+ZCMP(J))/SQRT(YCMP(J)) 60 CONTINUE CALL WSCTR(N,ZCMP,INDX,F) ENDIF C C ---------------- C CASE Y .GT. 8.0 C ---------------- C CALL WGT(M,Y,8.0E0,N,INDX) IF (N .GT. 0) THEN CALL WGTHR(N,Y,INDX,YCMP) DO 80 J=1,N TCMP(J) = 16.0E0/YCMP(J) - 1.0E0 80 CONTINUE CALL WCS(N,TCMP,AI02CS,NTAI02,ZCMP,B0,B1,B2) DO 90 J=1,N ZCMP(J) = EXP(YCMP(J))*(0.3750E0+ZCMP(J))/SQRT(YCMP(J)) 90 CONTINUE CALL WSCTR(N,ZCMP,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 ... ABS(X) SO LARGE I0 OVERFLOWS C 920 CONTINUE INFO = 2 INDX(1) = KEY GO TO 999 C 999 CONTINUE RETURN END