;Follows the algorithm described in Szapudi et al 2005, ;Apj, 631 L1 to compute the angular cross correlation function of two ;datasets ; As it is, works for image X image, and image X ascii catalog; it ; doesn't work on the case ascii catalog X ascii catalog (yet) ;INPUT: ; file1: 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. ; ; file2: same as file1 ; ; There is no preferred order for fits image or ascii catalog ; for file1 or file2. ; ; wtheta_file: output file for result. The output file contains ; 2 columns: theta and w(theta) ; ; mask_file1: 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'). ; ; mask_file2: same as mask_file1 ; ; ;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. ;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 ;get image or cat filename given input files ;if two images, cat is set to empty string ''; if one image and one ;cat, image2 is set to '' pro get_image_cat_files, file_name1, file_name2, image_name1, image_name2, $ cat_name ;check whether input files are images or cat (assumes for now that ;fits are images and all other catalogs) check_file1 = check_files(file_name1) check_file2 = check_files(file_name2) if(check_file1 eq 0 and check_file2 eq 0) then begin stop, 'One of the input file should be an image ...' endif if(check_file1 eq 1 and check_file2 eq 1) then begin image_name1 = file_name1 image_name2 = file_name2 cat_name = '' endif if(check_file1 eq 1 and check_file2 eq 0) then begin image_name1 = file_name1 image_name2 = '' cat_name = file_name2 endif if(check_file1 eq 0 and check_file2 eq 1) then begin image_name1 = file_name2 image_name2 = '' cat_name = file_name1 endif end function xwtheta_fft, file1, file2, wtheta_file, mask_file1, mask_file2, $ 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 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(custom_res)) then custom_res = 1. time_in = systime(1) ;get who's who (cat file or image file) get_image_cat_files, file1, file2, image_name1, image_name2, $ cat_name if(keyword_set(verbose)) then begin print, image_name1+' is considered as a fits image' if(strmatch(image_name2, '')) then begin print, cat_name+' is considered as an ascii catalog' if(not keyword_set(radec_cent)) then $ stop, 'In case of ascii catalog, need radec_cent input !' endif else begin print, image_name2+' is considered as a fits image' endelse endif if(keyword_set(verbose)) then print, 'Put data on grid' ;get first image fill_grid_image, image_name1, mask_file1, rmin, rmax, nbin, img1, $ ra_min, dec_min, ra_max, dec_max, spix, non_masked_img1, proj_corner, $ threshold = threshold, native_res = native_res, xstart = xstart, $ ystart = ystart, custom_res = custom_res, ra_cent = radec_cent[0], $ dec_cent = radec_cent[1] nxpix = n_elements(img1[*,0]) nypix = n_elements(img1[0,*]) ;get second image (from cat or image) if(strmatch(image_name2, '')) then begin fill_grid_cat_xwtheta, cat_name, radec_cent, mask_file2, rmin, rmax, $ nbin, img2, nxpix, nypix, proj_corner[0], proj_corner[1], ra_max, $ dec_max, spix, xstart, ystart, custom_res, non_masked_img2 endif else begin fill_grid_image, image_name2, mask_file2, rmin, rmax, nbin, img2, $ ra_min, dec_min, ra_max, dec_max, spix, non_masked_img2, proj_corner, $ threshold = threshold, native_res = native_res, xstart = xstart, $ ystart = ystart, custom_res = custom_res, ra_cent = radec_cent[0], $ dec_cent = radec_cent[1] endelse if(keyword_set(display)) then begin window, 0, xsize = 500, ysize = 500, retain = 2, title = 'img1' tv, bytscl(congrid(img1, 500, 500), 0, 10^(-4.)) window, 1, xsize = 500, ysize = 500, retain = 2, title = 'img2' tv, bytscl(congrid(img2, 500, 500), 0, 10^(-4.)) endif ;pixel coordinates of the center xpix_cent = -ra_min/spix ypix_cent = -dec_min/spix ;image means mean_img1 = mean(img1[non_masked_img1]) mean_img2 = mean(img2[non_masked_img2]) ;create random images if (keyword_set(verbose)) then print, 'Create random images' random_img1 = create_random(img1, non_masked_img1) random_img2 = create_random(img2, non_masked_img2) if(keyword_set(display)) then begin window, 2, xsize = 500, ysize = 500, retain = 2, title = 'random img1' tv, bytscl(congrid(random_img1, 500, 500), 0, 10^(-4.)) window, 3, xsize = 500, ysize = 500, retain = 2, title = 'random img2' tv, bytscl(congrid(random_img2, 500, 500), 0, 10^(-4.)) endif ;remove means img1[non_masked_img1] -= mean_img1 img2[non_masked_img2] -= mean_img2 ; get overdensity img1[non_masked_img1] /= mean_img1 img2[non_masked_img2] /= mean_img2 non_masked_img1 = -1 non_masked_img2 = -1 ;Add zero values around img1_zero = zero_around(img1) img2_zero = zero_around(img2) img1 = 0 img2 = 0 random_img1_zero = zero_around(random_img1) random_img2_zero = zero_around(random_img2) random_img1 = 0 random_img2 = 0 ;compute ffts ... if (keyword_set(verbose)) then print, 'Compute FFTs' ft_img1 = fft(img1_zero, /double) img1_zero = 0 ft_img2 = fft(img2_zero, /double) img2_zero = 0 ft_random_img1 = fft(random_img1_zero, /double) random_img1_zero = 0 ft_random_img2 = fft(random_img2_zero, /double) random_img2_zero = 0 d_r_prod = abs(ft_img1)*abs(ft_img2) ft_img1 = 0 ft_img2 = 0 r_prod = abs(ft_random_img1)*abs(ft_random_img2) ft_random_img1 = 0 ft_random_img2 = 0 xi_img = real_part(fft(d_r_prod, /inverse, /double)) dr_prod = 0 xi_random = real_part(fft(r_prod, /inverse, /double)) r_prod = 0 ;shift subscripts nxpix = n_elements(xi_img[*,0]) nypix = n_elements(xi_img[0,*]) 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 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) ;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)) wtheta = dblarr(nbin) den = dblarr(nbin) num = dblarr(nbin) 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 res = dblarr(2, nbin) r_arr_out = 10^(rmin_log +bin_log/2.D +bin_log*dindgen(nbin)) res[0,*] = r_arr_out res[1,*] = wtheta 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 openw, 44, wtheta_file printf, 44, res close, 44 time_out = systime(1) if (keyword_set(verbose)) then print, 'Done in '+$ strtrim(time_out-time_in, 2)+' s !' return, res END