DOUBLE PRECISION FUNCTION DSI (X)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 1998-10-29 DSI Krogh  Moved external statement up for mangle.
C>> 1995-11-03 DSI Krogh  Removed blanks in numbers for C conversion.
C>> 1994-11-11 DSI Krogh  Declared all vars.
c>> 1994-10-20 DSI Krogh  Changes to use M77CON
c>> 1989-03-14 DSI Original W. V. Snyder at JPL
C
C     COMPUTE THE SINE INTEGRAL OF X =
C     INTEGRAL FROM 0 TO X OF (SIN(T)/T DT).
C
C     FOR ABS(X)<16, USE A CHEBYSHEV SERIES WITH ARGUMENT 2*Z*Z-1 WHERE
C     Z=X/16 TO EVALUATE SI(X)/Z, THEN MULTIPLY THE RESULT BY Z.  THIS
C     AVOIDS STORING ZERO COEFFICIENTS FOR EVEN ORDERS, AND PRESERVES
C     ACCURACY FOR SMALL Z.
C
C     FOR 16.LE.ABS(X).LE.100, USE CHEBYCHEV SERIES WITH ARGUMENT
C     2*Z*Z-1, WHERE Z=16/X ARE USED TO COMPUTE F(X)/X AND G(X)/(X*X).
C     THEN SI(X)=0.5*PI*SIGN(X)-F(X)/X*COS(X)-G(X)/(X*X)*SIN(X).
C
C     WHEN X.GT.100, USE ASYMPTOTIC APPROXIMATIONS FOR F(X)/X AND
C     G(X)/(X*X) AND COMPUTE SI(X) AS ABOVE.
C
C     THIS ALGORITHM YIELDS AT MOST 15 DIGITS OF PRECISION.
C
C--D replaces "?": ?SI, ?CPVAL
C
      INTEGER N
      DOUBLE PRECISION X
      EXTERNAL DCPVAL
      DOUBLE PRECISION PI2,Z,ZW,FZ,GZ,DCPVAL
      DOUBLE PRECISION FT,GT
      DOUBLE PRECISION S(23),F(13),G(13)
      DATA PI2/1.57079632679489662D0/
      DATA S/
     * + 0.5D0,                 + 0.5D0,
     1 + 4.052926477680623D0, - 4.063980844911986D0,
     2 + 2.778756381742663D0, - 1.926565091150656D0,
     3 + 1.389308771171888D0, - 0.968322236987086D0,
     4 + 0.530148847916522D0, - 0.211263780976555D0,
     5 + 0.062033679432003D0, - 0.013867445589417D0,
     6 + 0.002436221404749D0, - 0.000345469155569D0,
     7 + 0.000040420271419D0, - 0.000003972908746D0,
     8 + 0.000000332988589D0, - 0.000000024100076D0,
     9 + 0.000000001522370D0, - 0.000000000084710D0,
     A + 0.000000000004185D0, - 0.000000000000185D0,
     B + 0.000000000000007D0/
      DATA F/
     * + 0.5D0,                 + 0.5D0,
     1 + 0.062263729028927D0, - 0.000233756041393D0,
     2 + 0.000002453755677D0, - 0.000000058670317D0,
     3 + 0.000000002356196D0, - 0.000000000136096D0,
     4 + 0.000000000010308D0, - 0.000000000000964D0,
     5 + 0.000000000000107D0, - 0.000000000000014D0,
     6 + 0.000000000000002D0/
      DATA G/
     * + 0.5D0,                 + 0.5D0,
     1 + 0.003862856096703D0, - 0.000042644182622D0,
     2 + 0.000000724995950D0, - 0.000000023468225D0,
     3 + 0.000000001169202D0, - 0.000000000079604D0,
     4 + 0.000000000006875D0, - 0.000000000000717D0,
     5 + 0.000000000000087D0, - 0.000000000000012D0,
     6 + 0.000000000000002D0/
C
      IF (ABS(X).LT.16.0) THEN
         Z = X/16.0
         ZW = Z*Z
         Z = Z*DCPVAL(S,20,ZW)
      ELSE
        IF (ABS(X).LE.100.0) THEN
C           16.LE.ABS(X).LE.100
            Z = 16.0/X
            ZW = Z*Z
            FZ = Z*DCPVAL(F,10,ZW)
            GZ = ZW*DCPVAL(G,10,ZW)
        ELSE
C           ABS(X).GT.100
            FZ = 1.0/X
            FT = FZ
            GZ = FZ/X
            GT = GZ
            Z = GZ
            DO 25 N = 2, 16, 2
               FT = -DBLE(N*(N-1))*Z*FT
               GT = -DBLE(N*(N+1))*Z*GT
               FZ = FZ + FT
               GZ = GZ + GT
25          CONTINUE
         END IF
         Z = SIGN(PI2,X) - FZ*COS(X) - GZ*SIN(X)
      END IF
C
      DSI = Z
      RETURN
      END

.