FUNCTION TSUM,X,Y,IMIN,IMAX, NAN=NAN              ;Trapezoidal summation
;+
; NAME:
;       TSUM
; PURPOSE:
;       Trapezoidal summation of the area under a curve. 
; EXPLANATION:
;       Adapted from the procedure INTEG in the IUE procedure library.  
;
; CALLING SEQUENCE:
;       Result = TSUM(y)
;              or
;       Result = TSUM( x, y, [ imin, imax, /nan ] )  
; INPUTS:
;       x = array containing monotonic independent variable.  If omitted, then
;               x is assumed to contain the index of the y variable.
;               x = lindgen( N_elements(y) ).
;       y = array containing dependent variable y = f(x)
;
; OPTIONAL INPUTS:
;       imin = scalar index of x array at which to begin the integration
;               If omitted, then summation starts at x[0].
;       imax = scalar index of x value at which to end the integration 
;               If omitted then the integration ends at x[npts-1].
;       nan: If set cause the routine to check for occurrences of the IEEE 
;                 floating-point values NaN or Infinity in the input data. 
;                 Elements with the value NaN or Infinity are treated as missing data
;
; OUTPUTS:
;       result = area under the curve y=f(x) between x[imin] and x[imax].
;
; EXAMPLE:
;       IDL> x = [0.0,0.1,0.14,0.3] 
;       IDL> y = sin(x)
;       IDL> print,tsum(x,y)    ===>  0.0445843
;       
;       In this example, the exact curve can be computed analytically as 
;       1.0 - cos(0.3) = 0.0446635     
; PROCEDURE:
;       The area is determined of individual trapezoids defined by x[i],
;       x[i+1], y[i] and y[i+1].
;
;       If the data is known to be at all smooth, then a more accurate
;       integration can be found by interpolation prior to the trapezoidal
;       sums, for example, by the standard IDL User Library int_tabulated.pro.
; MODIFICATION HISTORY:
;       Written, W.B. Landsman, STI Corp. May 1986
;       Modified so X is not altered in a one parameter call Jan 1990
;       Converted to IDL V5.0   W. Landsman   September 1997
;       Allow non-integer values of imin and imax  W. Landsman April 2001
;       Fix problem if only 1 parameter supplied W. Landsman June 2002
;       Added /nan keyword. Julio Castro/WL May 2014
;-
; Set default parameters
 On_error,2
 npar = N_params()
   
 if npar EQ 1 then begin
    npts = N_elements(x)
    yy = x
    xx = lindgen(npts)
    ilo = 0   & imin = ilo
    ihi = npts-1 & imax = ihi
 endif else begin

   if ( npar LT 3 ) then imin = 0
   npts = min( [N_elements(x), N_elements(y)] )
   if ( npar LT 4 ) then imax = npts-1
   ilo = long(imin)
   ihi = long(imax)
   xx = x[ilo:ihi]
   yy = y[ilo:ihi]
   npts = ihi - ilo + 1
 endelse   
; 
;  Remove NaN values
;
   if keyword_set(NaN) then begin 
   g = where(finite(yy),npts)
   yy = yy[g]
   xx = xx[g]
  endif          
;   
; Compute areas of trapezoids and sum result
;
  xdif = xx[1:*] - xx
  yavg =  ( yy[0:npts-2] + yy[1:npts-1] ) / 2.  
  sum = total( xdif*yavg ) 

; Now account for edge effects if IMIN or IMAX parameter are not integers

  hi = imax - ihi
  lo = imin - ilo
  if (ihi LT imax) then sum +=  (x[ihi+1]-x[ihi])*hi* $
              (y[ihi] + (hi/2.) *(y[ihi+1] - y[ihi]) )
  if (ilo LT imin) then sum -=  (x[ilo+1]-x[ilo])*lo* $
              (y[ilo] + (lo/2.) *(y[ilo+1] - y[ilo]) )
  return, sum

  end