Viewing contents of file '../idllib/ghrs/pro/blemish.pro'
;************************************************************************
;+
;*NAME:
;  BLEMISH     (General IDL Library 01) August 13, 1986
;
;*CLASS: 
;    Data Editing
;
;*PURPOSE:  
;    To interactively remove blemishes in a flux vector by linear
;    interpolation between the flux values on each side of the 
;    blemish. (Based on BLEMISH by D. Lindler.)
;
;*CALLING SEQUENCE: 
;    BLEMISH,WAVE,FLUX,BSIZE,FCOR
;
;*PARAMETERS:
;    WAVE   (REQ) (I) (1) (F D)
;           Required input vector containing the wavelength data associated
;           with the FLUX data which are to be edited.
;
;    FLUX   (REQ) (I) (1) (F D)
;           Required input vector containing the flux data which are to be
;           edited.
;
;    BSIZE  (OPT) (I) (0) (F D)
;           Optional input scalar parameter specifying the interval, in the same
;           units as WAVE, over which the FLUX data are to be averaged 
;           on either side of the blemish, before interpolating across the
;           blemish. 
;
;    FCOR   (REQ) (O) (1) (F D)
;           Required output vector giving the edited FLUX array.
;
;*EXAMPLES
;    To linearly interpolate across a blemish, using the adjacent good 
;    data points:
;    BLEMISH,WAVE,FLUX,FCOR
;
;    To linearly interpolate across a blemish using an average of the 
;    data points within BSIZE wavelength units of the edges of the blemish,
;    BLEMISH,WAVE,FLUX,BSIZE,FCOR
;
;*SYSTEM VARIABLES USED:
;    !ERR
;    !NOPRINT
;
;*INTERACTIVE INPUT:
;    Flux versus wavelength is plotted and the user is requested
;    to indicate with the verticle cross-hairs the extent of the
;    blemish.  The procedure then fills those points within the
;    blemish with values interpolated between the endpoints of the
;    blemish. 
;    This process is repeated until the user types 0 (zero) or presses
;    the right mouse button when asked to position the crosshairs.
;
;*SUBROUTINES CALLED:
;    BINS
;    TABINV
;    PARCHECK
;
;*FILES USED:
;    None.
;
;*SIDE EFFECTS:
;    Aborting the procedure may result in !ERR and !NOPRINT being 
;    altered. 
;    
;*RESTRICTIONS:
;    Device Dependent - This procedure requires a terminal equipped 
;                       with a graphics cursor.
;    (modified to run under unix sun idl version 1.2)
;
;*NOTES:
;        Order of specifying endpoints is not important.
; 
;                
;*PROCEDURE: 
;    For each pair of positions marked by the user the wavelengths
;    at those points are computed by CURSOR. The subscripts
;    associated with those wavelengths are computed with TABINV.
;    BINS is called to calculate binned flux if BSIZE is specified.
;    Let IMIN and IMAX be the pair of subscripts.  For each value
;    of I, IMIN < I < IMAX, a new flux is computed for FCOR(I):
;         
;    FCOR(I) = FCOR(IMIN) + S*(WAVE(I)-WAVE(IMIN))
;    where:
;    S = (FCOR(IMAX)-FCOR(IMIN)) / (WAVE(IMAX)-WAVE(IMIN))
;    or, if BSIZE is specified,
;    S = (MEAN(1) - MEAN(0))  / (WC(1) - WC(0))
;    where MEAN is the average flux values and WC is the bin center
;    positions.
;
;        
;*MODIFICATION HISTORY:
;                 D. Lindler  GSFC  initial version
;    Oct 31 1985  JKF         GSFC  DIDL compatible (indirect compilation)
;    Jul 25 1986  CLI         GSFC  add BSIZE parameter to take average flux
;                                   on either side of the specified points.
;    Aug 13 1986  RWT         GSFC  use N_ELEMENTS and CURSOR commands
;    Feb 19 1987  RWT         GSFC  VAX Mods: make BSIZE optional, account
;                                   references, @ for #, remove one GOTO
;                                   statement and DO loop, improve calculation
;                                   of subscripts, allow points to be 
;                                   specified in any order.
;
;    Mar 14 1987  RWT         GSFC  add PARCHECK
;    May  6 1987  RWT         GSFC  fix error with N_PARAMS(0)=3 mode
;    Mar  8 1988  CAG         GSFC  VAX RDAF-style prolog, add disabling
;                                   of tabular printout from BINS.
;    jan  1 1990  jtb         gsfc  modified for sun/unix idl.
;    Mar 4 1991      JKF/ACC    - moved to GHRS DAF (IDL Version 2)
;-
;********************************************************************
pro blemish,wave,flux,bsize,fcor
;
npar = n_params(0)
if npar eq 0 then begin
   print,' blemish,WAVE,FLUX,bsize,FCOR'
   retall & end
parcheck,npar,[3,4],'blemish'
pcheck,wave,1,010,0111
pcheck,flux,2,010,0111
;
fcor=flux
npts = long(n_elements(flux))
;
; check for optional bsize parameter (set to 0 if not specified)
;
if n_params(0) lt 4 then bs = 0.0 else bs = float(bsize)
wc = fltarr(2)
width = wc + bs
;
noprint=!noprint    ;save !noprint
!noprint=1          ;set noprint to suppress "bins" output
;
; loop to plot flux with directions for blemish removal
;
while 1 do begin            ;loop until interactive exit
;
plot,wave,fcor
;
if (!d.name ne 'TEK') then begin
   ;
   print,'Press the MIDDLE or LEFT buttons to select each POINT using'
   print,'the VERTICLE CROSSHAIR.'
   print,'Press the RIGHT mouse button to EXIT.'
   ;
   endif else begin
   ;
   print, 'Press the SPACE BAR to select each POINT using'
   print, 'the VERTICAL CROSSHAIR.'
   print, 'Type 0 (no  return) to EXIT.'
   ;
   endelse
   ;
if bs ne 0 then begin
   print,' '
   print,form="(1x,a,f7.2,a)", $
     'The interpolation will be made using mean fluxes in ',bs,' ang bins'
   print,' on either side of the points chosen.'
   end
   ;
; let user define limits to fix region using cursor
;
cursor,w1,f1,/down,/data
flush = get_kbrd(0)
print,string(7b)
if !err eq "60 or !err eq 4 then goto,end_loop
;
cursor,w2,f2,/down,/data
flush = get_kbrd(0)
print,string(7b)
if !err eq "60 or !err eq 4 then goto,end_loop
;
; compute subscripts for area to be corrected
;
tabinv,wave,w1,p1
tabinv,wave,w2,p2
imin = fix((p1<p2)+0.5) < (npts-2)
imax = fix((p1>p2)+0.5) > (imin+2)
;
; find average flux next to specified endpoints (for bs = 0.0)
; tabular printout is suppressed for this section (!noprint = 1)
;
if bs ne 0.0 then begin
   wc([0,1]) = [wave(imin)-bs/2.,wave(imax)+bs/2.]
   bins,wave,flux,wc,width,mean,sigma,np
   s = (mean(1)-mean(0))/(wc(1)-wc(0))
   f0 = mean(0) + s*bs/2.
   end else begin
      s = (fcor(imax)-fcor(imin))/(wave(imax)-wave(imin))
      f0 = fcor(imin)
      end
;    
; fix up bad region
;
fcor(imin+1) = f0 + s*(wave(imin+1:imax-1)-wave(imin))
endwhile
;
;
end_loop:
          if (n_params(0) eq 3) then bsize = fcor
          !noprint = noprint        ;restore !noprint
          ;
return
end