;
pro equiv, wave, flux, xstart, xend, flag
;+
; EQUIV_WIDTH:
;
; Quick and Dirty procedure to measure equivalent
;   widths of spectral lines using cursor input of
;   the left and right-hand side.
;
; Uses the function TSUM (from the Astronomy User's
;   Library) to calcluate the area in the line. 
;   Also uses TABINV from the Astron. User's Lib.
;
; INPUT wavelength and flux vectors
; OUTPUT equivalent width (in units of wavelength scale)
;
; Written by J.E. Neff 31 Jan 92.
;
; Status:
; Not tested under all circumstances.  Seems to work
;   ok for absorption lines from normalized spectra.
;   I haven't verified my memory of the equivalent
;   with algorithms yet, but the answers look sensible.
;-
np = n_params()
if np lt 2 then begin
  print,'call with wavelength and flux vectors as...'
  print,string(7B)
  print,'equiv_width,wavelength,flux,(xstart),(xend),(''w'')'
  print,string(7B)
  print,' '
  return
endif
;
if np eq 2 then begin
  xstart = 0
  xend = n_elements(wave) - 1
endif
if np eq 5 then begin
  if flag eq 'w' then begin
     tabinv,wave,xstart,i1
     tabinv,wave,xend,i2
     xstart=i1
     xend=i2
  endif
endif

!psym=10 & !linetype=0
set_xy,0,0,0,0
plot,wave(xstart:xend),flux(xstart:xend)
print,'Equivalent-Width Measurer -- press "s" and click to stop',string(7B)
repeat begin
    print,'Position cursor and click at continuum level on left side of line ',string(7B)
    cursor,x1,y1,1
    temp=get_kbrd(0)
    if (temp eq 's') or (temp eq 'S') then goto,DONE
    wait,1
    print,'Position cursor and click at continuum limit on  right side of line ', string(7B)
    cursor,x2,y2,1
    temp=get_kbrd(0)
    if (temp eq 's') or (temp eq 'S') then goto,DONE
    tabinv,wave,x1,i1
    tabinv,wave,x2,i2
;
    !psym=0 & !linetype=7
    oplot,[x1,x2],[y1,y2]
    !linetype=1
    for i=i1,i2 do oplot,[wave(i),wave(i)],[flux(i),(y1+y2)/2]
;
    i1=fix(i1) & i2=fix(i2)
    line_area=tsum(wave,flux,i1,i2)
    cont_level=mean([y1,y2])
    cont_area=cont_level*(x2-x1)
    ew=-(line_area-cont_area)/cont_level
    wbar = 0
    wsum = 0
    for i=i1,i2 do begin
	fc = y1 + ((y2 - y1)*(i-i1)/(i2-i1))
	weight = fc - flux(i)
	wbar = wbar + (wave(i) * weight)
	wsum = wsum + weight
    endfor
    wbar = wbar / wsum
    print,'Total Area under line =',line_area
    print,'Total Area under continuum =',cont_area
    print,'Mean continuum Flux level =',cont_level
    print,' '
    print,'Mean Wavelength =',wbar
    print,'The EQUIVALENT WIDTH =',ew
endrep until 0
DONE:
!linetype=0
return
end

FUNCTION TSUM,X,Y,IMIN,IMAX
;+
; NAME:
;       TSUM
; PURPOSE:
;       Trapezoidal tsumration of the area under a curve.   This procedure
;       was formerly know as INTEG
; CALLING SEQUENCE:
;       Result=TSUM(y)
;       Result=TSUM(x,y)
;       Result=TSUM(x,y,imin,imax)    ;Integrate between IMIN and IMAX
; INPUTS:
;       x = array containing independent variable.  If omitted, then
;           x is assumed to contain the index of the y variable.
;           x=indgen(n_elements(y)).
;       y = array containing dependent variable y=f(x)
;       imin = index of x array at which to begin the integration.  If
;              omitted, then tsumration starts at x(0).
;       imax = index of x value at which to end the integration.  If 
;              omitted then the integration ends at x(npts).
; OUTPUTS:
;       result = area under the curve y=f(x) between x(imin) and x(imax).
; PROCEDURE:
;       The area is determined of indivdual trapezoids defined by x(i),
;       x(i+1), y(i) and y(i+1).
; MODIFICATION HISTORY:
;       Written, W.B. Landsman, STI Corp. May 1986
;       Modified so X is not altered in a one parameter call Jan 1990
;-
; Set default parameters
NPAR = N_PARAMS(0)
IF NPAR EQ 1 THEN BEGIN
    NPTS = N_ELEMENTS(X)
    YY = X
    XX = INDGEN(NPTS)
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
   XX = X(IMIN:IMAX)
   YY = Y(IMIN:IMAX)
   NPTS = IMAX - IMIN + 1
ENDELSE
;
; Compute areas of trapezoids and sum result
XDIF = SHIFT(XX,-1) - XX
XDIF = XDIF(0:NPTS-2)
YAVG = (YY+SHIFT(YY,-1))/2.
YAVG = YAVG(0:NPTS-2)
RETURN,TOTAL(XDIF*YAVG)
END

PRO TABINV, XARR, X, IEFF
;+ 
; NAME:
;   TABINV     
; PURPOSE:  
;   To find the effective index of a function value in
;   an ordered vector.
; CALLING SEQUENCE:
;   TABINV, XARR, X, IEFF
; INPUTS:
;        XARR - the vector array to be searched, must be monotonic
;               increasing or decreasing
;        X    - the function value(s) whose effective
;               index is sought (scalar or vector)
; OUTPUT:
;        IEFF - the effective index or indices of X in XARR
;               real or double precision, same # of elements as X
; RESTRICTIONS:
;        TABINV will abort if XARR is not monotonic.  (Equality of 
;        neighboring values in XARR is allowed but results may not be
;        unique.)  This requirement may mean that input vectors with padded
;        zeroes could cause routine to abort.
; PROCEDURE:
;        A binary search is used to find the values XARR(I)
;        and XARR(I+1) where XARR(I) < X < XARR(I+1).
;        IEFF is then computed using linear interpolation 
;        between I and I+1.
;                IEFF = I + (X-XARR(I)) / (XARR(I+1)-XARR(I))
;        Let N = number of elements in XARR
;                if x < XARR(0) then IEFF is set to 0
;                if x > XARR(N-1) then IEFF is set to N-1
; EXAMPLE:
;        Set all flux values of a spectrum (WAVE vs FLUX) to zero
;        for wavelengths less than 1150 Angstroms.
;         
;           IDL> tabinv, wave, 1150.0, I
;           IDL> flux( 0:fix(I) ) = 0.                         
; FUNCTIONS CALLED:
;    ISARRAY
; REVISION HISTORY:
;   Adapted from the IUE RDAF                     January, 1988         
;   More elegant code  W. Landsman                August, 1989
;-               
On_error,2

if N_params() LT 3 then begin
     print,'Calling sequence- TABINV, XARR, X, I'
     return
endif
npoints = N_elements(xarr) & npt= npoints - 1
;
; Initialize binary search area and compute number of divisions needed
;
ileft = intarr( N_elements(x) ) & iright = ileft

ndivisions = fix( alog10(npoints) / alog10(2.0)+1.0 )
;
; Test for monotonicity 
;
i = xarr - shift( xarr,1)
a = where(i GE 0, N)

if ( N EQ npt) then $ ; Increasing array ?

    iright = iright + npt $

else begin

     a = where(i LE 0, N)  ; Test for decreasing array
     if ( N EQ npt ) then ileft = ileft + npt $    
     ELSE message, $
       'ERROR - First parameter must be a monotonic vector' 

endelse          
;
; Perform binary search by dividing search interval in
; half NDIVISIONS times
;
for i = 1,ndivisions do begin    

    idiv = (ileft + iright) / 2      ;Split interval in half
    xval = xarr(idiv)                ;Find function values at center
    greater = ( x GT xval )          ;Determine which side X is on
    less   = ( x LE xval )  
    ileft =  ileft*less + idiv*greater    ;Compute new search area
    iright = iright*greater + idiv*less     

endfor
;
; Interpolate between interval of width = 1
;
xleft =  xarr(ileft)              ;Value on left side
xright = xarr(iright)             ;Value on right side

ieff = (xright-x)*ileft + (x-xleft)*iright + ileft*( xright EQ xleft ) 

ieff = ieff / float( xright - xleft + ( xright EQ xleft ))    ;Interpolate

ieff = ieff > 0.0 < npt        ;Do not allow extrapolation beyond ends

if not ISARRAY(x) then ieff = ieff(0)  ;Make scalar if X was scalar

return
end