DECLARE SUB IntegralSimpson (kind%, a#, b#, n%, res#)
DECLARE FUNCTION F# (kind%, t#)
DECLARE FUNCTION sinintegral# (x#)
DECLARE FUNCTION cosintegral# (x#)
'**************************************************************
'* Program to demonstrate the use of functions sinintegral(x) *
'* or cosintegral(x) using the Simpson's method               *
'* ---------------------------------------------------------- *
'*                                                            *
'* SAMPLE RUNS:                                               *
'*                                                            *
'*         Kind (1 or 2)= 1                                   *
'*                                                            *
'* Test of Function sinintegral(x) by Simpson's method        *
'*                                                            *
'* Input value for x variable:                                *
'*                                                            *
'*         X = 2                                              *
'*                                                            *
'* Result:  1.605412976802747                                 *
'*                                                            *
'*         X = 5                                              *
'*                                                            *
'* Result:  1.549931244944637                                 *
'*                                                            *
'*         X = 80                                             *
'*                                                            *
'* Result:  1.572330886912482                                 *
'*                                                            *
'*                                                            *
'*         Kind (1 or 2) = 2                                  *
'*                                                            *
'* Test of Function cosintegral(x) by Simpson's method        *
'*                                                            *
'* Input value for x variable:                                *
'*                                                            *
'*         X = 2                                              *
'*                                                            *
'* Result:  .4229808287748087                                 *
'*                                                            *
'*         X = 5                                              *
'*                                                            *
'* Result: -.1900297496567194                                 *
'*                                                            *
'*         X = 80                                             *
'*                                                            *
'* Result: -1.240250115512912D-02                             *
'*                                                            *
'*                              Basic Version By J-P Moreau.  *
'*                                   (www.jpmoreau.fr)        *
'**************************************************************
'Program Test_Sinintegral
DEFDBL A-H, O-Z
DEFINT I-N

  CLS
  PRINT
  INPUT "         Kind (1 or 2) = ", kind
  PRINT
  IF kind = 1 THEN
    PRINT " Test of Function sinintegral(x) by Simpson's method"
  ELSE
    PRINT " Test of Function cosintegral(x) by Simpson's method"
  END IF
  PRINT
  PRINT " Input value for x variable:"
  PRINT
  INPUT "         X = ", x
  PRINT
  IF kind = 2 AND x <= 0# THEN
    PRINT " WARNING! X must be positive."
  ELSE
    IF kind = 1 THEN
      PRINT " Result: "; sinintegral(x)
    ELSE
      PRINT " Result: "; cosintegral(x)
    END IF
  END IF
  PRINT

END 'of main program

'*******************************************************
'* Calcualte Ci(x)= gamma + Ln(x) + I(x) with I(x) =   *
'* Integral of (cos(u)-1)/u from u=0 to u=x (x > 0).   *
'* Gamma = 0.57721566... (Euler's constant).           *
'*******************************************************
FUNCTION cosintegral (x)
  'x0     begin x value
  'x      end x value
  'nstep  number of integration steps
  'res    result of integral

  gamma = .5772156649015329#    'Euler's constant

  x0 = 0#
  nstep = INT(200 * ABS(x)) 'this should ensure about 14 exact digits

  IntegralSimpson 2, x0, x, nstep, res  'kind=2
  cosintegral = gamma + LOG(x) + res

END FUNCTION

'Given function to integrate
FUNCTION F (kind, t)
  IF kind = 1 THEN
    IF ABS(t) < .0000000001# THEN
      F = 1#
    ELSE
      F = SIN(t) / t              'for sinintegral
    END IF
  ELSE
    IF ABS(t) < .0000000001# THEN
      F = 0#
    ELSE
      F = (COS(t) - 1#) / t       'for cosintegral
    END IF
  END IF
END FUNCTION

'*******************************************************
'* Integral of a function F(X) by Simpson's method     *
'* --------------------------------------------------- *
'* INPUTS:                                             *
'*          kind   =1 for sinintegral                  *
'*                 =2 for cosintegral                  *
'*          a      begin value of x variable           *
'*          b      end value of x variable             *
'*          n      number of integration steps         *
'*                                                     *
'* OUTPUT:                                             *
'*          res    the integral of FUNC(X) from a to b *
'*                                                     *
'*******************************************************
SUB IntegralSimpson (kind, a, b, n, res)
  xstep = (b - a) / 2 / n
  r = F(kind, a)
  res = (r + F(kind, b)) / 2#
  FOR i = 1 TO 2 * n - 1
    r = F(kind, a + i * xstep)
    IF (i MOD 2) <> 0 THEN
      res = res + r + r
    ELSE
      res = res + r
    END IF
  NEXT i
  res = res * xstep * 2# / 3#
END SUB

'*******************************************************
'*     Si(x) = Integral of sin(u)/u from u=0 to u=x    *
'*******************************************************
FUNCTION sinintegral (x)
  'x0     begin x value
  'x      end x value
  'nstep  number of integration steps
  'res    result of integral

  x0 = 0#
  nstep = INT(200 * ABS(x)) 'this should ensure about 14 exact digits

  IntegralSimpson 1, x0, x, nstep, res  'kind=1
  sinintegral = res

END FUNCTION