SUBROUTINE VI1 (M, X, F, WORK, IWORK, INFO) C***BEGIN PROLOGUE VI1 C***PURPOSE Computes the hyperbolic Bessel function of the first kind C of order one (I1) for a vector of real arguments C***LIBRARY VFNLIB C***CATEGORY C10B1 C***TYPE SINGLE PRECISION (VI1-S, DVI1-D) C***KEYWORDS BESSEL FUNCTION,FIRST 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 VI1 computes the modified (hyperbolic) Bessel function of the C first kind of order one (I1) 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 abs(X(i)) so small I1 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 abs(X(i)) so big I1 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 BESI1 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 WI1 C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C 930203 Minor modifications to prologue. C***END PROLOGUE VI1 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 VI1 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 ... WI1 DOES ALL THE WORK C CALL WI1(M,X,F,WORK(IWY),WORK(IWTC),WORK(IWYC),WORK(IWZC), + IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO) C RETURN END SUBROUTINE WI1(M, X, F, Y, TCMP, YCMP, ZCMP, INDX, B0, B1, B2, + INFO) C***BEGIN PROLOGUE WI1 C***SUBSIDIARY C***PURPOSE Computes the hyperbolic Bessel function of the first kind C of order one (I1) 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 WI1 C C ---------- C PARAMETERS C ---------- C INTEGER LAI1, LAI12, LBI1 PARAMETER ( LAI1=21, LAI12=22, LBI1=11 ) 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 I, IWCS, J, KEY, N, NTI1, NTAI1, NTAI12, NTOT REAL AI1CS, AI12CS, BI1CS, EPMACH, EPS, R1MACH, + XMAX, XMIN, XSML C DIMENSION AI1CS(LAI1), AI12CS(LAI12), BI1CS(LBI1) C SAVE AI1CS, AI12CS, BI1CS, N, NTAI1, NTAI12, NTI1, + XMAX, XMIN, XSML 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 AI1 on the interval 1.25000D-01 to 3.33333D-01 C with weighted error 6.98E-17 C log weighted error 16.16 C significant figures required 14.53 C decimal places required 16.82 C DATA AI1 CS( 1) / -.0284674418 1881479E0 / DATA AI1 CS( 2) / -.0192295323 1443221E0 / DATA AI1 CS( 3) / -.0006115185 8579437E0 / DATA AI1 CS( 4) / -.0000206997 1253350E0 / DATA AI1 CS( 5) / .0000085856 1914581E0 / DATA AI1 CS( 6) / .0000010494 9824671E0 / DATA AI1 CS( 7) / -.0000002918 3389184E0 / DATA AI1 CS( 8) / -.0000000155 9378146E0 / DATA AI1 CS( 9) / .0000000131 8012367E0 / DATA AI1 CS(10) / -.0000000014 4842341E0 / DATA AI1 CS(11) / -.0000000002 9085122E0 / DATA AI1 CS(12) / .0000000001 2663889E0 / DATA AI1 CS(13) / -.0000000000 1664947E0 / DATA AI1 CS(14) / -.0000000000 0166665E0 / DATA AI1 CS(15) / .0000000000 0124260E0 / DATA AI1 CS(16) / -.0000000000 0027315E0 / DATA AI1 CS(17) / .0000000000 0002023E0 / DATA AI1 CS(18) / .0000000000 0000730E0 / DATA AI1 CS(19) / -.0000000000 0000333E0 / DATA AI1 CS(20) / .0000000000 0000071E0 / DATA AI1 CS(21) / -.0000000000 0000006E0 / C C---------------------------------------------------------------------- C C Series for AI12 on the interval 0. to 1.25000D-01 C with weighted error 3.55E-17 C log weighted error 16.45 C significant figures required 14.69 C decimal places required 17.12 C DATA AI12CS( 1) / .0285762350 1828014E0 / DATA AI12CS( 2) / -.0097610974 9136147E0 / DATA AI12CS( 3) / -.0001105889 3876263E0 / DATA AI12CS( 4) / -.0000038825 6480887E0 / DATA AI12CS( 5) / -.0000002512 2362377E0 / DATA AI12CS( 6) / -.0000000263 1468847E0 / DATA AI12CS( 7) / -.0000000038 3538039E0 / DATA AI12CS( 8) / -.0000000005 5897433E0 / DATA AI12CS( 9) / -.0000000000 1897495E0 / DATA AI12CS(10) / .0000000000 3252602E0 / DATA AI12CS(11) / .0000000000 1412580E0 / DATA AI12CS(12) / .0000000000 0203564E0 / DATA AI12CS(13) / -.0000000000 0071985E0 / DATA AI12CS(14) / -.0000000000 0040836E0 / DATA AI12CS(15) / -.0000000000 0002101E0 / DATA AI12CS(16) / .0000000000 0004273E0 / DATA AI12CS(17) / .0000000000 0001041E0 / DATA AI12CS(18) / -.0000000000 0000382E0 / DATA AI12CS(19) / -.0000000000 0000186E0 / DATA AI12CS(20) / .0000000000 0000033E0 / DATA AI12CS(21) / .0000000000 0000028E0 / DATA AI12CS(22) / -.0000000000 0000003E0 / C C---------------------------------------------------------------------- C DATA NTI1 / 0 / C C***FIRST EXECUTABLE STATEMENT WI1 C IF (M .LE. 0) GO TO 910 C IF (NTI1.EQ.0) THEN EPMACH = R1MACH(3) EPS = 0.10E0*EPMACH NTI1 = IWCS (BI1CS, LBI1, EPS) NTAI1 = IWCS (AI1CS, LAI1, EPS) NTAI12 = IWCS (AI12CS, LAI12, EPS) XMIN = 2.0E0*R1MACH(1) XSML = SQRT(8.0E0*EPMACH) XMAX = LOG(R1MACH(2)) ENDIF C NTOT = 0 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=0 OR Y TOO SMALL C --------------------------- C C NOTE -- I0 UNDERFLOWS FOR X .LE. XMIN C DO 15 I=1,M F(I) = 0.0E0 15 CONTINUE C C ---------------------------- C CASE XMIN .LT. Y .LE. XSML C ---------------------------- C CALL WGTLE(M,Y,XMIN,XSML,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL WGTHR(N,Y,INDX,YCMP) DO 20 J=1,N ZCMP(J) = 0.50E0*YCMP(J) 20 CONTINUE CALL WSCTR(N,ZCMP,INDX,F) ENDIF 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 NTOT = NTOT + N CALL WGTHR(N,Y,INDX,YCMP) DO 30 J=1,N TCMP(J) = YCMP(J)**2/4.50E0 - 1.0E0 30 CONTINUE CALL WCS(N,TCMP,BI1CS,NTI1,ZCMP,B0,B1,B2) DO 40 J=1,N ZCMP(J) = YCMP(J)*(0.8750E0 + ZCMP(J)) 40 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 NTOT = NTOT + N CALL WGTHR(N,Y,INDX,YCMP) DO 70 J=1,N TCMP(J) = (48.0E0/YCMP(J) - 11.0E0)/5.0E0 70 CONTINUE CALL WCS(N,TCMP,AI1CS,NTAI1,ZCMP,B0,B1,B2) DO 80 J=1,N ZCMP(J) = EXP(YCMP(J))*(0.3750E0 + ZCMP(J))/SQRT(YCMP(J)) 80 CONTINUE CALL WSCTR(N,ZCMP,INDX,F) ENDIF C C ---------------- C CASE 8.0 .LT. Y C ---------------- C CALL WGT(M,Y,8.0E0,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL WGTHR(N,Y,INDX,YCMP) DO 100 J=1,N TCMP(J) = 16.0E0/YCMP(J) - 1.0E0 100 CONTINUE CALL WCS(N,TCMP,AI12CS,NTAI12,ZCMP,B0,B1,B2) DO 110 J=1,N ZCMP(J) = EXP(YCMP(J))*(0.3750E0 + ZCMP(J))/SQRT(YCMP(J)) 110 CONTINUE CALL WSCTR(N,ZCMP,INDX,F) ENDIF C C ------------------------------------------- C REVERSE SIGN FOR NEGATIVE X (RESULT IS ODD) C ------------------------------------------- C DO 200 I=1,M IF (X(I) .LT. 0.0E0) F(I) = -F(I) 200 CONTINUE 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 ... ABS(X) SO LARGE I1 OVERFLOWS C 920 CONTINUE INFO = 2 INDX(1) = KEY GO TO 999 C 999 CONTINUE RETURN END