;Cut input sample into ncut independent subsamples; compute wtheta and
;derive errors on wtheta from the standard deviation of these results.

;N. B. input image should be the direct catalog, and NOT the
;overdensity !!!

pro get_wtheta_cuts_errors, ncuts, img, random, non_masked_img, $
                            rmin, rmax, nbin, spix, wtheta_j, wtheta_err
;create array for wtheta jackknife
wtheta_j = dblarr(ncuts, nbin)


;dimension of input image
nxpix = n_elements(img[*,0])
nypix = n_elements(img[0,*])

;get number of cuts along x and y directions
nxcut = long(sqrt(ncuts))
nycut = long(ncuts*1./nxcut)


;dimension of subimage to be considered at each iteration
nxpix_j = nxpix*1./nxcut
nypix_j = nypix*1./nycut

count = 0
;loop on subsamples
for ix = 0, nxcut-1 do begin
    ;x subscripts to be removed 
    x_j = ix*nxpix_j + lindgen(nxpix_j)
    ;put in matrix form
    x_j = cmreplicate(x_j, [nypix_j])
    for iy = 0, nycut-1 do begin
        ;y subscripts to be removed
        y_j = iy*nypix_j + lindgen(nypix_j)
    
        ;put in matrix form
        ;y_j = transpose(cmreplicate(reverse(y_j), [nxpix_j]))
        y_j = transpose(cmreplicate(y_j, [nxpix_j]))
        
                
        ;get one dimensional subscripts to keep
        xy_j = x_j+nxpix*y_j
        ;mask img and random
        img_j          = img[xy_j]
        
        random_j       = random[xy_j]
        
                                ;get non masked pixels in subsample
                                ;samples by intersecting with the xy_j
                                ;pixels from the non_masked_list
        non_masked_temp = setintersection(non_masked_img, xy_j)
                                ;get non masked pixels "good"
                                ;subscripts (given the dimension of
                                ;the subsample subimage)
        x_non_masked_temp = non_masked_temp mod nxpix
        y_non_masked_temp = long(non_masked_temp/nxpix)

        ;convert pixels to cut "frame"
        x_non_masked_j = x_non_masked_temp-ix*nxpix_j
        y_non_masked_j = y_non_masked_temp-iy*nypix_j
        
        ;get one dimensional subscripts
        non_masked_j = x_non_masked_j + nxpix_j*y_non_masked_j

        ;if(iy eq 0) then stop
        ;compute wtheta
        wtheta_j[ix + nxcut*iy, *] = $
          get_wtheta_img(img_j, random_j, non_masked_j, rmin, rmax, nbin, spix)
         ; window, count, xs = 500, ys = 500, retain = 2
;          tvscl, bytscl(congrid(random_j, 500, 500))

         count++
    endfor
endfor


;Compute the errors

wtheta_err = cmapply('user:stddev', wtheta_j, 1)


;stop
end