Gaussian source modeling

dictionary.md

@@ -188,6 +188,10 @@ This is a work in progress; please add keywords as you find them in alphabetical
   -*Turn off/on*: 0/1 <br />
   -*Default*: undefined (off) <br />
 
+**gaussian_source_models**:  uses SHAPE information provided in the sky catalog to build Gaussian models of extended sources. See Line et al. 2020 for more details on implementation. The models are only accurate to within ~10\%, and this is an ongoing issue (see Issue [\#211](https://github.com/EoRImaging/FHD/issues/211)). <br />
+  -*Turn off/on*: 0/1 <br />
+  -*Default*: undefined (off) <br />
+
 **include_catalog_sources**: Adds sources in the file specified by catalog_file_path to the source_list for simulations (used in vis_simulate.pro) <br />
    -*Default* : 0 <br />
    -*If set to zero, and no source_array is passed to vis_simulate, then no point sources will be included in the simulation. Originally, sim_versions_wrapper would load the source array directly and pass it along to vis_simulate, then additional sources could be included via catalog_file_path. This feature has been turned off, so this parameter alone specified whether or not point sources will be included in simulations.* !Q <br />

fhd_core/calibration/generate_source_cal_list.pro

@@ -47,7 +47,7 @@ IF N_Elements(beam_threshold) EQ 0 THEN beam_threshold=0.05
 no_restrict_sources_default = 1
 
 IF Keyword_Set(model_visibilities) THEN BEGIN
-    IF N_Elements(model_flux_threshold) GT 0 THEN flux_threshold=model_flux_threshold ELSE flux_threshold=0.
+    IF N_Elements(model_flux_threshold) GT 0 THEN flux_threshold=model_flux_threshold
     IF N_Elements(allow_sidelobe_model_sources) GT 0 THEN allow_sidelobe_sources=allow_sidelobe_model_sources ELSE allow_sidelobe_sources=0
     IF Keyword_Set(allow_sidelobe_sources) THEN IF N_Elements(no_restrict_model_sources) EQ 0 THEN no_restrict_sources=1 ELSE no_restrict_sources=no_restrict_model_sources $
         ELSE no_restrict_sources=no_restrict_sources_default
@@ -56,7 +56,7 @@ IF Keyword_Set(model_visibilities) THEN BEGIN
     IF Keyword_Set(allow_sidelobe_sources) THEN beam_threshold=0.01
     IF N_Elements(beam_model_threshold) GT 0 THEN beam_threshold=beam_model_threshold
 ENDIF ELSE BEGIN
-    IF N_Elements(calibration_flux_threshold) GT 0 THEN flux_threshold=calibration_flux_threshold ELSE flux_threshold=0.
+    IF N_Elements(calibration_flux_threshold) GT 0 THEN flux_threshold=calibration_flux_threshold
     IF N_Elements(allow_sidelobe_cal_sources) GT 0 THEN allow_sidelobe_sources=allow_sidelobe_cal_sources ELSE allow_sidelobe_sources=0
     IF Keyword_Set(allow_sidelobe_sources) THEN IF N_Elements(no_restrict_cal_sources) EQ 0 THEN no_restrict_sources=1 ELSE no_restrict_sources=no_restrict_cal_sources $
         ELSE no_restrict_sources=no_restrict_sources_default
@@ -116,8 +116,13 @@ ENDIF
 IF n_use GT 0 THEN BEGIN
     catalog=catalog[i_use]
     IF N_Elements(spectral_index) GT 1 THEN spectral_index=spectral_index[i_use]
-    source_list=source_comp_init(n_sources=n_use,freq=freq_use,ra=catalog.ra,dec=catalog.dec,$
+    if tag_exist(catalog, 'shape') then begin
+      source_list=source_comp_init(n_sources=n_use,freq=freq_use,ra=catalog.ra,dec=catalog.dec,$
+        alpha=spectral_index,extend=catalog.extend,shape=catalog.shape)
+    endif else begin
+      source_list=source_comp_init(n_sources=n_use,freq=freq_use,ra=catalog.ra,dec=catalog.dec,$
         alpha=spectral_index,extend=catalog.extend)
+    endelse
 
     apply_astrometry, obs, ra_arr=source_list.ra, dec_arr=source_list.dec, x_arr=x_arr, y_arr=y_arr, /ad2xy    
     source_list.x=x_arr
@@ -132,33 +137,31 @@ IF n_use GT 0 THEN BEGIN
 
     ;If flux_threshold is negative, assume that it is an UPPER bound, and only include the fainter sources
     flux_I_use = source_list.flux.I
-    IF flux_threshold LT 0 THEN flux_I_use = -flux_I_use
+    IF N_elements(flux_threshold) NE 0 THEN BEGIN
+        IF flux_threshold LT 0 THEN flux_I_use = -flux_I_use
+        src_use=where((x_arr GE fft_alias_range) AND (x_arr LE dimension-1-fft_alias_range) AND (y_arr GE fft_alias_range) $
+            AND (y_arr LE elements-1-fft_alias_range) AND (flux_I_use GT flux_threshold) AND (flux_I_use NE 0),n_src_use)
+    ENDIF ELSE BEGIN
+        src_use=where((x_arr GE fft_alias_range) AND (x_arr LE dimension-1-fft_alias_range) AND (y_arr GE fft_alias_range) $
+            AND (y_arr LE elements-1-fft_alias_range) AND (flux_I_use NE 0),n_src_use)
+    ENDELSE    
 
-    src_use=where((x_arr GE fft_alias_range) AND (x_arr LE dimension-1-fft_alias_range) AND (y_arr GE fft_alias_range) $
-        AND (y_arr LE elements-1-fft_alias_range) AND (flux_I_use GT flux_threshold) AND (flux_I_use NE 0),n_src_use)
-
     IF n_src_use EQ 0 THEN RETURN,source_comp_init(n_sources=0,freq=obs.freq_center);
     src_use2=where(beam_mask[Round(x_arr[src_use]),Round(y_arr[src_use])],n_src_use)
     IF n_src_use GT 0 THEN src_use=src_use[src_use2]
     source_list=source_list[src_use]
-    beam_list=Ptrarr(n_pol<2)
-
     influence=source_list.flux.I*beam[source_list.x,source_list.y]
-    
+
     order=Reverse(sort(influence))
     source_list=source_list[order]
     source_list.id=Lindgen(n_src_use)
     IF Keyword_Set(no_extend) THEN source_list.extend=Ptrarr(n_src_use) ELSE BEGIN
         extend_i=where(Ptr_valid(source_list.extend),n_extend)
         FOR ext_i=0L,n_extend-1 DO BEGIN
-            ex_spectral_index=source_list[extend_i[ext_i]].alpha
             extend_list=*source_list[extend_i[ext_i]].extend
             apply_astrometry, obs, ra_arr=extend_list.ra, dec_arr=extend_list.dec, x_arr=x_arr, y_arr=y_arr, /ad2xy
             extend_list.x=x_arr
             extend_list.y=y_arr
-            extend_list.alpha=ex_spectral_index
-            FOR i=0,7 DO extend_list.flux.(i)=extend_list.flux.(i)*(freq_use/catalog[extend_i[ext_i]].freq)^ex_spectral_index
-;            FOR pol_i=0,(n_pol<2)-1 DO extend_list.flux.(pol_i)=extend_list.flux.I*(*beam_list[pol_i])[extend_i[ext_i]]
             *source_list[extend_i[ext_i]].extend=extend_list
         ENDFOR
     ENDELSE

fhd_core/deconvolution/source_detection/source_comp_init.pro

@@ -1,5 +1,6 @@
 FUNCTION source_comp_init,source_comp_in,n_sources=n_sources,xvals=xvals,yvals=yvals,frequency=frequency,$
-    ra=ra,dec=dec,flux=flux,id=id,StoN=StoN,alpha=alpha,extend=extend,gain_factor=gain_factor,overwrite=overwrite,_Extra=extra
+    ra=ra,dec=dec,flux=flux,id=id,StoN=StoN,alpha=alpha,extend=extend,gain_factor=gain_factor,shape=shape,$
+    overwrite=overwrite,_Extra=extra
 ; NOTE: It is very important that the tags of the source structure be able to autocomplete to the keywords of this function
 ; For example, the current xvals is fine, because the .x tag can autocomplete to xvals. If a new x_error (for example)
 ; keyword or tag were introduced, this would break load_source_catalog.pro
@@ -9,9 +10,11 @@ IF N_Elements(n_sources) EQ 0 THEN $
 
 flux_struct={flux,xx:0.,yy:0.,xy:Complex(0.),yx:Complex(0.),I:0.,Q:0.,U:0.,V:0.}
 beam_struct={beam,xx:0.,yy:0.,xy:Complex(0.),yx:Complex(0.),I:0.,Q:0.,U:0.,V:0.}
+shape_struct={x:0.,y:0.,angle:0.} ;gets used if keyword gaussian_source_models is set
 ;flux order is 0-3: xx, yy, xy, yx in apparent brightness; 4-7: I, Q, U, V in sky brightness
 ;flag type codes are 0: no flag, 1: low confidence 2: sidelobe contamination
-struct_base={id:-1L,x:0.,y:0.,ra:0.,dec:0.,ston:0.,freq:100.,alpha:0.,gain:1.,flag:0,extend:Ptr_new(),flux:flux_struct, beam:beam_struct}
+struct_base={id:-1L,x:0.,y:0.,ra:0.,dec:0.,ston:0.,freq:100.,alpha:0.,gain:1.,flag:0,extend:Ptr_new(),$
+    flux:flux_struct,shape:shape_struct,beam:beam_struct}
 source_comp_new=Replicate(struct_base,n_sources>1)
 
 IF Keyword_Set(xvals) THEN source_comp_new.x=xvals
@@ -29,7 +32,7 @@ IF Keyword_Set(alpha) THEN source_comp_new.alpha=alpha ;spectral index
 IF Keyword_Set(extend) THEN source_comp_new.extend=extend ;extended source component list (not an extended source if Ptr_valid(source_comp[si].extend)=0)
 IF Keyword_Set(frequency) THEN IF Mean(frequency) GT 1E5 THEN source_comp_new.freq=frequency/1E6 ELSE source_comp_new.freq=frequency ;frequency in MHz
 IF Keyword_Set(gain_factor) THEN source_comp_new.gain=gain_factor
-
+IF Keyword_Set(shape) THEN source_comp_new.shape=shape ;gaussian model parameters
 IF Keyword_Set(source_comp_in) AND ~Keyword_Set(overwrite) THEN source_comp=[source_comp_in,source_comp_new] ELSE source_comp=source_comp_new
 RETURN,source_comp
 END  

fhd_core/fhd_main.pro

@@ -1,5 +1,5 @@
 PRO fhd_main, file_path_vis, status_str, export_images=export_images, cleanup=cleanup, recalculate_all=recalculate_all,$
-    mapfn_recalculate=mapfn_recalculate, grid_recalculate=grid_recalculate,$
+    mapfn_recalculate=mapfn_recalculate, grid_recalculate=grid_recalculate,beam_recalculate=beam_recalculate,$
     n_pol=n_pol, flag_visibilities=flag_visibilities, silent=silent, deconvolve=deconvolve, transfer_mapfn=transfer_mapfn,$
     tile_flag_list=tile_flag_list, diffuse_calibrate=diffuse_calibrate, diffuse_model=diffuse_model,$
     file_path_fhd=file_path_fhd, force_data=force_data, force_no_data=force_no_data, freq_start=freq_start, freq_end=freq_end,$
@@ -29,7 +29,7 @@ log_filepath=file_path_fhd+'_log.txt'
 IF Keyword_Set(!Journal) THEN journal
 
 data_flag=fhd_setup(file_path_vis,status_str,export_images=export_images,cleanup=cleanup,recalculate_all=recalculate_all,$
-    mapfn_recalculate=mapfn_recalculate,grid_recalculate=grid_recalculate,$
+    mapfn_recalculate=mapfn_recalculate,grid_recalculate=grid_recalculate,beam_recalculate=beam_recalculate,$
     n_pol=n_pol,flag_visibilities=flag_visibilities,deconvolve=deconvolve,transfer_mapfn=transfer_mapfn,$
     transfer_weights=transfer_weights,file_path_fhd=file_path_fhd,force_data=force_data,force_no_data=force_no_data,$
     calibrate_visibilities=calibrate_visibilities,transfer_calibration=transfer_calibration,$
@@ -74,10 +74,15 @@ IF data_flag LE 0 THEN BEGIN
     IF Keyword_Set(deproject_w_term) THEN vis_arr=simple_deproject_w_term(obs,params,vis_arr,direction=deproject_w_term)
 
     ;Read in or construct a new beam model. Also sets up the structure PSF
-    print,'Calculating beam model'
-    psf=beam_setup(obs,status_str,antenna,file_path_fhd=file_path_fhd,restore_last=0,silent=silent,timing=t_beam,no_save=no_save,_Extra=extra)
-    IF Keyword_Set(t_beam) THEN IF ~Keyword_Set(silent) THEN print,'Beam modeling time: ',t_beam
-    jones=fhd_struct_init_jones(obs,status_str,file_path_fhd=file_path_fhd,restore=0,mask=beam_mask,_Extra=extra)
+    IF keyword_set(beam_recalculate) THEN BEGIN
+        print,'Calculating beam model'
+        psf=beam_setup(obs,status_str,antenna,file_path_fhd=file_path_fhd,restore_last=0,silent=silent,timing=t_beam,no_save=no_save,_Extra=extra)
+        IF Keyword_Set(t_beam) THEN IF ~Keyword_Set(silent) THEN print,'Beam modeling time: ',t_beam
+        jones=fhd_struct_init_jones(obs,status_str,file_path_fhd=file_path_fhd,restore=0,mask=beam_mask,_Extra=extra)
+    ENDIF ELSE BEGIN
+        fhd_save_io,status_str,psf,file_path_fhd=file_path_fhd,var_name='psf',/restore
+        fhd_save_io,status_str,jones,file_path_fhd=file_path_fhd,var_name='jones',/restore
+    ENDELSE
 
     IF Keyword_Set(transfer_weights) THEN BEGIN
         flag_visibilities=0 ;
@@ -149,7 +154,9 @@ IF data_flag LE 0 THEN BEGIN
         vis_weights_update,vis_weights,obs,psf,params,_Extra=extra
         if keyword_set(cal_stop) then begin
           fhd_save_io,status_str,obs,var='obs',/compress,file_path_fhd=file_path_fhd,_Extra=extra ;need beam_integral for PS
-          message, "cal_stop initiated"
+          print, "cal_stop initiated. Returning."
+          run_report, start_mem, t0, silent=silent
+          RETURN
         endif
     ENDIF
 
@@ -272,12 +279,7 @@ ENDIF
 undefine_fhd,map_fn_arr,image_uv_arr,weights_arr,model_uv_arr,vis_arr,vis_weights,vis_model_arr
 undefine_fhd,obs,cal,jones,layout,psf,antenna,fhd_params,skymodel,skymodel_cal,skymodel_update
 
-end_mem = (MEMORY(/HIGHWATER) - start_mem)/1e9
-timing=Systime(1)-t0
-IF ~Keyword_Set(silent) THEN BEGIN
-  print,'Memory required (GB): ',Strn(Round(end_mem))
-  print,'Full pipeline time (minutes): ',Strn(Round(timing/60.))
-ENDIF
+run_report, start_mem, t0, silent=silent
 print,''
 !except=except
 

fhd_core/gridding/visibility_grid.pro

@@ -200,7 +200,7 @@ bin_n=histogram(xmin+ymin*dimension,binsize=1,reverse_indices=ri,min=0) ;should
 bin_i=where(bin_n,n_bin_use)
 
 ind_ref=indgen(max(bin_n))
-n_vis=Float(Total(bin_n))
+n_vis=(Total(double(bin_n)))
 FOR fi=0L,n_f_use-1 DO n_vis_arr[fi_use[fi]]=Total(Long(xmin[fi,*] GT 0))
 obs.nf_vis=n_vis_arr
 

fhd_core/setup_metadata/fhd_setup.pro

@@ -1,5 +1,5 @@
 FUNCTION fhd_setup,file_path_vis,status_str,export_images=export_images,cleanup=cleanup,recalculate_all=recalculate_all,$
-    mapfn_recalculate=mapfn_recalculate,grid_recalculate=grid_recalculate,$
+    mapfn_recalculate=mapfn_recalculate,grid_recalculate=grid_recalculate,beam_recalculate=beam_recalculate,$
     n_pol=n_pol,flag_visibilities=flag_visibilities,deconvolve=deconvolve,transfer_mapfn=transfer_mapfn,$
     transfer_weights=transfer_weights,healpix_recalculate=healpix_recalculate,$
     file_path_fhd=file_path_fhd,force_data=force_data,force_no_data=force_no_data,$
@@ -31,8 +31,8 @@ IF N_Elements(flag_visibilities) EQ 0 THEN flag_visibilities=0
 IF N_Elements(transfer_mapfn) EQ 0 THEN transfer_mapfn=0
 IF N_Elements(save_visibilities) EQ 0 THEN save_visibilities=1
 IF N_Elements(healpix_recalculate) EQ 0 THEN healpix_recalculate=recalculate_all
-;IF N_Elements(beam_recalculate) EQ 0 THEN IF status_str.psf EQ 0 THEN beam_recalculate=1 $
-;    ELSE beam_recalculate=recalculate_all
+IF N_Elements(beam_recalculate) EQ 0 THEN IF status_str.psf EQ 0 THEN beam_recalculate=1 $
+    ELSE beam_recalculate=recalculate_all
 
 IF Keyword_Set(n_pol) THEN n_pol1=n_pol ELSE BEGIN
     IF status_str.obs GT 0 THEN BEGIN
@@ -75,11 +75,11 @@ ENDIF ELSE IF N_Elements(mapfn_recalculate) EQ 0 THEN mapfn_recalculate=0
 IF mapfn_recalculate GT 0 THEN grid_recalculate=1
 IF grid_recalculate GT 0 THEN data_flag=0
 
-;IF Keyword_Set(beam_recalculate) THEN BEGIN
-;    status_str.psf=0
-;    status_str.antenna=0
-;    status_str.jones=0
-;ENDIF
+IF Keyword_Set(beam_recalculate) THEN BEGIN
+    status_str.psf=0
+    status_str.antenna=0
+    status_str.jones=0
+ENDIF
 status_str.hpx_cnv=0
 IF Keyword_Set(mapfn_recalculate) THEN grid_recalculate=1
 IF Keyword_Set(grid_recalculate) THEN BEGIN

fhd_core/setup_metadata/load_source_catalog.pro

@@ -1,13 +1,19 @@
-FUNCTION load_source_catalog, catalog_filepath, varname=varname
+FUNCTION load_source_catalog, catalog, varname=varname
     ; Restore a source catalog from a .sav file, being careful about changed structure definitions
 
     IF N_Elements(varname) EQ 0 THEN varname="catalog"
 
     ;define structure BEFORE restoring, in case the definition has changed
     cat_init=source_comp_init(n_sources=0)
-    source_cat = getvar_savefile(catalog_filepath, varname, /compatibility_mode)
+    IF size(catalog,/type) EQ 8 THEN BEGIN
+        source_cat = catalog ;If a valid structure is supplied, use that 
+        catalog_filepath = ''
+    ENDIF ELSE BEGIN
+        catalog_filepath=catalog
+        source_cat = getvar_savefile(catalog_filepath, varname, /compatibility_mode)
+    ENDELSE
 
-    extend_i=where(Ptr_valid(source_cat.extend),n_ext,complement=point_i,ncomp=n_point)
+    extend_i=where(Ptr_valid(source_cat.extend),n_ext)
     ignore_tags = ['EXTEND']
 
     ; If the source list structure definition has changed we have to update the structure ourselves
@@ -42,16 +48,21 @@ FUNCTION load_source_catalog, catalog_filepath, varname=varname
     ENDFOR
     IF n_missing + n_extra GT 0 THEN BEGIN
         IF n_missing GT 0 THEN BEGIN
-            print,"WARNING! Mis-matched source list structure definition for file", catalog_filepath
+            print,"WARNING! Mis-matched source list structure definition for file ", catalog_filepath
             print,"Missing tags replaced with default values:", missing_tags
         ENDIF
 
         IF n_extra GT 0 THEN BEGIN
-            print,"WARNING! Mis-matched source list structure definition for file", catalog_filepath
+            print,"WARNING! Mis-matched source list structure definition for file ", catalog_filepath
             print,"Extraneous tags have been removed:", extra_tags
         ENDIF
+
         source_cat_mod = fill_source_list(source_cat, good_tag_names, good_tag_ids)
-        FOR ext_i=0,n_ext-1 DO source_cat_mod[extend_i[ext_i]].extend = Ptr_new(fill_source_list(*source_cat[extend_i[ext_i]].extend, good_tag_names, good_tag_ids))
+        FOR ext_i=0,n_ext-1 DO BEGIN
+	    extend_struct = load_source_catalog(*source_cat[extend_i[ext_i]].extend)
+            source_cat_mod[extend_i[ext_i]].extend = Ptr_new(extend_struct)
+        ENDFOR 
+
         RETURN, source_cat_mod
     ENDIF ELSE RETURN, source_cat