;+
; NAME:
;       FMUXBLEED_BCD
;
; PURPOSE:
;       Fixes the 'muxbleed' effect in IRAC images, works on BCDs
;
; INPUTS:
;       A 2D image array
;
; KEYWORD PARAMETERS:
;
;       TITLE      - prefix to output fits name; default='mux_'
;       MUXLEV     - flux level at which to stop using default correction parameters
;	PRELAUNCH  - switch; use prelaunch muxbleed correction parameters
;       DISPLAY    - switch; display resulting image  (not implemented)
;       WRITE      - switch; write out resulting fits image+header
;       HELP       - switch; returns quick keyword summary
;       VERBOSE    - switch; give extra information
;
; OUTPUTS:
;       Returns muxbleed-corrected 2D arrays, optionally writing out
;       corrected fits image, and/or displaying to window.
;
; COMMON BLOCKS:
;       None.
;
; RESTRICTIONS:
;       Does not work for Channel 2 (yet), returns -1 in that case.
;	(Not intended for Channels 3 and 4).
;
; PROCEDURE:
;       fim = fmuxbleed(image.fits, TITLE = title, MUXLEV = muxlev, $
;             /PRELAUNCH, /DISPLAY, /WRITE, /HELP, /VERBOSE)
;
; COMMENTS:
; 	Uses a hybrid approach to correct for muxbleed:
;
;	1.  10k < counts in offending pixel < muxlev [DEFAULT = 30k]:
;	   --> use counts in offending pixel to predict and subtract
;	       level of muxbleed correction required, either using 
;	       prelaunch or irac-R coefficients
;	
;	2.  muxlev < counts in offending pixel < 50k:
;	   --> use image to derive level of muxbleed correction required
;
; PROCEDURES USED:
;	LITERSTAT
;       MUXBLEED_LOOKUP 
;	MUXBLEED_SF
;
; MODIFICATION HISTORY:
;       DS '03nov14 - written
;       DS '04may25 - fixed problem with offending pixels in last two rows
;       SJC 04sep09 - commented out display option and getcolor() call
;       SJC 04sep09 - added conversion of data to counts, uses BCDs
;-
function fmuxbleed_bcd, file,                    $
                    MUXLEV = muxlev,         $
                    PRELAUNCH = prelaunch,   $
; Commented out display option 09 Sep 2004 SJC
;                    DISPLAY = display,       $
                    WRITE = write,           $
                    TITLE = title,           $
                    HELP = help,             $
                    VERBOSE = verbose

     if n_elements(help) ne 0 then begin
         print,' '
         print,'fim = fmuxbleed_bcd("image.fits", TITLE = title, MUXLEV = muxlev, $'
;         print,'      /PRELAUNCH, /DISPLAY, /WRITE, /HELP, /VERBOSE)'
         print,'      /PRELAUNCH, /WRITE, /HELP, /VERBOSE)'
         print,' '
         return,-1
     endif

     if n_elements(file) eq 0 then begin
         print,'% FMUXBLEED: No fits image provided'
         return,-1
     endif
     if size(file,/type) ne 7 then begin
         print,'% FMUXBLEED: Please provide a fits image'
         return,-1
     endif

     if n_elements(PRELAUNCH)  eq 0 then PRELAUNCH = 0
     if n_elements(MUXLEV)     eq 0 then MUXLEV = 30000.
     if muxlev ge 50000. or muxlev le 10000. then begin
	 print,'% FMUXBLEED: Requires that 10k < MUXLEV < 50k'
	 return,-1
     endif


;    prep work:
; Commented out 09 Sep 2004 SJC
;     colors=getcolor(/load)

;    read in image and split into 4 readouts:
     print,'% FMUXBLEED: Doing muxbleed correction on '+file
     image=readfits(file,hd)
; Convert image to counts
     scale = sxpar(hd, 'GAIN') * sxpar(hd, 'EXPTIME') / sxpar(hd, 'FLUXCONV') 
     ptr = where(image eq image)
     image[ptr] = image[ptr] * scale

     if fix(sxpar(hd,'CHNLNUM')) eq 2 then begin
         print,'% FMUXBLEED:  This routine is not yet set up for Channel 2'
         return,-1
     endif else if fix(sxpar(hd,'ACHANID')) gt 2 then begin
         print,'% FMUXBLEED:  This routine is not intended for Channels 3 or 4'
         return,-1
     endif
     readout_in = fltarr(4, 258.*256./4.)
     k = 0.
     for j = 0, 255 do begin
       for i = 0, 255, 4 do begin
         for n = 0, 3 do begin
	   readout_in[n,k] = image[i+n,j]
	 endfor
         k = k + 1
       endfor
     endfor
     readout_out = readout_in

;    cycle through the 4 readout channels:
     for m = 0,3 do begin

;      determine sky level:
       literstat,readout_in[m,*],readstat,/silent
       if n_elements(VERBOSE) ne 0 then $
         print,'% FMUXBLEED: Sky level in readout'+strcompress(m)+ $
	   ' is'+strcompress(readstat.median)+' +/-'+strcompress(readstat.sigma)

;      identify pixels with 10k < counts < 50k: 
       highct = where(readout_in[m,*] gt 10000 and readout_in[m,*] lt 50000)
       if n_elements(VERBOSE) ne 0 then $
         print,'% FMUXBLEED: Number pixels in readout'+strcompress(m)+ $
           ' with 10k < counts < 50k ='+strcompress(size(highct,/dimen))

;      cycle through high-count, offending pixels:
       nmax = size(highct,/dimen) 
       if nmax[0] ne -1 then begin
       for n = 0, nmax[0]-1 do begin

;	 setting up files:
         tmparr = readout_in[m,highct[n]:highct[n]+127]
	 corr_lowct = findgen(128)
	 corr_highct = findgen(64)
 
; 	 for muxlev < highct < 50k, determine scale factor for lookup table:
;        - do statistics on pixels 11:53 past each high pixel on original array
;        - median of lookup table in same region is 13.15 +/- 3.55 
;	 - however, using empirically derived scale of 16.5 below
         tmparr_highct = readout_in[m,highct[n]+11:highct[n]+53] 
         literstat,tmparr_highct,readstat_highct,/silent
	 scale_highct = ((readstat_highct.median - readstat.median) / 16.5) > 0

	 ; hybrid correction: 
	 ; 10k < high counts < muxlev:
	 if tmparr[0] lt muxlev then begin
            if n_elements(VERBOSE) ne 0 then $
	        print,'% FMUXBLEED: '+strcompress(n+1)+' - counts less than'+ $
                strcompress(muxlev)+' for pixel'+strcompress(highct[n])+ $
                ' in readout'+strcompress(m)
	    for i = 0, 127 do begin
	       if prelaunch eq 1 then $
                  corr_lowct[i] = muxbleed_lookup(i) / muxbleed_sf(tmparr[0],/prelaunch)
	       if prelaunch ne 1 then $
		  corr_lowct[i] = muxbleed_lookup(i) / muxbleed_sf(tmparr[0])
	       readout_out[m,highct[n]+i] = readout_in[m,highct[n]+i] - corr_lowct[i]
	    endfor
	 endif else begin
	 ; muxlev < high counts < 50k:
         if n_elements(VERBOSE) ne 0 then $
                print,'% FMUXBLEED: '+strcompress(n+1)+' - counts greater than'+ $
                strcompress(muxlev)+' for pixel'+strcompress(highct[n])+ $
                ' in readout'+strcompress(m)
	    for i = 0, 63 do begin
	       corr_highct[i] = (muxbleed_lookup(i)*scale_highct) > 0
	       readout_out[m,highct[n]+i] = readout_in[m,highct[n]+i] - corr_highct[i]
	    endfor
	 endelse

;    end of cycling through all the high-count pixels in a single readout:
     endfor
     endif
;    end of cycling through the 4 readouts:
     endfor


;    turn 4 readouts back into images:
     image_out = fltarr(256., 256.)
     k = 0.
     for j = 0, 255 do begin
       for i = 0, 255, 4 do begin
         for n = 0,3 do begin
           image_out[i+n,j] = readout_out[n,k]
         endfor
         k = k + 1
       endfor
     endfor

;   possibly display resultant image:
; Commented out 09 Sep 2004 SJC
;    if n_elements(display) ne 0 then begin
;         literstat,image_out,imstat,/silent
;         loadct,0
;         tvlct,r,g,b,/get
;         tvlct,reverse(r),reverse(g),reverse(b)
;         plotimage,bytscl(image_out,$
;             min=imstat.median-   imstat.sigma,$
;             max=imstat.median+5.*imstat.sigma),/preserve
;     endif

; Convert image back to MJy/st
     ptr = where(image_out eq image_out)
     image_out[ptr] = image_out[ptr] / scale

;    possibly write resultant image:
     if n_elements(write) ne 0 then begin
         sxaddhist,'muxbleed corrected - FMUXBLEED.PRO',hd
         if n_elements(title) eq 0 then title='mux_'
         outfile=strmid(file,0,strpos(file,'/',/reverse_search)+1)+title+$
                 strmid(file,strpos(file,'/',/reverse_search)+1,100)
         writefits,outfile,image_out,hd
     endif

     return,image_out
     end