; NAME:
;      BADPIXLOC
;
; PURPOSE: 
;      Returns a bad pixel mask
;
; CALLING SEQUENCE:
;      Result = BADPIXLOC (nx, ny, nbad)
;
; INPUTS:
;      Nx   = Size in x direction
;      Ny   = Size in y direction
;      Nbad = Number of bad pixel events
;
; OUTPUT:
;      Returns a 2D array with bad pixel values set to 0, and good
;      values set to 1
;      
; PROCEDURE:
;      Make a mask with simple values.  Then create a 3 trees for more
;      complicated values: one for the number, one for the sum of it and
;      previous events in its generation, and one that contains the
;      actual locations of those events.  The tree is then filled, and
;      at the locations set the mask to 0
;
; MODIFICATION HISTORY:
;      created Jun 14 2003 by John Dermody

function badpixloc, nx, ny, nbad, count=count

; make a wider mask
mx = nx + 3
my = ny + 3
mask = bytarr(mx, my) + 1
count = 0

; constants for frequecy of continuation
f = [0.3, 0.5, 0.5] * 0.5 ; correction for 2 possible directions

if (nbad eq 0) then return, mask [0 : nx - 1, 0 : ny - 1]

badloc = floor(randomu(seed, nbad) * (nx)) + $
         floor(randomu(seed, nbad) * (ny)) * mx

numtree = intarr(4, 9)
numtree[0,1] = nbad
sumtree = numtree
loctree = lonarr(4, nbad)
loctree[0,*] = badloc

mask[loctree[0, *]] = 0
for i = 0, 2 do begin
  for j = 1, 2 ^ i do begin
    s0 = sumtree[i+1, j*2-2]
    if (numtree[i, j] gt 0) then begin
      n1 = randomu(seed, binomial=[numtree[i,j], f[i]])
    endif else begin
      n1 = 0
    endelse
    s1 = s0 + n1
    if (n1 ne 0) then loctree[i+1, s0] = loctree[i, s0 : s1 - 1] + 1
    numtree[i+1, j*2-1] = n1
    sumtree[i+1, j*2-1] = s1
    if ((numtree[i, j] - n1) gt 0) then begin
      n2 = randomu(seed, binomial=[(numtree[i,j] - n1), f[i] / (1.- f[i])])
    endif else begin
      n2 = 0
    endelse
    s2 = s1 + n2
    if (n2 ne 0) then loctree[i+1, s1] = loctree[i, s1 : s2 - 1] + mx
    numtree[i+1, j*2] = n2
    sumtree[i+1, j*2] = s2
  endfor
  if (s2 ne 0) then mask[loctree[i+1, 0 : s2 - 1]] = 0
endfor

mask = mask[0: nx-1, 0:ny-1]
count = n_elements(where(mask eq 0))

return, mask
end