pro get_astr_st, spix, ra_cent, dec_cent, x_cent, y_cent, astr

tagnames    = ['cd', 'cdelt', 'crpix', 'crval', 'ctype', 'longpole', $
               'latpole','pv2']
tagdescript = 'd(2,2),d(2),f(2),d(2),a(2),d,d,d(2)'
create_struct_mod, astr, '', tagnames, tagdescript

astr.cd[0,0]  = 1.
astr.cd[1,1]  = 1.
astr.cd[1,0]  = -0.
astr.cdelt[0] = -spix
astr.cdelt[1] = spix
astr.crpix[0] = x_cent
astr.crpix[1] = y_cent
astr.crval[0] = ra_cent
astr.crval[1] = dec_cent
astr.ctype[0] = 'RA---TAN'
astr.ctype[1] = 'DEC--TAN'
astr.longpole = 180.
astr.latpole  = 90.
astr.pv2      = [0., 0.]

end


pro fill_grid, infile, radec_cent, mask_file, rmin, rmax, nbin, img, $
               ra_min, dec_min, ra_max, dec_max, spix, non_masked_img,$
               mask, mask_hdr, xstart = xstart, ystart = ystart, display = display
;read cat file
cat = read_file(infile)

;project coordinates at equator
radec = radec2deg(radec_cent[0], radec_cent[1], cat[0,*], cat[1,*])


 ;Define grid
ra_max   = max(radec[0,*], min = ra_min)
dec_max  = max(radec[1,*], min = dec_min)

ra_max  += 0.075D
dec_max += 0.075D
ra_min  -= 0.075D
dec_min -= 0.075D

    
delta_ra = ra_max-ra_min
delta_dec = dec_max-dec_min
    
rmin_log  = alog10(rmin)
rmax_log  = alog10(rmax)
bin_log   = (rmax_log-rmin_log)/nbin
r_arr     = 10^(rmin_log +bin_log*dindgen(nbin))

dr        = r_arr[1] - r_arr[0]



npix_ra  = long(delta_ra/dr) + 1L
npix_dec = long(delta_dec/dr) + 1L
  
;actual cell size
spix  = delta_ra/npix_ra

img = dblarr(npix_ra, npix_dec)

;compute pix coordinates for each object
xpix = long((radec[0,*] - ra_min)/spix)
ypix = long((radec[1,*] - dec_min)/spix)

;define a pseudo astrometry structure to handle masking
;coordinates of center in pixels
x_cent = long(-ra_min/spix)
y_cent = long(-dec_min/spix)
get_astr_st, spix, radec_cent[0], radec_cent[1], x_cent, y_cent, astr

;get cell numbers
cell_number = xpix+ ypix*npix_ra



;count_max = -1000L
ncells = npix_ra*npix_dec


;get the number of objects in each cell
histo_cells  = histogram(cell_number, max = ncells-1L, min = 0L)

;fill the image
img[lindgen(ncells)] = histo_cells
if(keyword_set(display)) then begin
   window, 1, xs = 500, ys = 500, retain = 2, title = 'input img'
   tv, bytscl(congrid(img, 500, 500))
endif
;apply mask
;non_masked_img = mask_cat_image(mask_file, img, radec_cent, ra_min, $
;                                dec_min, spix, xstart, ystart, mask, $
;                                astrometry =  astr)
mask_cat_image, mask_file, img, radec_cent, ra_min, dec_min, spix, xstart, $
                ystart, mask, mask_hdr, non_masked_img, astrometry =  astr
;print, 'in fill grid'
;help, mask
;help, mask_hdr
img_temp = dblarr(npix_ra, npix_dec)
img_temp[non_masked_img] = img[non_masked_img]
img = img_temp

;stop
END