;get jackknife errors for w(theta) given an input catalog (already
;loaded into an image img, and as well for the random "catalog"
;random). Create njackknife samples (by removing 1/njackknife for
;each jackknife sample) and compute w(theta) for each of those. The
;errors are then given by the ~ dispersion between the measures.
;N. B. input image should be the direct catalog, and NOT the
;overdensity !!!
pro get_wtheta_jackknife_errors, njackknife, img, random, non_masked_img, $
                                 rmin, rmax, nbin, spix, wtheta_j, wtheta_err

;create array for wtheta jackknife
wtheta_j = dblarr(njackknife, 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(njackknife))
nycut = long(njackknife*1./nxcut)


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

count = 0

;loop on jackknife 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 be removed
        xy_j = x_j+nxpix*y_j
        ;mask img and random
        img_j          = img
        img_j[xy_j]    = 0
        
        random_j       = random
        random_j[xy_j] = 0
        
                                ;get non masked pixels in jackknife
                                ;samples by removing the xy_j pixels
                                ;from the non_masked_list
        non_masked_j = setdifference(non_masked_img, xy_j)

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

;Compute the errors

wtheta_err = sqrt((njackknife-1))*cmapply('user:stddev', wtheta_j, 1)
;stop
end