LOFAR beam

fhd_core/beam_modeling/beam_image_hyperresolved.pro

@@ -13,28 +13,46 @@ function beam_image_hyperresolved,antenna1,antenna2,ant_pol1,ant_pol2,freq_i,$
   ;; Defined only at non-zero pixels for memory reduction
   pix_use=*antenna1[0].pix_use
   psf_image_dim = antenna1.psf_image_dim
-
-  ;; Jones pointers and antenna beams
-  Jones1=antenna1.Jones[*,*,freq_i]
-  Jones2=antenna2.Jones[*,*,freq_i]
-  beam_ant1=DComplex(*(antenna1.response[ant_pol1,freq_i]))
-  beam_ant2=DComplex(Conj(*(antenna2.response[ant_pol2,freq_i])))
-
-  ;; Amplitude of the response from ant1 is Sqrt(|J1[0,pol1]|^2 + |J1[1,pol1]|^2)
-  ;; Amplitude of the response from ant2 is Sqrt(|J2[0,pol2]|^2 + |J2[1,pol2]|^2)
-  ;; Amplitude of the baseline response is the product of the antenna responses
-  power_zenith_beam=Sqrt((abs(*Jones1[0,ant_pol1])^2.+abs(*Jones1[1,ant_pol1])^2.)*$
-    (abs(*Jones2[0,ant_pol2])^2.+abs(*Jones2[1,ant_pol2])^2.))
-
   ;; Create only one full-scale array for memory reduction
   image_power_beam = Dcomplexarr(psf_image_dim,psf_image_dim)
 
-  ;; Co-opt this array to calculate the power at zenith
-  image_power_beam[pix_use] = power_zenith_beam
-  power_zenith=Interpolate(image_power_beam,zen_int_x,zen_int_y,cubic=-0.5)
-
-  ;; Redefine the array to be the image power beam, normalized to zenith
-  image_power_beam[pix_use] = (power_zenith_beam*beam_ant1*beam_ant2)/power_zenith
+
+  if ~tag_exist(antenna1,'F_mag') and ~tag_exist(antenna1,'F_phase') and $
+    ~tag_exist(antenna2,'F_mag') and ~tag_exist(antenna2,'F_phase') then begin
+    ;; Jones pointers and antenna beams
+    Jones1=antenna1.Jones[*,*,freq_i]
+    Jones2=antenna2.Jones[*,*,freq_i]
+    beam_ant1=DComplex(*(antenna1.response[ant_pol1,freq_i]))
+    beam_ant2=DComplex(Conj(*(antenna2.response[ant_pol2,freq_i])))
+
+    ;; Amplitude of the response from ant1 is Sqrt(|J1[0,pol1]|^2 + |J1[1,pol1]|^2)
+    ;; Amplitude of the response from ant2 is Sqrt(|J2[0,pol2]|^2 + |J2[1,pol2]|^2)
+    ;; Amplitude of the baseline response is the product of the antenna responses
+    power_zenith_beam=Sqrt((abs(*Jones1[0,ant_pol1])^2.+abs(*Jones1[1,ant_pol1])^2.)*$
+      (abs(*Jones2[0,ant_pol2])^2.+abs(*Jones2[1,ant_pol2])^2.))
+
+    ;; Co-opt this array to calculate the power at zenith
+    image_power_beam[pix_use] = power_zenith_beam
+    power_zenith=Interpolate(image_power_beam,zen_int_x,zen_int_y,cubic=-0.5)
+
+    ;; Redefine the array to be the image power beam, normalized to zenith
+    image_power_beam[pix_use] = (power_zenith_beam*beam_ant1*beam_ant2)/power_zenith
+  endif else begin
+    ; F_mag is the diagonal amplitude response of the antenna to unpolarized radiation, 
+    ; ordered as [n_pol, n_freq]. 
+    ; F_phase is the diagonal time delay response of the antenna 
+    ; to unpolarized radiation, ordered as [n_pol, n_freq].
+    power_zenith_beam = (*antenna1.F_mag[ant_pol1,freq_i]) * (*antenna2.F_mag[ant_pol2,freq_i])
+
+    ;; Co-opt this array to calculate the power at zenith
+    image_power_beam[pix_use] = power_zenith_beam
+    power_zenith=Interpolate(image_power_beam,zen_int_x,zen_int_y,cubic=-0.5)
+
+    ;; Redefine the array to be the image power beam, normalized to zenith
+    image_power_beam[pix_use] = power_zenith_beam * $
+      (*antenna1.F_phase[ant_pol1,freq_i]) * Conj(*antenna2.F_phase[ant_pol2,freq_i]) / power_zenith
+
+  endelse
 
   return, image_power_beam
 

fhd_core/beam_modeling/beam_setup.pro

@@ -236,7 +236,11 @@ FUNCTION beam_setup,obs,status_str,antenna,file_path_fhd=file_path_fhd,restore_l
         bi_use2=Reform(rebin((ant_1_arr+1),ant_1_n,ant_2_n)+$
           Rebin(Transpose(ant_2_arr+1),ant_1_n,ant_2_n)*baseline_mod,baseline_group_n)
         bi_use = [bi_use, bi_use2]
-        IF Max(bi_use) GT bi_max THEN bi_use=bi_use[where(bi_use LE bi_max)]
+        IF Max(bi_use) GT bi_max THEN BEGIN
+          ;Check that the group is not larger than the maximum baseline index and that the baselines exist
+          inds = where(bi_use LE bi_max,n_count)
+          if n_count GT 0 then bi_use=bi_use[inds] else continue
+        ENDIF
         bi_use_i=where(bi_hist0[bi_use],n_use)
         IF n_use GT 0 THEN bi_use=bi_use[bi_use_i]
         baseline_group_n=N_Elements(bi_use)

fhd_core/beam_modeling/fhd_struct_init_antenna.pro

@@ -98,10 +98,22 @@ FOR fi=0L,nfreq_bin-1 DO BEGIN
 ENDFOR
 
 ;initialize antenna structure
-antenna_str={n_pol:n_ant_pol,antenna_type:instrument,names:ant_names,model_version:beam_model_version,freq:freq_center,nfreq_bin:nfreq_bin,$
-    n_ant_elements:0,Jones:Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin),coupling:Ptrarr(n_ant_pol,nfreq_bin),gain:Ptrarr(n_ant_pol),coords:Ptrarr(3),$
-    delays:Ptr_new(),size_meters:0.,height:0.,response:Ptrarr(n_ant_pol,nfreq_bin),group_id:Lonarr(n_ant_pol)-1,pix_window:Ptr_new(),pix_use:Ptr_new(),$
-    psf_image_dim:0.,psf_scale:0.}
+; If the Jones matrix is going to be decomposed into separate magnitude, delay, and projection contributions, then add those tags.
+; Keep a combined Jones matrix for backwards compatibility
+If instrument EQ 'lofar' then begin
+  antenna_str={n_pol:n_ant_pol,antenna_type:instrument,names:ant_names,model_version:beam_model_version,$
+      freq:freq_center,nfreq_bin:nfreq_bin,n_ant_elements:0,Jones:Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin),$
+      coupling:Ptrarr(n_ant_pol,nfreq_bin),gain:Ptrarr(n_ant_pol),coords:Ptrarr(3),delays:Ptr_new(),$
+      size_meters:0.,height:0.,rotation:0,response:Ptrarr(n_ant_pol,nfreq_bin),group_id:Lonarr(n_ant_pol)-1,$
+      pix_window:Ptr_new(),pix_use:Ptr_new(),psf_image_dim:0.,psf_scale:0.,K_proj: Ptrarr(n_ant_pol, n_ant_pol, nfreq_bin), $
+      F_phase: Ptrarr(n_ant_pol, nfreq_bin), F_mag: Ptrarr(n_ant_pol, nfreq_bin)}
+endif else begin
+  antenna_str={n_pol:n_ant_pol,antenna_type:instrument,names:ant_names,model_version:beam_model_version,$
+      freq:freq_center,nfreq_bin:nfreq_bin,n_ant_elements:0,Jones:Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin),$
+      coupling:Ptrarr(n_ant_pol,nfreq_bin),gain:Ptrarr(n_ant_pol),coords:Ptrarr(3),delays:Ptr_new(),$
+      size_meters:0.,height:0.,rotation:0,response:Ptrarr(n_ant_pol,nfreq_bin),group_id:Lonarr(n_ant_pol)-1,$
+      pix_window:Ptr_new(),pix_use:Ptr_new(),psf_image_dim:0.,psf_scale:0.}
+endelse
 
 ;update structure with instrument-specific values, and return as a structure array, with an entry for each tile/antenna
 ;first, update to include basic configuration data
@@ -228,8 +240,10 @@ if N_elements(instrument) GT 1 then begin
   endfor  
 endif else antenna=Call_function(tile_gain_fn,obs,antenna,za_arr=za_arr,az_arr=az_arr,psf_image_dim=psf_image_dim,Jdate_use=Jdate_use,_Extra=extra) ;mwa_beam_setup_gain
 
-;Finally, update antenna structure to include the response of each antenna
-antenna=general_antenna_response(obs,antenna,za_arr=za_arr,az_arr=az_arr,psf_image_dim=psf_image_dim,_Extra=extra)
+if ~tag_exist(antenna, 'K_proj') and ~tag_exist(antenna, 'F_mag') and ~tag_exist(antenna, 'F_phase') then begin
+  ;Finally, update antenna structure to include the response of each antenna
+  antenna=general_antenna_response(obs,antenna,za_arr=za_arr,az_arr=az_arr,psf_image_dim=psf_image_dim,_Extra=extra)
+endif
 
 timing=Systime(1)-t0
 RETURN,antenna

fhd_core/polarization/rotate_jones_matrix.pro

@@ -8,7 +8,11 @@ psf_image_dim = antenna.psf_image_dim
 for pol_i = 0, antenna.n_pol-1 do begin
     for pol_j=0, antenna.n_pol-1 do begin
         Jones[pol_i,pol_j] = Ptr_new(DComplexarr(psf_image_dim,psf_image_dim))
-        (*Jones[pol_i,pol_j])[pix_use] = (*antenna.Jones[pol_i,pol_j])
+        if ~tag_exist(antenna1,'K_proj') then begin
+            (*Jones[pol_i,pol_j])[pix_use] = (*antenna.Jones[pol_i,pol_j])
+        endif else begin
+            (*Jones[pol_i,pol_j])[pix_use] = (*antenna.K_proj[pol_i,pol_j])
+        endelse
     endfor
 endfor
 

fhd_core/setup_metadata/fhd_struct_init_meta.pro

@@ -118,6 +118,7 @@ ENDIF ELSE BEGIN
     IF missing_n GT 0 THEN tile_flag0[missing_i]=1
     tile_flag=Ptrarr(n_pol) & FOR pol_i=0,n_pol-1 DO tile_flag[pol_i]=Ptr_new(tile_flag0)
     date_obs=hdr.date_obs
+    base_gain=0
 
     IF ~Keyword_Set(time_offset) THEN time_offset=0d
     time_offset/=(24.*3600.)
@@ -173,10 +174,10 @@ projection_slant_orthographic,astr=astr,degpix=degpix2,obsra=obsra,obsdec=obsdec
     epoch=2000.,JDate=JD0
 
 Eq2Hor,obsra,obsdec,JD0,obsalt,obsaz,lat=lat,lon=lon,alt=Mean(alt)
-meta={obsra:Float(obsra),obsdec:Float(obsdec),zenra:Float(zenra),zendec:Float(zendec),phasera:Float(phasera),phasedec:Float(phasedec),$
+meta={obsra:Double(obsra),obsdec:Double(obsdec),zenra:Double(zenra),zendec:Double(zendec),phasera:Double(phasera),phasedec:Double(phasedec),$
     epoch:Float(epoch),tile_names:tile_names,lon:Float(lon),lat:Float(lat),alt:Float(alt),JD0:Double(JD0),Jdate:Double(Jdate),$
     astr:astr,obsx:Float(obsx),obsy:Float(obsy),zenx:Float(zenx),zeny:Float(zeny),obsaz:Float(obsaz),obsalt:Float(obsalt),$
-    delays:beamformer_delays,tile_height:Float(tile_height),tile_flag:tile_flag,orig_phasera:Float(orig_phasera),orig_phasedec:Float(orig_phasedec),$
+    delays:beamformer_delays,tile_height:Float(tile_height),tile_flag:tile_flag,orig_phasera:Double(orig_phasera),orig_phasedec:Double(orig_phasedec),$
     time_res:Float(time_res),base_gain:base_gain}
 
 RETURN,meta

fhd_core/source_modeling/vis_model_transfer.pro

@@ -1,29 +1,39 @@
 function vis_model_transfer,obs,params,model_transfer
 
-  data_nbaselines = obs.nbaselines
+  data_uniq_time_inds = [0,uniq(params.time)+1]
+  ;data_uniq_time_inds = data_uniq_time_inds[0:n_elements(data_uniq_time_inds)-2]
   data_n_time = obs.n_time
 
-  ;; Option to transfer pre-made and unflagged model visbilities
-  vis_model_arr=PTRARR(obs.n_pol,/allocate)
+  if strlowcase(strmid(model_transfer, strlen(model_transfer)-7, 7)) eq '.uvfits' then begin
+    ; model_transfer is a UVFITS file
+    if file_test(model_transfer) then begin
+      uvfits_read,hdr,params_model,layout,vis_model_arr,vis_weights,file_path_vis=model_transfer,n_pol=obs.n_pol,silent=silent,error=error
+    endif else message, model_transfer + ' not found during model transfer.'
 
-  for pol_i=0, obs.n_pol-1 do begin
-    transfer_name = model_transfer + '/' + obs.obsname + '_vis_model_'+obs.pol_names[pol_i]+'.sav'
-    if ~file_test(transfer_name) then $
-      message, transfer_name + ' not found during model transfer.'
-    vis_model_arr[pol_i] = getvar_savefile(transfer_name,'vis_model_ptr')
-    print, "Model visibilities transferred from " + transfer_name
-  endfor
-
-  ;; Get the params file associated with the model visibilities
-  if file_test(model_transfer+'/'+obs.obsname+'_params.sav') then begin
-    params_model = getvar_savefile(model_transfer+'/'+obs.obsname+'_params.sav','params')
   endif else begin
-    if file_test(file_dirname(model_transfer)+'/metadata/'+obs.obsname+'_params.sav') then begin
-      params_model = getvar_savefile(file_dirname(model_transfer)+'/metadata/'+obs.obsname+'_params.sav','params')
-    endif else message, 'No params file found in model transfer directory.'
-  endelse 
-  model_n_time = n_elements(uniq(params_model.time))
-  model_nbaselines = (uniq(params_model.time))[0]+1
+    ; model_transfer is a directory of saved visibility model files
+    vis_model_arr=PTRARR(obs.n_pol,/allocate)
+
+    for pol_i=0, obs.n_pol-1 do begin
+      transfer_name = model_transfer + '/' + obs.obsname + '_vis_model_'+obs.pol_names[pol_i]+'.sav'
+      if ~file_test(transfer_name) then $
+        message, transfer_name + ' not found during model transfer.'
+      vis_model_arr[pol_i] = getvar_savefile(transfer_name,'vis_model_ptr')
+      print, "Model visibilities transferred from " + transfer_name
+    endfor
+
+    ;; Get the params file associated with the model visibilities
+    if file_test(model_transfer+'/'+obs.obsname+'_params.sav') then begin
+      params_model = getvar_savefile(model_transfer+'/'+obs.obsname+'_params.sav','params')
+    endif else begin
+      if file_test(file_dirname(model_transfer)+'/metadata/'+obs.obsname+'_params.sav') then begin
+        params_model = getvar_savefile(file_dirname(model_transfer)+'/metadata/'+obs.obsname+'_params.sav','params')
+      endif else message, 'No params file found in model transfer directory.'
+    endelse 
+  endelse
+
+  model_uniq_time_inds = [0,uniq(params_model.time)+1]
+  model_n_time = n_elements(model_uniq_time_inds)
 
   if model_n_time NE data_n_time then begin
   ;; Exclude flagged times from the model visibilities if the number of time steps do not match
@@ -93,19 +103,19 @@ function vis_model_transfer,obs,params,model_transfer
   ;Initialize matched model pointer array
   matched_model = PTRARR(obs.n_pol,/allocate)
   for pol_i=0, obs.n_pol-1 do begin
-    if size(*vis_model_arr[pol_i], /type) eq 6 then *matched_model[pol_i] = complex(FLTARR(obs.n_freq, data_nbaselines * data_n_time))
-    if size(*vis_model_arr[pol_i], /type) eq 9 then *matched_model[pol_i] = dcomplex(DBLARR(obs.n_freq, data_nbaselines * data_n_time))
+    if size(*vis_model_arr[pol_i], /type) eq 6 then *matched_model[pol_i] = complex(FLTARR(obs.n_freq, obs.nbaselines * data_n_time))
+    if size(*vis_model_arr[pol_i], /type) eq 9 then *matched_model[pol_i] = dcomplex(DBLARR(obs.n_freq, obs.nbaselines * data_n_time))
   endfor
 
   ; In each matched timestep, match the baselines and fill a new, matched model array 
   ; Unmatched times in the transferred model are not included in the matched model array 
   for time_i=0, data_n_time-1 do begin
     ;suba is the subset of indices in the first array that are also in the second.
     ;subb is the subset of indices in the second array that are also in the first.
-    match, data_baseline_index[time_i*data_nbaselines:(time_i+1)*data_nbaselines-1], $
-      model_baseline_index[matched_times[time_i]*model_nbaselines:(matched_times[time_i]+1)*model_nbaselines-1], $
+    match, data_baseline_index[data_uniq_time_inds[time_i]:data_uniq_time_inds[time_i+1]-1], $
+      model_baseline_index[model_uniq_time_inds[matched_times[time_i]]:model_uniq_time_inds[matched_times[time_i]+1]-1], $
       suba, subb
-    for pol_i=0, obs.n_pol-1 do (*matched_model[pol_i])[*,suba+time_i*(data_nbaselines)] = (*vis_model_arr[pol_i])[*,subb+matched_times[time_i]*model_nbaselines]
+    for pol_i=0, obs.n_pol-1 do (*matched_model[pol_i])[*,suba+data_uniq_time_inds[time_i]] = (*vis_model_arr[pol_i])[*,subb+model_uniq_time_inds[matched_times[time_i]]]
   endfor
 
   for pol_i=0, obs.n_pol-1 do vis_model_arr[pol_i] = Pointer_copy(matched_model[pol_i])

fhd_core/visibility_manipulation/vis_flag_basic.pro

@@ -144,6 +144,27 @@ ENDFOR
 (*obs.baseline_info).time_use = time_use
 ;;; End time-based flagging
 
+;;; Baseline-based flagging
+IF instrument EQ 'lofar' THEN BEGIN
+    ; Intra-cabinet baselines need to be flagged due to cabinet RFI 
+    ; Pairs of sequential tiles are in the same cabinet,
+    ; i.e. tiles 0 and 1, tiles 2 and 3, etc. for just the core stations (1 through 48)
+    tile_A_i=(*obs.baseline_info).tile_A-1
+    tile_B_i=(*obs.baseline_info).tile_B-1
+    n_core_stations = 48
+
+    FOR pairs_i=0, n_core_stations/2-1 DO BEGIN
+        inds_tile_a = where(tile_A_i EQ (pairs_i*2),n_count)
+        IF n_count GT 0 THEN BEGIN 
+            inds_tile_b = where(tile_B_i[inds_tile_a] EQ (pairs_i*2+1),n_count)
+        ENDIF
+        IF n_count GT 0 THEN BEGIN
+            FOR pol_i=0,n_pol-1 DO (*vis_weight_ptr[pol_i])[*,inds_tile_a[inds_tile_b]]=0
+        ENDIF
+    ENDFOR
+ENDIF
+;;; End baseline-based flagging
+
 IF Keyword_Set(unflag_all) THEN BEGIN
     tile_use[*]=1
     freq_use[*]=1

fhd_core/visibility_manipulation/vis_header_extract.pro

@@ -68,16 +68,18 @@ IF found_baseline NE 1 THEN BEGIN
         found_baseline = 1 ; 'BASELINE' only needs to be present if the antennas are not included separately
     ENDIF ELSE print,"WARNING: Group parameter BASELINE not found within uvfits header PTYPE keywords" 
 ENDIF
-uu_i=(where(Strmatch(param_list,'UU', /fold_case),found_uu))[0] & IF found_uu NE 1 THEN print,"WARNING: Group parameter UU not found within uvfits header PTYPE keywords"
-vv_i=(where(Strmatch(param_list,'VV', /fold_case),found_vv))[0] & IF found_vv NE 1 THEN print,"WARNING: Group parameter VV not found within uvfits header PTYPE keywords"
-ww_i=(where(Strmatch(param_list,'WW', /fold_case),found_ww))[0] & IF found_ww NE 1 THEN print,"WARNING: Group parameter WW not found within uvfits header PTYPE keywords"
+
+uu_i=(where(Strmatch(param_list,'UU', /fold_case),found_uu))[0] & IF found_uu LT 1 THEN print,"WARNING: Group parameter UU not found within uvfits header PTYPE keywords"
+vv_i=(where(Strmatch(param_list,'VV', /fold_case),found_vv))[0] & IF found_vv LT 1 THEN print,"WARNING: Group parameter VV not found within uvfits header PTYPE keywords"
+ww_i=(where(Strmatch(param_list,'WW', /fold_case),found_ww))[0] & IF found_ww LT 1 THEN print,"WARNING: Group parameter WW not found within uvfits header PTYPE keywords"
 date_i=where(Strmatch(param_list,'DATE', /fold_case),found_date) & IF found_date LT 1 THEN print,"WARNING: Group parameter DATE not found within uvfits header PTYPE keywords"
-IF (found_baseline NE 1) OR (found_uu NE 1) OR (found_vv NE 1) OR (found_ww NE 1) OR (found_date LT 1) THEN error=1 
+
+IF (found_baseline NE 1) OR (found_uu LT 1) OR (found_vv LT 1) OR (found_ww LT 1) OR (found_date LT 1) THEN error=1 
 
 IF found_date GT 1 THEN BEGIN
-    Jdate0=Double(sxpar(header,String(format='("PZERO",I1)',date_i[0]+1)))
+    Jdate0=Double(sxpar(header,String(format='("PZERO",I0)',date_i[0]+1)))
 ENDIF ELSE BEGIN
-    Jdate0=Double(sxpar(header,String(format='("PZERO",I1)',date_i+1)))
+    Jdate0=Double(sxpar(header,String(format='("PZERO",I0)',date_i+1)))
 ENDELSE
 ;Jdate_extract=sxpar(header,String(format='("PZERO",I1)',date_i+1),count=found_jd0)
 ;IF Keyword_Set(found_jd0) THEN IF Jdate_extract GT 2.4E6 THEN Jdate0=Jdate_extract