      SUBROUTINE VJ1 (M, X, F, WORK, IWORK, INFO)
C***BEGIN PROLOGUE  VJ1
C***PURPOSE  Computes the Bessel function of the first kind
C            of order one (J1) for a vector of real arguments
C***LIBRARY   VFNLIB
C***CATEGORY  C10A1
C***TYPE      SINGLE PRECISION (VJ1-S, DVJ1-D)
C***KEYWORDS  BESSEL FUNCTION, FIRST KIND, SPECIAL FUNCTION,
C             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   VJ1 computes the  Bessel function of the first kind of order
C   one (J1) 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
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 J1 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 that no precision
C                     possible in computing J1.  The index of the
C                     first offending argument is returned in IWORK(1).
C
C *********************************************************************
C   This routine is a modification of the function BESJ1 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  WJ1
C***REVISION HISTORY  (YYMMDD)
C   910226  DATE WRITTEN
C   930203  Minor modifications to prologue.
C***END PROLOGUE  VJ1
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  VJ1
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     ... WJ1 DOES ALL THE WORK
C
      CALL WJ1(M,X,F,WORK(IWY),WORK(IWTC),WORK(IWYC),WORK(IWZC),
     +         IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO)
C
      RETURN
      END
      SUBROUTINE WJ1(M, X, F, Y, TCMP, YCMP, ZCMP, INDX, B0, B1, B2,
     +   INFO)
C***BEGIN PROLOGUE  WJ1
C***SUBSIDIARY
C***PURPOSE  Computes the Bessel function of the first kind
C            of order one (J1) 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,
C                    WLE, WSCTR, WCS
C***REVISION HISTORY  (YYMMDD)
C   910226  DATE WRITTEN
C***END PROLOGUE  WJ1
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 LBJ1, LBM1, LBM12, LBTH1, LBT12
      PARAMETER (LBJ1=12, LBM1=21, LBM12=2, LBTH1=24, LBT12=2)
C
      INTEGER I, IWCS, J, JH, K, KEY, N, NA, NB, NTJ1, NTM1, NTM12,
     +        NTTH1, NTT12, NTOT
      REAL    BJ1CS, BM1CS, BM12CS, BTH1CS, BT12CS, C1, C2,
     +        EPMACH, EPS, PI4, R1MACH, TPI4, XMAX, XMIN, XSML
C
      DIMENSION BJ1CS(LBJ1), BM1CS(LBM1), BM12CS(LBM12), BTH1CS(LBTH1),
     +          BT12CS(LBT12)
C
      SAVE BJ1CS, BM1CS, BM12CS, BTH1CS, BT12CS, NTJ1, NTM1, NTM12,
     +      NTTH1, NTT12, PI4, TPI4, XMIN, XSML, XMAX
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
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
      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
      DATA NTJ1 / 0 /
C
C***FIRST EXECUTABLE STATEMENT  WJ1
C
      IF (M .LE. 0) GO TO 910
C
      IF (NTJ1.EQ.0) THEN
         EPMACH = R1MACH(3)
         EPS = 0.10E0*EPMACH
         NTJ1  = IWCS(BJ1CS, LBJ1, EPS)
         NTM1  = IWCS(BM1CS, LBM1, EPS)
         NTM12 = IWCS(BM12CS, LBM12, EPS)
         NTTH1 = IWCS(BTH1CS, LBTH1, EPS)
         NTT12 = IWCS(BT12CS, LBT12, EPS)
         XMIN = 2.0E0*R1MACH(1)
         XSML = SQRT(8.0E0*EPMACH)
         XMAX = 1.0E0/R1MACH(4)
         PI4 = ATAN(1.0E0)
         TPI4 = 3.0E0*PI4
      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 --- J1 UNDERFLOWS FOR  0 .LT. Y .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. 4.0
C  ---------------------------
C
      CALL WGTLE(M,Y,XSML,4.0E0,N,INDX)
      IF (N .GT. 0) THEN
         NTOT = NTOT + N
         CALL WGTHR(N,Y,INDX,YCMP)
         DO 30 J=1,N
            TCMP(J) = 0.1250E0*YCMP(J)**2 - 1.0E0
  30     CONTINUE
         CALL WCS(N,TCMP,BJ1CS,NTJ1,ZCMP,B0,B1,B2)
         DO 40 K=1,N
            ZCMP(K) = YCMP(K)*(0.250E0 + ZCMP(K))
  40     CONTINUE
         CALL WSCTR(N,ZCMP,INDX,F)
      ENDIF
C
C  ----------------
C  CASE  Y .GT. 4.0
C  ----------------
C
      CALL WGT(M,Y,4.0E0,N,INDX)
      IF (N .GT. 0) THEN
         NTOT = NTOT + N
         CALL WGTHR(N,Y,INDX,YCMP)
         DO 50 J=1,N
            TCMP(J) = 32.0E0/YCMP(J)**2 - 1.0E0
  50     CONTINUE
         CALL WCS(N,TCMP,BM1CS,NTM1,Y,B0,B1,B2)
         CALL WCS(N,TCMP,BTH1CS,NTTH1,ZCMP,B0,B1,B2)
         DO 60 J=1,N
            Y(J) = (0.750E0 + Y(J)) / SQRT(YCMP(J))
            ZCMP(J) = ( YCMP(J) - TPI4 ) + ZCMP(J)/YCMP(J)
            ZCMP(J) = Y(J) * COS(ZCMP(J))
  60     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     ... NO PRECISION BECAUSE X BIG
C
  920 CONTINUE
      INFO = 2
      INDX(1) = KEY
      GO TO 999
C
  999 CONTINUE
      RETURN
      END