;Follows the algorithm described in Szapudi et al 2005, ;Apj, 631 L1 to compute the angular correlation function of a dataset ;INPUT: ; infile: full path to ascii catalog or image. Fits files are assumed ; to be images only; WCS information need to be present in fits ; header. For catalogs, only the first two columns (supposedly ; ra, dec) are used. ; ; wtheta_file: output file for result. The output file contains ; 2 columns: theta and w(theta) ; ; mask_file: full path to ascii of fits mask. For fits masks, ; non masked pixels are assumed to be set to 0. Fits masks need ; to have WCS information in fits header. Ascii masks can ; contain region description for convexs (CONVEX), polygons ; (POLY) or circles (CIRCLE). The format for a mask file line is: ; - for convexes: CONVEX x1 y1 z1 c1 x2 y2 ..... xn yn zn cn ; with n arbitrary. (xi, yi, zi, ci) defines a halfspace. ; - for polygons: POLY ra1 dec1 ra2 dec2 .... ran decn where n ; arbitrary. The rai, deci are the vertices of the polygons ; - for circles: CIRCLE ra_c dec_c radius where ra_c, dec_c ; are the coordinates of the center. ; All coordinates (masks and input catalogs) are assumed to be ; decimal degrees (on sky). In case of ascii catalog, the ; coordinates will be projected on equator to use the small ; angles approximation (using the input of the radec_cent ; keyword). ; An ascii mask file NEEDS to contain at least one inclusion ; mask, that can be CONVEX, POLY or CIRCLE. The inclusion mask ; is specified by adding '_IN' to the mask region description ; ('CONVEX_IN', 'POLY_IN', 'CIRCLE_IN'). ; ; ;OUTPUT: ; A 2D array with output[0,*] = theta and output[1,*] = w_theta ; ;KEYWORDS: ; radec_cent: set this keyword to a 2 elements one-dimension ; array [ra_cent, dec_cent] with the coordinates of the center of ; the field corresponding to the ascii catalog. Needs to be set ; only in case of input ascii catalog, and not in case of input ; fits image. ; ; rmin: minimum value of theta for w(theta) computation, in ; degrees. Default: 0.001 degrees ; ; rmax: maximum value of theta for w(theta) computation, in ; degrees. Default: 0.5 degrees ; ; nbin: number of bins for w(theta) computation. Default: 10 ; ; threshold: additional threshold to be set on images (only for ; fits images). Mask will be updated accordingly. ; ; native_res: if keyword set, fits images are used with their ; native resolution (Default) ; ; custom_res: if set, fits image resolution will be reduced by a ; custom_res factor on each direction (ex. custom_res = 1.25 will ; reduced image x size and y size by 1.25, so a gain in virtual ; memory of 1.5625 compared to native resolution) ; ; xstart: set this keyword to a pixel number, in order to remove, ; on each side of the fits image, xstart pixels. ex. on a 2000 by ; 2000 image, setting xstart to 100 will internally use a 1800 by ; 2000 image, removing 100 pixels on the left and right edges. ; ; ystart: similar to xstart, but in the y direction ; ; display: if set, will display plots of images, random catalogs ; and results ; ; verbose: if set, will print a few details about computation. ; ; native_res, custom_res, xstart, and ystart apply only on fits image ; ; nbootstrap: if set no N, will generate N bootstrap samples, and ; compute errors on the correlation function ; ; njackknife: if set no N, will generate N jackniffe samples and ; get errors on correlation function from them ; ; ncut: if set to N, will cut the sample in N subsamples and get ; errors on correlation from dispersion ; ; cut_out: if set, will output the w(theta) computed for each cut ; (so generating ncut files) ;returns 1 if file is a fits file function check_files, file_name fits_open, file_name, fcb, /no_abort, message = message if(strmatch(message, '')) then begin ;fits file res = 1 fits_close, fcb endif else begin ;no fits file res = 0 endelse return, res end function wtheta_fft, infile, wtheta_file, mask_file, radec_cent = radec_cent, $ rmin = rmin, rmax = rmax, nbin = nbin,threshold = $ threshold, native_res = native_res, custom_res = $ custom_res, xstart = xstart, ystart = ystart, display = $ display, verbose = verbose, nbootstrap = nbootstrap, $ njackknife = njackknife, ncut = ncut, cut_out = cut_out if (not keyword_set(rmin)) then rmin = 0.001D if (not keyword_set(rmax)) then rmax = 0.5D if (not keyword_set(nbin)) then nbin = 10 if (not keyword_set(xstart)) then xstart = 0 if (not keyword_set(ystart)) then ystart = 0 time_in = systime(1) if (keyword_set(verbose)) then print, 'Put data on a grid' ;check whether input file is fits file check_fits = check_files(infile) if (check_fits) then begin if (keyword_set(verbose)) then $ print, 'Input file is considered as a fits image' fill_grid_image, infile, mask_file, rmin, rmax, nbin, img, $ ra_min, dec_min, ra_max, dec_max, spix, non_masked_img, proj_corner, $ threshold = threshold, native_res = native_res, xstart = xstart, $ ystart = ystart, custom_res = custom_res endif else begin if (keyword_set(verbose)) then $ print, 'Input file is considered as an ascii catalog' if(not keyword_set(radec_cent)) then $ stop, 'In case of ascii catalog, need radec_cent input !' 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 endelse nxpix = n_elements(img[*,0]) nypix = n_elements(img[0,*]) ;pixel coordinates of the center xpix_cent = -ra_min/spix ypix_cent = -dec_min/spix ;create random image if (keyword_set(verbose)) then print, 'Create random image' random = create_random(img, non_masked_img) ;get jackknife errors if needed if(n_elements(njackknife) ne 0) then begin ;stop get_wtheta_jackknife_errors, njackknife, img, random, non_masked_img, $ rmin, rmax, nbin, spix, wtheta_jackknife, wtheta_jackknife_err endif ;get errors from independent subsamples if needed if(n_elements(ncut) ne 0) then begin get_wtheta_cuts_errors, ncut, img, random, non_masked_img, $ rmin, rmax, nbin, spix, wtheta_cut, wtheta_cut_err endif ;Get mean mean_img = mean(img[non_masked_img]) ;remove mean and normalize by mean (= compute the overdensity) img[non_masked_img] -= mean_img img[non_masked_img] /= mean_img ;Add zero values around img_zero = zero_around(img) random_zero = zero_around(random) if(keyword_set(display)) then begin window, 0, xsize = 500, ysize = 500, retain = 2, title = 'img' tv, bytscl(congrid(img, 500, 500), 0, 10^(-4.5D)) window, 1, xsize = 500, ysize = 500, retain = 2, title = 'random' tv, bytscl(congrid(random, 500, 500)) endif ;stop ;compute ffts ... if (keyword_set(verbose)) then print, 'Compute FFTs' ft_img = fft(img_zero, /double) ft_random = fft(random_zero, /double) img_zero = 0 random_zero = 0 square_mod_img = (abs(ft_img))^2 ft_img = 0 square_mod_random = (abs(ft_random))^2 ft_random = 0 xi_img = real_part(fft(square_mod_img, /inverse, /double)) xi_random = real_part(fft(square_mod_random, /inverse, /double)) ;get xi(r) nxpix = n_elements(xi_img[*,0]) nypix = n_elements(xi_img[0,*]) xi_r = dblarr(nxpix, nypix) use = where(xi_random ne 0, count) xi_r[use] = xi_img[use]/xi_random[use] ;shift subscripts xi_r = shift(xi_r,nxpix/2 -1L, nypix/2 -1L) xi_img = shift(xi_img,nxpix/2 -1L, nypix/2 -1L) xi_random = shift(xi_random,nxpix/2 -1L, nypix/2 -1L) ;get position array mx = indgen(nxpix) - (nxpix/2 -1) rx = mx*spix my = indgen(nypix) - (nypix/2 -1) ry = my*spix ;Compute distances (in order to integrate xi(r) over r) dist = dblarr(nxpix, nypix) for iy = 0L, nypix -1L do begin dist[*, iy] = ((rx)^2 + (ry[iy])^2) endfor dist = sqrt(dist) ;Compute w(theta) ;Set up log binning rmin_log = alog10(rmin) rmax_log = alog10(rmax) bin_log = (rmax_log-rmin_log)/nbin r_arr = 10^(rmin_log +bin_log*dindgen(nbin +1L)) num = dblarr(nbin) den = dblarr(nbin) ;Loop over bins in r and integrate 'xi(r)' if (keyword_set(verbose)) then print, 'Integrate to get w(theta)' for ibin = 0, nbin-1 do begin search = where(dist GT r_arr[ibin] and dist LE r_arr[ibin +1]) dr = r_arr[ibin +1] - r_arr[ibin] dr2 = dist[search]*dr*2.*!dpi num[ibin] = total(xi_img[search]*dr2) den[ibin] = total(xi_random[search]*dr2) endfor wtheta = num/den ;get bootstrap error if needed if(n_elements(nbootstrap) ne 0) then begin seed = -150l print, 'youhou !!!!!!!!!!!!' help, mask get_wtheta_bootstrap_errors, nbootstrap, seed, dist, r_arr, $ xi_random, infile, radec_cent, mask_file, $ mask, mask_hdr, rmin, rmax, nbin, img, $ ra_min, $ dec_min, ra_max, dec_max, spix, $ non_masked_img, wtheta_err, $ xstart = xstart, ystart = ystart endif xi_img = 0 xi_random = 0 xi_r = 0 r_arr_out = 10^(rmin_log +bin_log/2.D +bin_log*dindgen(nbin)) if(keyword_set(display)) then begin ymin = 1.25*min(wtheta, max = ymax) ymax *= 1.25 window, 2, xsize = 500, ysize = 500, retain = 2 plot, [0.],[0.], xrange = [rmin, rmax], xs = 1, yrange = [10^(-5.), 1.], $ ys = 1,/nodata, /xlog, /ylog,charsize = 1.25, xtitle = '!7h!3',$ ytitle = '!3w(!7h)!3' vsym, 24, /fill oplot, r_arr_out, wtheta, psym = 8, symsize= 1. endif ;number of output columns n_col_out = 2 ;modify size of output array depending on keywords if(n_elements(nbootstrap) ne 0 ) then n_col_out++ if(n_elements(njackknife) ne 0) then n_col_out++ if(n_elements(ncut) ne 0) then n_col_out++ res = dblarr(n_col_out, nbin) res[0,*] = r_arr_out res[1,*] = wtheta col_out = 1 if(n_elements(nbootstrap) ne 0) then res[++col_out,*] = wtheta_err if(n_elements(njackknife) ne 0) then res[++col_out,*] = wtheta_jackknife_err if(n_elements(ncut) ne 0) then res[++col_out,*] = wtheta_cut_err openw, 44, wtheta_file printf, 44, res close, 44 ;output cuts if needed wtheta_file_cp = wtheta_file if(keyword_set(cut_out)) then begin ;get output files root wtheta_str_arr_pos = strsplit(wtheta_file_cp, '.', length = length) ;position of last period in namee ('.'); assumed to be the file extension nsubstring = n_elements(wtheta_str_arr_pos) extension = strmid(wtheta_file_cp, $ wtheta_str_arr_pos[nsubstring-1], $ length[nsubstring-1]) beginning = gettok(wtheta_file_cp, '.'+extension) for icut = 0, ncut-1 do begin wtheta_cut_file = beginning+'_cut'+strtrim(icut+1, 2)+'.'+extension out_arr = dblarr(2, nbin) out_arr[0,*] = r_arr_out out_arr[1,*] = wtheta_cut[icut, *] openw, 44, wtheta_cut_file printf, 44, out_arr close, 44 endfor endif time_out = systime(1) if (keyword_set(verbose)) then print, 'Done in '+$ strtrim(time_out-time_in, 2)+' s !' ;stop return, res END