; $Id: //depot/Release/ENVI52_IDL84/idl/idldir/lib/int_tabulated.pro#1 $
;
; Copyright (c) 1995-2014, Exelis Visual Information Solutions, Inc. All
;       rights reserved. Unauthorized reproduction is prohibited.
;+
; NAME:
;       INT_TABULATED
;
; PURPOSE:
;       This function integrates a tabulated set of data { x(i) , f(i) },
;       on the closed interval [min(X) , max(X)].
;
; CATEGORY:
;       Numerical Analysis.
;
; CALLING SEQUENCE:
;       Result = INT_TABULATED(X, F)
;
; INPUTS:
;       X:  The tabulated X-value data. This data may be irregularly
;           gridded and in random order. If the data is randomly ordered
;	    you must set the SORT keyword to a nonzero value.
;           Duplicate x values will result in a warning message.
;       F:  The tabulated F-value data. Upon input to the function
;           X(i) and F(i) must have corresponding indices for all
;	    values of i. If X is reordered, F is also reordered.
;
;       X and F must be of floating point or double precision type.
;
; KEYWORD PARAMETERS:
;       SORT:   A zero or non-zero scalar value.
;               SORT = 0 (the default) The tabulated x-value data is
;                        already in ascending order.
;               SORT = 1 The tabulated x-value data is in random order
;                        and requires sorting into ascending order. Both
;			 input parameters X and F are returned sorted.
;       DOUBLE: If set to a non-zero value, computations are done in
;               double precision arithmetic.
;
; OUTPUTS:
;       This fuction returns the integral of F computed from the tabulated
;	data in the closed interval [min(X) , max(X)].
;
; RESTRICTIONS:
;       Data that is highly oscillatory requires a sufficient number
;       of samples for an accurate integral approximation.
;
; PROCEDURE:
;       INT_TABULATED.PRO constructs a regularly gridded x-axis with a
;	number of segments as an integer multiple of four. Segments
;	are processed in groups of four using a 5-point Newton-Cotes
;	integration formula.
;       For 'sufficiently sampled' data, this algorithm is highly accurate.
;
; EXAMPLES:
;       Example 1:
;       Define 11 x-values on the closed interval [0.0 , 0.8].
;         x = [0.0, .12, .22, .32, .36, .40, .44, .54, .64, .70, .80]
;
;       Define 11 f-values corresponding to x(i).
;         f = [0.200000, 1.30973, 1.30524, 1.74339, 2.07490, 2.45600, $
;              2.84299,  3.50730, 3.18194, 2.36302, 0.231964]
;
;       Compute the integral.
;         result = INT_TABULATED(x, f)
;
;       In this example, the f-values are generated from a known function,
;       (f = .2 + 25*x - 200*x^2 + 675*x^3 - 900*x^4 + 400*x^5)
;
;       The Multiple Application Trapazoid Method yields;  result = 1.5648
;       The Multiple Application Simpson's Method yields;  result = 1.6036
;				INT_TABULATED.PRO yields;  result = 1.6232
;         The Exact Solution (4 decimal accuracy) yields;  result = 1.6405
;
;	Example 2: 
;       Create 30 random points in the closed interval [-2 , 1].
;         x = randomu(seed, 30) * 3.0 - 2.0
;
;       Explicitly define the interval's endpoints.
;         x(0) = -2.0  &  x(29) = 1.0
;
;       Generate f(i) corresponding to x(i) from a given function.
;         f = sin(2*x) * exp(cos(2*x))
;
;       Call INT_TABULATED with the SORT keyword.
;         result = INT_TABULATED(x, f, /sort)
;
;       In this example, the f-values are generated from the function,
;       f = sin(2*x) * exp(cos(2*x))
;
;       The result of this example will vary because the x(i) are random.
;       Executing this example three times gave the following results:
;		               INT_TABULATED.PRO yields;  result = -0.0702
;		               INT_TABULATED.PRO yields;  result = -0.0731
;		               INT_TABULATED.PRO yields;  result = -0.0698
;        The Exact Solution (4 decimal accuracy) yields;  result = -0.0697
;
; MODIFICATION HISTORY:
;           Written by:  GGS, RSI, September 1993
;           Modified:    GGS, RSI, November  1993
;                        Use Numerical Recipes cubic spline interpolation 
;                        function NR_SPLINE/NR_SPLINT. Execution time is 
;                        greatly reduced. Added DOUBLE keyword. The 'sigma' 
;                        keyword is no longer supported.
;           Modified:    GGS, RSI, April  1995
;                        Changed cubic spline calls from NR_SPLINE/NR_SPLINT
;                        to SPL_INIT/SPL_INTERP. Improved double-precision
;                        accuracy.
;           Modified:    GGS, RSI, April 1996
;                        Replaced WHILE loop with vector operations. 
;                        Check for duplicate points in x vector.  
;                        Modified keyword checking and use of double precision.
;-

FUNCTION Int_Tabulated, X, F, Double = Double, Sort = Sort

  ;Return to caller if an error occurs.
  ON_ERROR, 2 

  TypeX = SIZE(X)
  TypeF = SIZE(F)

  ;Check F data type.
  if TypeF[TypeF[0]+1] ne 4 and TypeF[TypeF[0]+1] ne 5 then $
    MESSAGE, "F values must be float or double."

  ;Check length.
  if TypeX[TypeX[0]+2] ne TypeF[TypeF[0]+2] then $
    MESSAGE, "X and F arrays must have the same number of elements."

  ;Check duplicate values.
  if TypeX[TypeX[0]+2] ne N_ELEMENTS(UNIQ(X[SORT(X)])) then $
    MESSAGE, "X array contains duplicate points."

  ;If the DOUBLE keyword is not set then the internal precision and
  ;result are identical to the type of input.
  if N_ELEMENTS(Double) eq 0 then $
    Double = (TypeX[TypeX[0]+1] eq 5 or TypeF[TypeF[0]+1] eq 5) 

  Xsegments = TypeX[TypeX[0]+2] - 1L

  ;Sort vectors into ascending order.
  if KEYWORD_SET(Sort) ne 0 then begin
    ii = SORT(x)
    X = X[ii]
    F = F[ii]
  endif

  while (Xsegments MOD 4L) ne 0L do $
    Xsegments = Xsegments + 1L

  Xmin = MIN(X)
  Xmax = MAX(X)

  ;Uniform step size.
    h = (Xmax+0.0 - Xmin) / Xsegments
  ;Compute the interpolates at Xgrid.
    ;x values of interpolates >> Xgrid = h * FINDGEN(Xsegments + 1L) + Xmin
    z = SPL_INTERP(X, F, SPL_INIT(X, F, Double = Double), $
                   h * FINDGEN(Xsegments + 1L) + Xmin, Double = Double)
  ;Compute the integral using the 5-point Newton-Cotes formula.
    ii = (LINDGEN((N_ELEMENTS(z) - 1L)/4L)+1) * 4
    if Double eq 0 then $
      RETURN, FLOAT(TOTAL(2.0 * h * (7.0 * (z[ii-4] + z[ii]) + $
                    32.0 * (z[ii-3] + z[ii-1]) + 12.0 * z[ii-2]) / 45.0)) $
    else $
      RETURN, TOTAL(2D * h * (7D * (z[ii-4] + z[ii]) + $
                    32D * (z[ii-3] + z[ii-1]) + 12D * z[ii-2]) / 45D, /DOUBLE)

END