function debleed, imgs, bpar0, bpar1, $
                  nfix = nfix, nowrap = nowrap, $
                  badpix = bp, masks = masks, $
                  silent = silent, help = help

;+
; Correct the bleeding in a NIRC image assuming exponential decay.
; For a pixel at location s0 with flux f(s0), the decay model assumes
; this pixel contributes to a pixel at location s1 which is
; subsequently readout:
;
;       f_bleed(s1) = f(s0) * c0 * exp[c1*(s1-s0)]
; where
;       c0    = normalization
;      -1./c1 = decay length scale
;       c2    = residual bias level in the data
;
; Note that NIRC array is readout from left to right and bottom to top.
;
; Can have different coeffs for each readout or just use same coeffs
; for all 4 readouts (latter runs much faster)
;
; By default, the debleeding wraps from the top of image to the
; bottom of the image (set /nowrap to disable). Technically, this is
; incorrect, but in practice, this is useful to prevent bottom edge of
; the debled images from having slightly more counts than the rest of
; the images.
;
; INPUTS
;   imgs    images to be debleeded
;   bpar0   decay coeff c0 
;             a scalar or a 4-element vector (for each readout)
;   bpar1   decay coeff c1 for decay (in units of every 4th
;           pixel, i.e. per readout) 
;             a scalar or a 4-element vector (for each readout)
;
; OPTIONAL KEYWORD INPUTS
;   nfix    length of bleeding to correct (per readout, 
;             i.e., nfix=64 means correct 256 pix length in image) 
;   nowrap  do not wrap from top of image to bottom when debleeding
;   badpix  bad pixel mask (0=bad, !0=good)
;   masks   masks for each image (0=bad, !0=good)
; 
; RETURNS
;   outimgs debled images
;
; USES
;   strc
;
; HISTORY: Written by M. Liu (UCB) 11/12/96 
; 02/03/97 (MCL): much faster (vectorized) algorithm 
;                 able to handle multiple images
; 02/10/97 (MCL): accepts coeffs as either scalar or 4-vectors
;                 wraps from top to bottom of the image (added /nowrap)
; 02/19/97 (MCL): bad & masked pixels are unchanged in output image
; 09/29/00 (MCL): prints 4 updates per image (instead of 8)
;
; Please send comment/questions to mliu@astro.berkeley.edu
;-

BADVAL = -1e6
n0 = n_elements(bpar0)
n1 = n_elements(bpar1)
if (n0 eq 0) then bpar0 = 0.0027
if (n1 eq 0) then bpar1 = -1.0/16.31
if n_params() lt 1 or n_elements(imgs) eq 0 or keyword_set(help) then begin
    print, 'function debleed(imgs,{bpar0},{bpar1},'
    print, '                 [nfix=],[nowrap],[badpix=],[masks=],[inlist=],[silent],[help])'   
    print, '     default c0 = ', strc(bpar0)
    print, '          -1/c1 = ', strc(-1/bpar1), ' pix/readout'
    retall
endif


; initialize output
sz = size(imgs)
if (sz(0) eq 3) then npics=sz(3) else npics=1 
xs = sz(1)
ys = sz(2)
np = sz(4)
np1 = np-1
outimgs = imgs


; sanity checks of masks
if not(keyword_set(bp)) then bp =  bytarr(xs, ys) + 1B $
else begin
    szb = size(bp)
    if total((sz-szb)(1:2)) ne 0 then begin
        message, 'input and badpix images not the same size', /info
        retall
    endif
endelse
if keyword_set(masks) then begin
    szm = size(masks)
    if total((sz-szm)(0:2)) ne 0 then begin
        message, 'input images and masks not the same size', /info
        retall
    endif
endif


; initialize cofficients: 
; check if correction will be done per readout (rdflag=1) or 
;   if treat all readouts the same (rdflag=0)
rdflag = 0
if (n0 ne 1) and (n1 ne 1) then begin
    if (total(bpar0 - bpar0(0) + bpar1 - bpar1(0)) ne 0) then begin
        bscl = fltarr(sz(1), sz(2))
        rdflag = 1 
    endif else begin
        bpar0 = bpar0(0)
        bpar1 = bpar1(0)
    endelse
    if not(keyword_set(silent)) then $
      message, 'correcting each readout separately', /info
endif else if not(keyword_set(silent)) then $
  message, 'treating all readouts the same', /info
if (bpar1(0) ge 0) then begin
    message, 'invalid c1 = '+strc(bpar1(0)), /info
    retall
endif


; establish number of pixels to fix PER READOUT
if not(keyword_set(nfix)) then nfix =  xs/4
xx =  findgen(nfix)
xx4 = xx * 4.
if not(keyword_set(silent)) then begin
  message, 'correcting ' + strc(nfix) + ' pixel bleeding per readout ('$
    + strc(nfix*4) +' pix total)', /info
  message, strc(npics)+' images to debleed', /info
endif


; loop over images and correct
for k = 0, npics-1 do begin

    if not(keyword_set(silent)) then $
      print, form = '($,"image ",A,": ")', strc(k)

    ; copy of original image with bad pixels set to 0
    ; neglect bad pixels in the bleeding correction since
    ; their bleeding will cancel during sky subtraction
    out = imgs(*, *, k)
    if keyword_set(masks) then begin
        immask = (bp * masks(*, *, k)) ne 0
        orig = out * immask
    endif else begin
	immask = (bp ne 0)*(out ne BADVAL)
        orig = out * immask
    endelse
      
    ; do correction by shifting the entire image in steps of 4 pixels (for
    ; the 4 readouts), multipying by the appropriate bleed factor for
    ; each shift, and then subtracting from the data
    for i = 1, nfix do begin
        
        ; shift the image
        ; if /nowrap, mask out the pixels at the bottom
        imsh = shift(orig, i*4)
        if keyword_set(nowrap) then $
          imsh(indgen(4*i)) =  0

        ; determine bleed factor for this shift for each readout
        if (rdflag eq 1) then $
          for j = 0, 3 do $
          bscl(xx4+j, *) = bpar0(j) * exp(bpar1(j)*i) $
        else $
          bscl = bpar0 * exp(bpar1*i)

        ; subtract the bleed for pixels with this shift
        out = out - bscl * imsh

        if (i mod (nfix/4)) eq 0 and not(keyword_set(silent)) then $
          print, form = '($,A," ")', strc(4*i)

    endfor
    if not(keyword_set(silent)) then print

    ; leave the bad and/or masked pixels unchanged in output image
    wbad =  where(immask eq 0, nbad)
    if (nbad gt 0) then $
      out(wbad) = (imgs(*, *, k))(wbad)

    ; save to output array
    outimgs(*, *, k) =  out

endfor



;----------------------------------------------------------------------
; OLD SCHEME: more intuitive, much much slower!    
;   treat image as a 1-d array for correction
;   note: don't need to correct the last 4 pixels
;----------------------------------------------------------------------
; calculate (normalized) bleed for each readout
;nbleed = fltarr(nfix, 4)
;for i = 0, 3 do $
;  nbleed(*, i) = bpar0(i) * exp(bpar1(i)*(xx+1))
;wgood =  where(bp ne 0, ngood)
;if not(keyword_set(silent)) then $
;  message, 'pixels to debleed = '+strc(ngood), /info
;for i = 0L, (ngood-1)-4 do begin

;    ; determine which readout line
;    rr =  wgood(i) mod 4

;    ; get the pixel
;    pix = img(wgood(i))
    
;    ; make the bleed
;    bleed = pix * nbleed(*, rr)
    
;    ; subtract the bleed
;    bb = wgood(i) + xx4 + 4
;    wg =  where(bb le np1, ng)
;    outimg(bb(wg)) =  temporary(outimg(bb(wg))) - bleed(wg)

;    if (i mod 1000) eq 0 and not(keyword_set(silent)) then $
;      print, form = '($,A," ")', strc(i)

;endfor
;if not(keyword_set(silent)) then print


return, outimgs


end