;
pro pixel,wave,flux,xstart,xend,pix
;+
; PIXEL:
;
;  Procedure to measure pixel index
;  of spectral lines using cursor input of center
;  of emission or absorbtion line.
;
;   Uses TABINV from the Astron. User's Lib.
;
; INPUT 2-dimensional vector of index and flux
; OUTPUT pixel index (in units of x-axis scale)
;
; Written by J.E. Neff 31 Jan 92. Modified from equiv_width.pro
; by James Lynch 25 Jan 93
; Freeze-dried and reconstituted by Robin Ciardullo, Groundhog's day 1995
;
; 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()
sz = size(wave)
if np eq 0 or sz(0) eq 0 then goto,ERRMES
if sz(0) eq 2 and sz(1) gt 1 then begin
    if np eq 1 then begin
	xstart = 0
	xend = sz(2)-1
    endif
    pix = ' '
    if np eq 4 then begin
	if xend eq 'w' then pix='w'
    endif
    if np ge 3 then begin
	xend = xstart
	xstart = flux
    endif
    if pix eq 'w' then begin
	tabinv,wave(0,*),xstart,i1
	xstart=i1
	tabinv,wave(0,*),xend,i2
	xend=i2
    endif
    ymax = max(wave(1,xstart:xend))
    ymin = min(wave(1,xstart:xend))
endif else begin
    sz(0) = 1
    if np eq 2 then begin
	xstart = 0
	xend = n_elements(wave)-1
	pix = ' '
    endif
    if pix eq 'w' then begin
        tabinv,wave,xstart,i1
	xstart = i1
        tabinv,wave,xend,i2
	xend = i2
    endif
    ymax = max(flux(xstart:xend))
    ymin = min(flux(xstart:xend))
endelse
;
ydif = ymax - ymin
ymin = ymin - (0.05 *ydif)
ymax = ymax + (0.05 *ydif)
!psym=10 & !linetype=0

if sz(0) eq 2 then begin
    plot,wave(0,xstart:xend),wave(1,xstart:xend),yrange=[ymin,ymax],ystyle=1 
endif else begin
    plot,wave(xstart:xend),flux(xstart:xend),yrange=[ymin,ymax],ystyle=1 
endelse

print, 'Position cursor at center of feature and click left mouse button'
print, string(7B)
repeat begin
    cursor,x1,y1,1
    temp=get_kbrd(0)
    if (temp eq 's') or (temp eq 'S') then goto,DONE
    if sz(0) eq 2 then begin
	tabinv,wave(0,*),x1,i1
        pix=wave(0,i1)
    endif else begin
	tabinv,wave,x1,i1
        pix=wave(i1)
    endelse
;
    print,' '
    print,'The PIXEL INDEX of this feature is:',pix
    print,' '
    print, 'Next feature; otherwise type S, then click to STOP',string(7B)
endrep until 0
DONE:
return
ERRMES:
print,'ERROR! Call with a 2-dimensional array or 2 1-dimension arrays'
print,string(7B)
return
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