diff --git a/fhd_core/HEALPix/healpix_cnv_generate.pro b/fhd_core/HEALPix/healpix_cnv_generate.pro index 08e18f09e..3dd8c7bbc 100644 --- a/fhd_core/HEALPix/healpix_cnv_generate.pro +++ b/fhd_core/HEALPix/healpix_cnv_generate.pro @@ -20,6 +20,7 @@ elements=obs.elements IF N_Elements(hpx_radius) EQ 0 THEN BEGIN IF Keyword_Set(mask) THEN BEGIN + xv_arr=meshgrid(dimension,elements,1) yv_arr=meshgrid(dimension,elements,2) ;set /ignore_refraction for speed, since we don't need to be exact @@ -32,15 +33,14 @@ ENDIF ELSE radius=hpx_radius ;all angles in DEGREES ;uses RING index scheme - ;check if a string, if it is assume it is a filepath to a save file with the desired indices ; (will NOT be over-written with the indices) IF Keyword_Set(restrict_hpx_inds) AND (size(restrict_hpx_inds,/type) NE 7) THEN restrict_hpx_inds=observation_healpix_inds_select(obs) IF size(restrict_hpx_inds,/type) EQ 7 THEN BEGIN file_path_use=restrict_hpx_inds - IF file_test(file_path_use) EQ 0 THEN file_path_use=filepath(file_path_use,root=Rootdir('fhd'),subdir='Observations') + IF file_test(file_path_use) EQ 0 THEN file_path_use=filepath(file_path_use,root=Rootdir('fhd'),subdir='Observations') - IF file_test(file_path_use) THEN BEGIN + IF file_test(file_path_use) THEN BEGIN hpx_inds=getvar_savefile(file_path_use,'hpx_inds') nside_test=getvar_savefile(file_path_use,names=sav_contents) IF Max(strmatch(StrLowCase(sav_contents),'nside')) EQ 1 THEN nside=getvar_savefile(file_path_use,'nside') ELSE BEGIN @@ -52,6 +52,7 @@ ENDIF IF ~Keyword_Set(nside) THEN BEGIN pix_sky=4.*!Pi*!RaDeg^2./Product(Abs(obs.astr.cdelt)) Nside=2.^(Ceil(ALOG(Sqrt(pix_sky/12.))/ALOG(2))) ;=1024. for 0.1119 degrees/pixel + ENDIF npix=nside2npix(nside) @@ -60,9 +61,18 @@ IF Keyword_Set(divide_pixel_area) THEN BEGIN ENDIF ELSE pixel_area_cnv=1. ;turn this off for now IF N_Elements(hpx_inds) GT 1 THEN BEGIN - pix2vec_ring,nside,hpx_inds,pix_coords - vec2ang,pix_coords,pix_dec,pix_ra,/astro - apply_astrometry, obs, ra_arr=pix_ra, dec_arr=pix_dec, x_arr=xv_hpx, y_arr=yv_hpx, /ad2xy + ;set up orthoslant grid + x_vals = repmat(findgen(obs.dimension), 1, obs.dimension) + y_vals = Transpose(x_vals) + ;shape them correctly + x_vals = reform(x_vals, obs.dimension^2) + y_vals = reform(y_vals, obs.dimension^2) + ;transform to ra/dec, back to orthoslant. For different observations obs would be different + ;between the two + apply_astrometry, obs, x_arr=x_vals, y_arr=y_vals, ra_arr=ra_vals, dec_arr=dec_vals, /xy2ad + apply_astrometry, obs, ra_arr=ra_vals, dec_arr=dec_vals, x_arr=xv_hpx, y_arr=yv_hpx, /ad2xy + + ENDIF ELSE BEGIN ang2vec,obs.obsdec,obs.obsra,cen_coords,/astro Query_disc,nside,cen_coords,radius,hpx_inds0,ninds,/deg @@ -90,27 +100,29 @@ ENDIF ELSE BEGIN ENDELSE ; Test for pixels past the horizon. We don't need to be precise with this, so turn off precession, etc.. -Eq2Hor,pix_ra, pix_dec, obs.JD0, alt_arr, az_arr, nutate=0,precess=0,aberration=0, refract=0, lon=obs.lon, alt=obs.alt, lat=obs.lat -horizon_i = where(alt_arr LE 0, n_horizon, complement=h_use) -IF n_horizon GT 0 THEN BEGIN - print,String(format='("Cutting ",A, " HEALPix pixels that were below the horizon.")',Strn(n_horizon)) - xv_hpx = xv_hpx[h_use] - yv_hpx = yv_hpx[h_use] - hpx_inds = hpx_inds[h_use] -ENDIF +;Eq2Hor,pix_ra, pix_dec, obs.JD0, alt_arr, az_arr, nutate=0,precess=0,aberration=0, refract=0, lon=obs.lon, alt=obs.alt, lat=obs.lat +;horizon_i = where(alt_arr LE 0, n_horizon, complement=h_use) +;IF n_horizon GT 0 THEN BEGIN +; print, "im here" +; print,String(format='("Cutting ",A, " HEALPix pixels that were below the horizon.")',Strn(n_horizon)) +; xv_hpx = xv_hpx[h_use] +; yv_hpx = yv_hpx[h_use] +; hpx_inds = hpx_inds[h_use] +;ENDIF x_frac=1.-(xv_hpx-Floor(xv_hpx)) y_frac=1.-(yv_hpx-Floor(yv_hpx)) min_bin=Min(Floor(xv_hpx)+dimension*Floor(yv_hpx))>0L -max_bin=Max(Ceil(xv_hpx)+dimension*Ceil(yv_hpx))<(dimension*elements-1L) +max_bin=Max(Ceil(xv_hpx)+dimension*Ceil(yv_hpx))<(dimension*elements-1L) +; give hist some data, min and max for bins, binsize=1? h00=histogram(Floor(xv_hpx)+dimension*Floor(yv_hpx),min=min_bin,max=max_bin,/binsize,reverse_ind=ri00) h01=histogram(Floor(xv_hpx)+dimension*Ceil(yv_hpx),min=min_bin,max=max_bin,/binsize,reverse_ind=ri01) h10=histogram(Ceil(xv_hpx)+dimension*Floor(yv_hpx),min=min_bin,max=max_bin,/binsize,reverse_ind=ri10) h11=histogram(Ceil(xv_hpx)+dimension*Ceil(yv_hpx),min=min_bin,max=max_bin,/binsize,reverse_ind=ri11) + htot=h00+h01+h10+h11 inds=where(htot,n_img_use) - n_arr=htot[inds] i_use=inds+min_bin @@ -151,7 +163,9 @@ FOR i=0L,n_img_use-1L DO BEGIN *ija[i]=ija0 ENDFOR - +;sa=weighting function? +;set hpx inds to be the right thing for the orthoslant +hpx_inds = findgen(obs.dimension^2) hpx_cnv={nside:nside,ija:ija,sa:sa,i_use:i_use,inds:hpx_inds} IF tag_exist(obs,'healpix') THEN BEGIN IF N_Elements(restrict_hpx_inds) NE 1 THEN ind_list="UNSPECIFIED" ELSE ind_list=restrict_hpx_inds diff --git a/fhd_core/HEALPix/healpix_snapshot_cube_generate.pro b/fhd_core/HEALPix/healpix_snapshot_cube_generate.pro index ad9c1b3c6..4cd755f20 100644 --- a/fhd_core/HEALPix/healpix_snapshot_cube_generate.pro +++ b/fhd_core/HEALPix/healpix_snapshot_cube_generate.pro @@ -3,7 +3,7 @@ PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,v ps_kbinsize=ps_kbinsize,ps_kspan=ps_kspan,ps_beam_threshold=ps_beam_threshold,ps_nfreq_avg=ps_nfreq_avg,$ rephase_weights=rephase_weights,n_avg=n_avg,vis_weights=vis_weights,split_ps_export=split_ps_export,$ restrict_hpx_inds=restrict_hpx_inds,hpx_radius=hpx_radius,cmd_args=cmd_args,save_uvf=save_uvf,save_imagecube=save_imagecube,$ - obs_out=obs_out,psf_out=psf_out,ps_tile_flag_list=ps_tile_flag_list,_Extra=extra + obs_out=obs_out,psf_out=psf_out,ps_tile_flag_list=ps_tile_flag_list,beam_ptr=beam_ptr,_Extra=extra t0=Systime(1) @@ -184,31 +184,61 @@ PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,v undefine_fhd,weights_arr1,variance_arr1,residual_arr1,dirty_arr1,model_arr1 ;free memory for beam_arr later! IF iter EQ n_iter-1 THEN undefine_fhd,beam_arr + ;set up the pointers for later arrays + + dirty_img_arr=Ptrarr(n_pol,n_freq,/allocate) + weights_img_arr=Ptrarr(n_pol,n_freq,/allocate) + variance_img_arr=Ptrarr(n_pol,n_freq,/allocate) + model_img_arr=Ptrarr(n_pol,n_freq,/allocate) + + dirty_uv_arr=Ptrarr(n_pol,n_freq,/allocate) + weights_uv_arr=Ptrarr(n_pol,n_freq,/allocate) + variance_uv_arr=Ptrarr(n_pol,n_freq,/allocate) + model_uv_arr=Ptrarr(n_pol,n_freq,/allocate) + FOR pol_i=0,n_pol-1 DO BEGIN IF dirty_flag THEN BEGIN dirty_cube=fltarr(n_hpx,n_freq_use) ;write index in much more efficient memory access order - FOR fi=Long64(0),n_freq_use-1 DO dirty_cube[n_hpx*fi]=Temporary(*dirty_hpx_arr[pol_i,fi]) + ;remove temporary to avoid losing arrays + FOR fi=Long64(0),n_freq_use-1 DO dirty_cube[n_hpx*fi]=*dirty_hpx_arr[pol_i,fi] ENDIF IF model_flag THEN BEGIN model_cube=fltarr(n_hpx,n_freq_use) - FOR fi=Long64(0),n_freq_use-1 DO model_cube[n_hpx*fi]=Temporary(*model_hpx_arr[pol_i,fi]) + FOR fi=Long64(0),n_freq_use-1 DO model_cube[n_hpx*fi]=*model_hpx_arr[pol_i,fi] ENDIF IF residual_flag THEN BEGIN res_cube=fltarr(n_hpx,n_freq_use) - FOR fi=Long64(0),n_freq_use-1 DO res_cube[n_hpx*fi]=Temporary(*residual_hpx_arr[pol_i,fi]) + FOR fi=Long64(0),n_freq_use-1 DO res_cube[n_hpx*fi]=*residual_hpx_arr[pol_i,fi] ENDIF weights_cube=fltarr(n_hpx,n_freq_use) - FOR fi=Long64(0),n_freq_use-1 DO weights_cube[n_hpx*fi]=Temporary(*weights_hpx_arr[pol_i,fi]) + FOR fi=Long64(0),n_freq_use-1 DO weights_cube[n_hpx*fi]=*weights_hpx_arr[pol_i,fi] variance_cube=fltarr(n_hpx,n_freq_use) - FOR fi=Long64(0),n_freq_use-1 DO variance_cube[n_hpx*fi]=Temporary(*variance_hpx_arr[pol_i,fi]) + FOR fi=Long64(0),n_freq_use-1 DO variance_cube[n_hpx*fi]=*variance_hpx_arr[pol_i,fi] beam_squared_cube=fltarr(n_hpx,n_freq_use) - FOR fi=Long64(0),n_freq_use-1 DO beam_squared_cube[n_hpx*fi]=Temporary(*beam_hpx_arr[pol_i,fi]) + FOR fi=Long64(0),n_freq_use-1 DO beam_squared_cube[n_hpx*fi]=*beam_hpx_arr[pol_i,fi] + + ;now do the FFTs back and normalize + + FOR fi=Long64(0),n_freq_use-1 DO *dirty_img_arr[pol_i,fi]=reform(*dirty_hpx_arr[pol_i,fi], SQRT(n_hpx), SQRT(n_hpx)) + FOR fi=Long64(0),n_freq_use-1 DO *dirty_uv_arr[pol_i,fi]=fft_shift(fft(fft_shift((*dirty_img_arr[pol_i,fi])*(obs_out.degpix*!DtoR)^2), /inverse)) + + + FOR fi=Long64(0),n_freq_use-1 DO *weights_img_arr[pol_i,fi]=reform(*weights_hpx_arr[pol_i,fi], SQRT(n_hpx), SQRT(n_hpx)) + FOR fi=Long64(0),n_freq_use-1 DO *weights_uv_arr[pol_i,fi]=fft_shift(fft(fft_shift((*weights_img_arr[pol_i,fi])*(obs_out.degpix*!DtoR)^2), /inverse)) + + + FOR fi=Long64(0),n_freq_use-1 DO *variance_img_arr[pol_i,fi]=reform(*variance_hpx_arr[pol_i,fi], SQRT(n_hpx), SQRT(n_hpx)) + FOR fi=Long64(0),n_freq_use-1 DO *variance_uv_arr[pol_i,fi]=fft_shift(fft(fft_shift((*variance_img_arr[pol_i,fi])*(obs_out.degpix*!DtoR)^2), /inverse)) + + FOR fi=Long64(0),n_freq_use-1 DO *model_img_arr[pol_i,fi]=reform(*model_hpx_arr[pol_i,fi], SQRT(n_hpx), SQRT(n_hpx)) + FOR fi=Long64(0),n_freq_use-1 DO *model_uv_arr[pol_i,fi]=fft_shift(fft(fft_shift((*model_img_arr[pol_i,fi])*(obs_out.degpix*!DtoR)^2), /inverse)) + ;call fhd_save_io first to obtain the correct path. Will NOT update status structure yet fhd_save_io,status_str,file_path_fhd=file_path_fhd,var=cube_name[iter],pol_i=pol_i,path_use=path_use,/no_save,_Extra=extra @@ -219,9 +249,21 @@ PRO healpix_snapshot_cube_generate,obs_in,status_str,psf_in,cal,params,vis_arr,v fhd_save_io,status_str,file_path_fhd=file_path_fhd,var=cube_name[iter],pol_i=pol_i,/force,_Extra=extra dirty_cube=(model_cube=(res_cube=(weights_cube=(variance_cube=(beam_squared_cube=0))))) ENDFOR + + image_filepath = file_path_fhd+ '_' + uvf_name[iter] +'_interped_images.sav' + save, filename = image_filepath, dirty_img_arr, weights_img_arr, variance_img_arr, model_img_arr, obs_out, /compress + + uvf_filepath = file_path_fhd+ '_' + uvf_name[iter] +'_gridded_uvf.sav' + save, filename = uvf_filepath, dirty_uv_arr, weights_uv_arr, variance_uv_arr, model_uv_arr, obs_out, /compress + undefine_fhd,dirty_hpx_arr,model_hpx_arr,weights_hpx_arr,variance_hpx_arr,beam_hpx_arr,residual_hpx_arr ENDFOR + + + obs_out=obs ;for return + + Ptr_free,vis_weights_use timing=Systime(1)-t0 IF ~Keyword_Set(silent) THEN print,'HEALPix cube export timing: ',timing,t_split,t_hpx diff --git a/fhd_utils/FFT/fft_shift_cube.pro b/fhd_utils/FFT/fft_shift_cube.pro new file mode 100644 index 000000000..5f5e012d9 --- /dev/null +++ b/fhd_utils/FFT/fft_shift_cube.pro @@ -0,0 +1,9 @@ +FUNCTION fft_shift_cube, image_cube +dimension=(size(image_cube,/dimension))[0] +elements=(size(image_cube,/dimension))[1] +n_slice=(size(image_cube,/dimension))[2] +shifted_image_cube = Fltarr(dimension, elements, n_slice) +FOR s=0,n_slice-1 DO $ + shifted_image_cube[*, *, s] = fft_shift(Reform(image_cube[*, *, s])) +RETURN, shifted_image_cube +END \ No newline at end of file diff --git a/idl_startup.pro b/idl_startup.pro new file mode 100644 index 000000000..e69de29bb