      FUNCTION SI(X)
C
C     WRITTEN BY D.E. AMOS AND S.L. DANIEL, NOVEMBER,1975.
C     MODIFIED BY A.H. MORRIS
C
C     REFERENCES
C         SAND76-006C
C         THE SPECIAL FUNCTIONS AND THEIR APPROXIMATIONS, VOL. II, BY
C         Y.L. LUKE. ACADEMIC PRESS, NEW YORK, 1969.
C
C     ABSTRACT
C         SI COMPUTES THE INTEGRAL OF SIN(T)/T ON (0,X) BY MEANS
C         OF CHEBYSHEV EXPANSIONS ON (0,5) AND (5,INFINITY).
C
C     DESCRIPTION OF ARGUMENTS
C
C         INPUT
C           X      - LIMIT OF INTEGRATION, UNRESTRICTED
C
C         OUTPUT
C           SI     - VALUE OF THE INTEGRAL
C
C
      DIMENSION BB(16),CC(46)
C     -------------------
      DATA N1,N2,M1,M2/16,46,14,21/
      DATA PIO2/1.5707963267949/
C     -------------------
      DATA BB(1) / 6.84101190850653E-01/, BB(2) /-3.74538448460062E-01/,
     1     BB(3) /-2.82656062651542E-02/, BB(4) / 3.06078454012071E-02/,
     2     BB(5) /-8.99242948380352E-04/, BB(6) /-1.09884251456048E-03/,
     3     BB(7) / 5.81151604367358E-05/, BB(8) / 2.28802638122969E-05/,
     4     BB(9) /-1.35078982929539E-06/, BB(10)/-3.13213946132892E-07/,
     5     BB(11)/ 1.86619586786257E-08/, BB(12)/ 3.03991719607226E-09/,
     6     BB(13)/-1.76437788946489E-10/, BB(14)/-2.20236421792690E-11/,
     7     BB(15)/ 1.22710107703240E-12/, BB(16)/ 1.23680681116783E-13/
C     -------------------
      DATA CC(1) / 9.76155271128712E-01/, CC(2) / 8.96845854916423E-02/,
     1     CC(3) /-3.04656658030696E-02/, CC(4) / 8.50892472922945E-02/,
     2     CC(5) /-5.78073683148386E-03/, CC(6) /-5.07182677775691E-03/,
     3     CC(7) / 8.38643256650893E-04/, CC(8) /-3.34223415981738E-04/,
     4     CC(9) /-2.15746207281216E-05/, CC(10)/ 1.28560650086065E-04/,
     5     CC(11)/-1.56456413510232E-05/, CC(12)/-1.52025513597262E-05/,
     6     CC(13)/ 4.04001013843204E-06/, CC(14)/-5.95896122752160E-07/,
     7     CC(15)/-4.34985305974340E-07/, CC(16)/ 7.13472533530840E-07/,
     8     CC(17)/-5.34302186061100E-08/, CC(18)/-1.76003581156610E-07/,
     9     CC(19)/ 3.85028855125900E-08/, CC(20)/ 1.92576544441700E-08/,
     1     CC(21)/-1.00735358217200E-08/, CC(22)/ 3.36359194377000E-09/,
     2     CC(23)/ 1.28049619406000E-09/, CC(24)/-2.42546870827000E-09/,
     3     CC(25)/ 1.86917288950000E-10/, CC(26)/ 7.13431298340000E-10/,
     4     CC(27)/-1.70673483710000E-10/, CC(28)/-1.14604070350000E-10/,
     5     CC(29)/ 5.88004411500000E-11/, CC(30)/-6.78417843000000E-12/,
     6     CC(31)/-1.21572380900000E-11/, CC(32)/ 1.26561248700000E-11/,
     7     CC(33)/ 4.74814180000000E-13/, CC(34)/-5.32309477000000E-12/,
     8     CC(35)/ 9.05903810000000E-13/, CC(36)/ 1.40046450000000E-12/,
     9     CC(37)/-5.00968320000000E-13/, CC(38)/-1.80458040000000E-13/
      DATA CC(39)/ 1.66162910000000E-13/, CC(40)/-5.02616400000000E-14/,
     1     CC(41)/-3.48453600000000E-14/, CC(42)/ 4.60056600000000E-14/,
     2     CC(43)/ 5.74000000000000E-16/, CC(44)/-1.95310700000000E-14/,
     3     CC(45)/ 3.68837000000000E-15/, CC(46)/ 5.62862000000000E-15/
C     -------------------
C     ****** AMAX IS A MACHINE DEPENDENT CONSTANT. IT IS ASSUMED THAT
C            SIN(X) AND COS(X) ARE DEFINED FOR  ABS(X) .LE. AMAX, AND
C            THAT PIO2 - (1 + 1/X)/X = PIO2 FOR X .GT. AMAX.
C
      AMAX = 1.0/SPMPAR(1)
C
C
      AX=ABS(X)
      IF (AX.GT.5.0) GO TO 20
      J=N1
      BX=0.40*AX-1.0
      TX=BX+BX
      B1=BB(J)
      B2=0.
      DO 10 I=1,M1
      J=J-1
      TEMP=B1
      B1=TX*B1-B2+BB(J)
   10 B2=TEMP
      SI=(BX*B1-B2+BB(1))*X
      RETURN
C
   20 IF (AX.GT.AMAX) GO TO 50
      BX=10./AX-1.
      TX=BX+BX
      J=N2
      B1=CC(J)
      B2=0.0
      DO 30 I=1,M2
      J=J-2
      TEMP=B1
      B1=TX*B1-B2+CC(J)
   30 B2=TEMP
      AIC=BX*B1-B2+CC(2)
C
      J=N2-1
      B1=CC(J)
      B2=0.0
      DO 40 I=1,M2
      J=J-2
      TEMP=B1
      B1=TX*B1-B2+CC(J)
   40 B2=TEMP
      RC=BX*B1-B2+CC(1)
C
      SI=(RC*COS(AX)+AIC*SIN(AX))/AX
      SI=PIO2-SI
      IF (X.LT.0.0) SI=-SI
      RETURN
C
   50 SI=SIGN(PIO2,X)
      RETURN
      END
