;+
;NAME: ivm_sub [_4.0]
;
;PURPOSE: substitute good-seeing frames for bad ones, correcting
;	  for intensity shifts caused by local velocity signals
;
;CATEGORY: IVM data reduction
;
;CALLING SEQUENCE: ivm_sub, oblk, rmssee, filein, fileout, indat=idatps, $
;			outdat=idatws, outfile=datws, infile=datw
;
;INPUTS:
;	hstr = oblk header structure
;	seearr = array of seeing values for the images (as from seerms.pro)
;	nn = npix per side (assumed square)
; 	nw = number of wavelength points in the spectral line, incl.continuum
;	nmod = number of modulation sets per spectral point, usually 6 or 4
;	nr = #repeats - number of times spectra was scanned, same as in rivm
; 	hstr = oblk header structure
; 	filein  = If set, use 'infile' file input, else use 'indat' array
; 	fileout = If set, use 'outfile' file output, else use 'outdat' array
;
;OPTIONAL KEYWORD PARAMETERS:
;	nr = Number of repeats - use to restrict substitution to first 'nr' only
;	infile = name of data-camera file to be subbed
;	outfile = Output file
;	indat = data-camera array to be subbed
;	outdat = Resultant array of subbed data
;	verbose = Verbosity
;
;OUTPUTS:   outfile = file of images written to disk.
;
;COMMON BLOCKS: none
;
;SIDE EFFECTS: Writes a file to disk, if fileoutput set.
;
;RESTRICTIONS: Images which have a seeing value (nominally from
;		seerms.pro, but from whatever parameter you wish to use)
;		less than median-stdev are the ones which are substituted.
;		This is hard-wired, but could be changed in the future.
;
;PROCEDURE:
;	warning - makes corrections to all three nr sets of nmod
;		at once and writes them to a temporary file,
;		then lastly reads things back in & writes them out in the
;		correct order.
;
;MODIFICATION HISTORY:
; 	written 12/97 KDL, categorizing the seeing-quality in to 'good',
;		'mediocre' or 'bad', and constructing intensity-offsets
;		only between images of like-seeing.  
;	12/15 - the 3-category algorithm was a bit too restrictive.
;		changed to construct intensity offsets from all images, no matter
;		what the seeing-quality-match was. rewrote for 'case'
;		statements vs. if/then.  
;	KDL 10/02 - small debug for when no frames with seeing lt med-stdev...,
;		made more standard w.r.t. using header structure.
;	ELW 20091216 - Moved to the realm of memory, with file I/O optional.
;	ELW 20091217 - Use if/else rather than case; eliminate reordering loop.
;	ELW 20100528 - Added 'filein' parameter, comments
;	ELW 20110718 - Added nr, verbose parameters
;	ELW 20110719 - Added KDL diagnostics
;-

pro ivm_sub, hstr, seearr, filein, fileout, nr=nr, indat=indat, infile=infile, $
				outfile=outfile, outdat=outdat, verbose=verbose

pn='ivm_sub: '
vb = (keyword_set(verbose)) ? verbose : 0
doprint= (vb ge 4) ? 1 : 0 	; Set for extreme (>4) verbosity only
if (vb gt 0) then print, 'ivm_sub:  start - ', systime()
nn=hstr.pps
nmod=hstr.nframes	; number of modulation sets per spectral point
nw=hstr.n_image_sets	; number of wavelength points in the spectral line
if (keyword_set(nr) && (vb gt 1)) then $
    print, 'ivm_sub:  nr keyword set to ', strtrim(nr,2)
if (~keyword_set(nr)) then $ 
    nr=hstr.ntimes	; number of repeats = how many times spectra was scanned

; Open & read input file if 'infile' is set, and array keyword is not
if (filein) then begin
    if ~keyword_set(infile) then begin
        print, 'ivm_sub: ERROR! Filename "infile" must be set.  Stopping.'
        stop
    endif
    indat = fltarr(nn,nn,nmod, nw*nr)
    openr,lun1,infile,/get
    readu, lun1, indat
    close,lun1 & free_lun,lun1
endif else if (~keyword_set(indat)) then begin
    print, 'ivm_sub: ERROR! indat array must be set.  Stopping.'
    stop
endif

med=median(seearr)
st=stdev(seearr)
bad=where(seearr lt med-st,badcount)
if (badcount lt 1) then begin
    print, 'ivm_sub:  Skewed seeing distribution:  no frames with seeing lt median-stdev, exiting.'
    outdat=indat
    if (fileout eq 1) then begin
	openw,lun2,outfile,/get_lun
        writeu,lun2,indat
	close,lun2 & free_lun,lun2
    endif
    return
endif

; Images are categorized as acceptable (1) or bad (0)
okarr=intarr(nmod*nw*nr)+1
okarr(bad) = 0 ; these definitions can be changed...
okarr  =reform(okarr, nmod*nw, nr) ; should be e.g. 0-79,80-159,160-239
rseearr=reform(seearr,nmod*nw, nr)
datarr=fltarr(nn,nn,nmod,nr)
outarr=fltarr(nn,nn,nmod,nr)
okm=intarr(nmod,nr)
seem=fltarr(nmod,nr)

; Intensity-correction matrix  - intensity offsets are spatial due to 
;   oscillations.  they are between nrepeats, i.e. the diff between repeat 0,1
;   at the wavelength k that we're at
ioff=fltarr(nn,nn,nr,nr)

jarr=indgen(nr)		; used for indexing below

outdat = fltarr(nn,nn,nmod, nw*nr)	 ; Output array
subarr=intarr(nw,nmod,nr) 	; diagnostic of which got substituted where

for k=0,nw -1 do begin	;loop over spectral positions
   for j=0,nr-1 do begin    ; loop over spectral repeats
      datarr(0,0,0,j)=indat(*,*,*,j*nw+k)  ; load in conformal data from each repeat - was (0,0,0,j)
      okm(*,j)=(okarr([k*nmod+indgen(nmod)],j)) 	; and the "quality" in the
      seem(*,j)=(rseearr([k*nmod+indgen(nmod)],j)) 	;  same format
   endfor

   replicate_inplace, ioff, 0.		; reinitialize

   ; Calculate entries for the intensity-offsets
   for j=0,nr-1 do begin
      ; ioffs are average differences between the nmod images in each nr,
      ; smoothed to remove small-scale intense features.  Tried
      ; median-smoothing, but it left some sharp-edged features.
      ioff(0,0,j,(shift(jarr,-1))(j))=smoothe((total(datarr(*,*,*,j),3)-$
	     total(datarr(*,*,*,(shift(jarr,-1))(j)),3))/float(nmod),3)
      ; corresponding diagonal element is (-1)* original, i.e. 2->1 is (-1)*1->2
      ioff(0,0,(shift(jarr,-1))(j),j)=(-1.)*ioff(*,*,j,(shift(jarr,-1))(j))
   endfor

   ; Begin substituting and ioff-applying when needed 
   for j=0,nr-1 do begin
     subarr(k,*,j)=j		; diagnostic
      for i=0,nmod-1 do begin
	 targetj=(reverse(sort(seem(i,*))))(0) ; order the images' seeing and
	 targetoj=(okm(i,targetj))	    ; whether they're considered good.
	 if (okm(i,j) ge 1) then begin
	    outarr(0,0,i,j)=datarr(*,*,i,j)    ; if it's good, just copy it
	 endif else if (targetoj gt 0) then begin
	    ; if image is bad and there's a corresponding good image (i.e.
	    ; a good targetj) then sub in targetj image and ioff-correct it
	    outarr(0,0,i,j)=datarr(*,*,i,targetj)+ioff(*,*,j,targetj) 
	    subarr(k,i,j)=targetj		; diagnostic
	 endif else begin
	    ; if image is bad but targetj image is also bad, just use
	    ; use original bad image (since targetj is a sort of 
	    ; the seeing-parameter, next in line won't help any...
	    outarr(0,0,i,j)=datarr(*,*,i,j)
	 endelse
      endfor ; nmod loop

      outdat(0,0,0,k+j*nw) = outarr(*,*,*,j)	; Set output array cube

   endfor ; nr loop
   if doprint then begin
     print,'Substitution matricies:'
     print, 'I    Repeat Number'
     print, 'Q'
     print, 'U'
     print, 'V'
     print,reform(subarr[k,*,*]) 
     print,'-----'
   endif
endfor ; end of nw loop

; Tally up for warnings:
badtot=fltarr(nr)
for k=0,nw-1 do begin
   for j=0,nr-1 do begin
       bb=where(subarr[k,*,j] ne j,count)
       badtot[j]=badtot[j]+count
   endfor
endfor
; hard limit for warning is more than of points substituted:
for j=0,nr-1 do begin 
   if (vb ge 1) then print,pn,'Repeat #: ',strtrim(j,2),' ,',strtrim(reform(badtot[j]),2)
   if reform(badtot[j])/(nw*nmod) gt float(nw)/(2.*nw*nmod) then print,pn,'WARNING: TOO MANY in repeat #',strtrim(j,2)
endfor

if (fileout && keyword_set(outfile)) then begin
   openw,lun2,outfile,/get
   writeu,lun2,outdat
   close,lun2 & free_lun,lun2
endif

end