diff --git a/fhd_core/beam_modeling/beam_image_hyperresolved.pro b/fhd_core/beam_modeling/beam_image_hyperresolved.pro index 14a22501..b7eb1304 100644 --- a/fhd_core/beam_modeling/beam_image_hyperresolved.pro +++ b/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 diff --git a/fhd_core/beam_modeling/beam_setup.pro b/fhd_core/beam_modeling/beam_setup.pro index c09aa4f3..20157db6 100644 --- a/fhd_core/beam_modeling/beam_setup.pro +++ b/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) diff --git a/fhd_core/beam_modeling/fhd_struct_init_antenna.pro b/fhd_core/beam_modeling/fhd_struct_init_antenna.pro index 19d5bfa2..a4739cb9 100644 --- a/fhd_core/beam_modeling/fhd_struct_init_antenna.pro +++ b/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 diff --git a/fhd_core/polarization/rotate_jones_matrix.pro b/fhd_core/polarization/rotate_jones_matrix.pro index 80d2dfba..6a75f6c1 100644 --- a/fhd_core/polarization/rotate_jones_matrix.pro +++ b/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 diff --git a/fhd_core/setup_metadata/fhd_struct_init_meta.pro b/fhd_core/setup_metadata/fhd_struct_init_meta.pro index 7e124137..084b31d7 100644 --- a/fhd_core/setup_metadata/fhd_struct_init_meta.pro +++ b/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 diff --git a/fhd_core/source_modeling/vis_model_transfer.pro b/fhd_core/source_modeling/vis_model_transfer.pro index d590d2c0..4ce23fb0 100644 --- a/fhd_core/source_modeling/vis_model_transfer.pro +++ b/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,8 +103,8 @@ 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 @@ -102,10 +112,10 @@ function vis_model_transfer,obs,params,model_transfer 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]) diff --git a/fhd_core/visibility_manipulation/vis_flag_basic.pro b/fhd_core/visibility_manipulation/vis_flag_basic.pro index 7708ae98..497db605 100644 --- a/fhd_core/visibility_manipulation/vis_flag_basic.pro +++ b/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 diff --git a/fhd_core/visibility_manipulation/vis_header_extract.pro b/fhd_core/visibility_manipulation/vis_header_extract.pro index 94c9fe48..e4adb1b5 100644 --- a/fhd_core/visibility_manipulation/vis_header_extract.pro +++ b/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 diff --git a/instrument_config/lofar_beam_setup_gain.pro b/instrument_config/lofar_beam_setup_gain.pro new file mode 100644 index 00000000..d589ffba --- /dev/null +++ b/instrument_config/lofar_beam_setup_gain.pro @@ -0,0 +1,225 @@ +FUNCTION rot_cubic_complex, img, angle_deg, x_interp, y_interp, center_ind + + theta = -angle_deg * !dpi / 180.0 + cos_t = COS(theta) + sin_t = SIN(theta) + + x_rot = (x_interp-center_ind)*cos_t - (y_interp-center_ind)*sin_t + center_ind + y_rot = (x_interp-center_ind)*sin_t + (y_interp-center_ind)*cos_t + center_ind + + ; Interpolate all at once + out = INTERPOLATE(img, x_rot, y_rot, CUBIC=-0.5) + RETURN, out +END + +FUNCTION lofar_beam_setup_gain,obs,antenna,file_path_fhd=file_path_fhd,$ + za_arr=za_arr,az_arr=az_arr,psf_image_dim=psf_image_dim,debug_flip=debug_flip,$ + _EXTRA = EXTRA + + n_ant_pol=Max(antenna.n_pol) + nfreq_bin=Max(antenna.nfreq_bin) + IF N_Elements(file_path_fhd) EQ 0 THEN file_path_fhd='' + n_tile=obs.n_tile + beam_model_version=Max(antenna.model_version) + xvals_instrument=za_arr*Sin(az_arr*!DtoR) + yvals_instrument=za_arr*Cos(az_arr*!DtoR) + horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) + freq_center=antenna[0].freq ;all need to be identical, so just use the first + speed_light=299792458. ;speed of light, in meters/second + icomp=DComplex(0,1) + + ; Read the instrumental pol Jones matrix in az za + if (min((*obs.baseline_info).freq) GT 189.06250e6) OR (max((*obs.baseline_info).freq) LT 128.85550e6) then begin + message, 'No LOFAR beam model found for this frequency range' + endif else begin + if (min((*obs.baseline_info).freq) GE 159.59950e6) then begin + file_path_J_matrix=filepath('lofar_hamaker_jones_highband.h5',root=rootdir('FHD'),sub='instrument_config') + endif else begin + if (max((*obs.baseline_info).freq) LE 159.59950e6) then begin + file_path_J_matrix=filepath('lofar_hamaker_jones_lowband.h5',root=rootdir('FHD'),sub='instrument_config') + endif else begin + message, 'Frequency range spans both LOFAR high and low band. Please split the observation or develop the code.' + endelse + endelse + endelse + + file_id = H5F_OPEN(file_path_J_matrix) ;Get file ID + + ; Define the dataset names + dataset_names = ['gx', 'gy', 'Dx', 'Dy', 'gx_element', 'gy_element', 'Dx_element', 'Dy_element', $ + 'az', 'freqs', 'nside', 'radec_reso'] + ; also present in file but not required at this stage: 'za' (zenith angle), 'station' (station index) + + ; Read the LOFAR h5 file + FOR i = 0, N_ELEMENTS(dataset_names)-1 DO BEGIN + dataset_id = H5D_OPEN(file_id, dataset_names[i]) + data = H5D_READ(dataset_id) + H5D_CLOSE, dataset_id + + ; Assign to appropriately named variable + ; Reconstruct complex values + ; "_element" refers to a single element beam, used for phase flip information in D calculation + CASE dataset_names[i] OF + 'gx': gx = (complex(data.R, data.I)) + 'gy': gy = (complex(data.R, data.I)) + 'Dx': Dx = (complex(data.R, data.I)) + 'Dy': Dy = (complex(data.R, data.I)) + 'gx_element': gx_element = (complex(data.R, data.I)) + 'gy_element': gy_element = (complex(data.R, data.I)) + 'Dx_element': Dx_element = (complex(data.R, data.I)) + 'Dy_element': Dy_element = (complex(data.R, data.I)) + 'az': az = data + ;'za': za = data ;not used + 'freqs': freqs = data + 'nside': nside = data + 'radec_reso': radec_reso = data + ENDCASE + ENDFOR + + ; Get rotation angles of each LOFAR station compared to the first LOFAR station + rotation_filepath=filepath('lofar_hba_rotation_angles.txt',root=rootdir('FHD'),subdir='instrument_config') + textfast,rotation_angle,/read,file_path=rotation_filepath + + nside2 = ulong(nside)^2 + + ; Find pixels used at a single frequency + nan_inds = where(~finite(az),nan_count, complement=pix_use_input) + ; Remove NaN values to make interpolation easier across all freqs + nan_inds = where(~finite(gx),nan_count) + if nan_count GT 0 then begin + gx[nan_inds] = 0 + Dx[nan_inds] = 0 + Dy[nan_inds] = 0 + gy[nan_inds] = 0 + endif + + ; The inputs from the hd5 file were created with python, and thus have + ; row-col major instead of IDL's column-row major. + ; We need to transpose the dimensions of the image to match IDL's expectations. + n_freq_input = N_Elements(gx[0,*]) + gx = reform(transpose(reform(gx,nside,nside,n_freq_input),[1, 0, 2]),nside2,n_freq_input) + gy = reform(transpose(reform(gy,nside,nside,n_freq_input),[1, 0, 2]),nside2,n_freq_input) + Dx = reform(transpose(reform(Dx,nside,nside,n_freq_input),[1, 0, 2]),nside2,n_freq_input) + Dy = reform(transpose(reform(Dy,nside,nside,n_freq_input),[1, 0, 2]),nside2,n_freq_input) + ; Transpose image for element beams, and leave in 2D since no freq interp required + gx_element = transpose(reform(gx_element,nside,nside)) + gy_element = transpose(reform(gy_element,nside,nside)) + Dx_element = transpose(reform(Dx_element,nside,nside)) + Dy_element = transpose(reform(Dy_element,nside,nside)) + + ; Get the astrometry structure for the beam pointed at the phase center + hor2eq,obs.phasedec,obs.phasera,obs.jd0,beam_center_ra,beam_center_dec,ha_out,lat=obs.lat,lon=obs.lon,/precess,/nutate + projection_slant_orthographic,astr=beam_astr,degpix=radec_reso,obsra=beam_center_ra,obsdec=beam_center_dec,zenra=beam_center_ra,zendec=beam_center_dec,$ + dimension=nside,elements=nside,obsx=obsx,obsy=obsy,zenx=beam_zenx,zeny=beam_zeny,phasera=beam_center_ra,phasedec=beam_center_dec,$ + epoch=2000.,JDate=obs.JD0 + + ; Using this astrometry, find the x,y positions of the instrumental beam pixels + Hor2Eq,90.-za_arr[pix_use],az_arr[pix_use],obs.jd0,ra_use,dec_use,lat=obs.lat,lon=obs.lon,alt=obs.alt,precess=1,/nutate, refract=0 + ad2xy,ra_use,dec_use,beam_astr,x_interp,y_interp + + ; Find indices where a flip in phase occurs, given a 0.1 rad tolerance. + ; To define a positive delay from the dipole arm, these J matrix indices will be multiplied by -1. + delay_element = atan(imaginary(float(gx_element[pix_use_input])), real_part(float(gx_element[pix_use_input]))) + flip_inds_J00 = where(delay_element GT (!pi-0.1)) + delay_element = atan(imaginary(float(Dx_element[pix_use_input])), real_part(float(Dx_element[pix_use_input]))) + flip_inds_J01 = where(delay_element GT (!pi-0.1)) + delay_element = atan(imaginary(float(Dy_element[pix_use_input])), real_part(float(Dy_element[pix_use_input]))) + flip_inds_J10 = where(delay_element GT (!pi-0.1)) + delay_element = atan(imaginary(float(gy_element[pix_use_input])), real_part(float(gy_element[pix_use_input]))) + flip_inds_J11 = where(delay_element GT (!pi-0.1)) + + ; Calculate center index of input for rotation algorithm + center_ind = (ulong(nside)) / 2 + + ; initialize K matrix for instrumental pixels + K_proj_inst = PTRARR(2, 2) + + ; Loop over instrumental frequencies + FOR freq_i = 0, nfreq_bin-1 DO BEGIN + + J00 = complex(fltarr(N_elements(pix_use_input))) + J11 = complex(fltarr(N_elements(pix_use_input))) + J01 = complex(fltarr(N_elements(pix_use_input))) + J10 = complex(fltarr(N_elements(pix_use_input))) + + for pix_i=0L, N_elements(pix_use_input)-1 do begin + J00[pix_i] = interpol(gx[pix_use_input[pix_i],*],freqs, freq_center[freq_i],/quadratic) + J11[pix_i] = interpol(gy[pix_use_input[pix_i],*],freqs, freq_center[freq_i],/quadratic) + J01[pix_i] = interpol(Dx[pix_use_input[pix_i],*],freqs, freq_center[freq_i],/quadratic) + J10[pix_i] = interpol(Dy[pix_use_input[pix_i],*],freqs, freq_center[freq_i],/quadratic) + endfor + + ; Solve for the magnitude of F, the diagonal amplitude, hence only two elements defined + F_mag = [[sqrt(abs(J00)^2 + abs(J01)^2)], $ + [sqrt(abs(J10)^2 + abs(J11)^2)]] + + ; Get the combination of F_phase and K (pseudo, or old K), + ; which is easily solvable because the inverse of F is the reciprocal + pseudo_K = [[[J00 / F[*,0]], [J01 / F[*,0]]], $ + [[J10 / F[*,1]], [J11 / F[*,1]]]] + + J00[flip_inds_J00] *= -1 + J01[flip_inds_J01] *= -1 + J10[flip_inds_J10] *= -1 + J11[flip_inds_J11] *= -1 + + ; F_phase is the diagonal time delay matrix, the combined phase received by p and q components + F_phase = [[exp(icomp * atan(imaginary(J00+J01), real_part(J00+J01)))], $ + [exp(icomp * atan(imaginary(J10+J11), real_part(J10+J11)))]] + + ; K_proj is the polarisation-dependent (projection) response which should be tile-independent + K_proj = [[[ pseudo_K[*,0,0] / D[*, 0] ], [ pseudo_K[*,0,1] / D[*, 0] ]], $ + [[ pseudo_K[*,1,0] / D[*, 1] ], [ pseudo_K[*,1,1] / D[*, 1] ]]] + + ; Dummy matrix for intermediate steps + input_matrix = Dcomplexarr(nside,nside) + + ;Interpolate K_proj to the instrumental pixel grid + for i = 0,1 do begin + for j = 0,1 do begin + input_matrix[pix_use_input] = K_proj[*, i, j] + K_proj_inst[i, j] = Ptr_new(INTERPOLATE(input_matrix, x_interp, y_interp)) + endfor + endfor + + antenna[*].K_proj[*,*,freq_i] = K_proj_inst + + for tile_i=0, n_tile-1 do begin + ; For each tile, rotate the F_mag and F_phase matrices. + ; This is way around having a seperate beam model file for each LOFAR station. + ; Generalised input_matrix and output_matrix used to reduce RAM overhead. + + input_matrix[pix_use_input] = F_mag[*,0] + output_matrix = rot_cubic_complex(input_matrix, rotation_angle[tile_i], x_interp, y_interp, center_ind) + (antenna[tile_i].F_mag[0, freq_i]) = ptr_new(output_matrix) + + input_matrix[pix_use_input] = F_phase[*,0] + output_matrix = rot_cubic_complex(input_matrix, rotation_angle[tile_i], x_interp, y_interp, center_ind) + (antenna[tile_i].F_phase[0, freq_i]) = ptr_new(output_matrix) + + input_matrix[pix_use_input] = F_mag[*,1] + output_matrix = rot_cubic_complex(input_matrix, rotation_angle[tile_i], x_interp, y_interp, center_ind) + (antenna[tile_i].F_mag[1, freq_i]) = ptr_new(output_matrix) + + input_matrix[pix_use_input] = F_phase[*,1] + output_matrix = rot_cubic_complex(input_matrix, rotation_angle[tile_i], x_interp, y_interp, center_ind) + (antenna[tile_i].F_phase[1, freq_i]) = ptr_new(output_matrix) + + ; Update combined Jones matrix ( F D K ) to reflect rotation of delays and amplitude + (antenna[tile_i].jones[0, 0, freq_i]) = Ptr_new( (*antenna[tile_i].F_mag[0, freq_i]) * (*antenna[tile_i].F_phase[0, freq_i]) * (*antenna[tile_i].K_proj[0, 0, freq_i]) ) + (antenna[tile_i].jones[0, 1, freq_i]) = Ptr_new( (*antenna[tile_i].F_mag[0, freq_i]) * (*antenna[tile_i].F_phase[0, freq_i]) * (*antenna[tile_i].K_proj[0, 1, freq_i]) ) + (antenna[tile_i].jones[1, 0, freq_i]) = Ptr_new( (*antenna[tile_i].F_mag[1, freq_i]) * (*antenna[tile_i].F_phase[1, freq_i]) * (*antenna[tile_i].K_proj[1, 0, freq_i]) ) + (antenna[tile_i].jones[1, 1, freq_i]) = Ptr_new( (*antenna[tile_i].F_mag[1, freq_i]) * (*antenna[tile_i].F_phase[1, freq_i]) * (*antenna[tile_i].K_proj[1, 1, freq_i]) ) + + endfor + + endfor + + ; Group_id identifies which tiles have the same F D. + ; By definition, all lofar stations have a unique F D. + for tile_i=0, n_tile-1 do begin + antenna[tile_i].group_id = [tile_i, tile_i] + endfor + + RETURN,antenna +END diff --git a/instrument_config/lofar_beam_setup_init.pro b/instrument_config/lofar_beam_setup_init.pro new file mode 100644 index 00000000..90540318 --- /dev/null +++ b/instrument_config/lofar_beam_setup_init.pro @@ -0,0 +1,47 @@ +FUNCTION lofar_beam_setup_init, obs, antenna_str, antenna_size=antenna_size, dead_dipole_list=dead_dipole_list, $ + dipole_mutual_coupling_factor=dipole_mutual_coupling_factor, antenna_spacing=antenna_spacing, flag_dead_dipoles=flag_dead_dipoles + +; LOFAR HBA tile: 4x4 dipoles, ~1.25m spacing +n_ant_pol = antenna_str.n_pol +n_tiles = obs.n_tile +n_dipoles = 16 +nfreq_bin = antenna_str.nfreq_bin + +IF N_Elements(dipole_mutual_coupling_factor) EQ 0 THEN dipole_mutual_coupling_factor = 1 +IF N_Elements(antenna_size) EQ 0 THEN antenna_size = 5. ; meters (approximate tile size) +IF NOT Keyword_Set(antenna_spacing) THEN antenna_spacing = 1.25 ; meters (LOFAR HBA tile spacing) +antenna_height = 0.5 ; meters (typical HBA dipole height) +speed_light = 299792458. +base_delay_unit = 4.35E-10 + +; Dipole positions (centered) +xc_arr0 = Reform((meshgrid(4,4,1)) * antenna_spacing, 16) +xc_arr = xc_arr0 - Mean(xc_arr0) +yc_arr0 = Reform(Reverse(meshgrid(4,4,2),2) * antenna_spacing, 16) +yc_arr = yc_arr0 - Mean(yc_arr0) +zc_arr = Fltarr(16) + antenna_height + +antenna_coords = Ptrarr(3) +antenna_coords[0] = Ptr_new(xc_arr) +antenna_coords[1] = Ptr_new(yc_arr) +antenna_coords[2] = Ptr_new(zc_arr) + +; Read rotation angles from file +tile_names = (*obs.baseline_info).tile_names +file_path_rotations=filepath('lofar_hba_rotation_angles.txt',root=rootdir('FHD'),sub='instrument_config') +textfast,rotation_angles,/read,file_path=file_path_rotations + +; Set up antenna structure +antenna_str.n_ant_elements = n_dipoles +antenna_str.size_meters = antenna_size +antenna_str.coords = antenna_coords +antenna_str.height = antenna_height + +antenna = replicate(antenna_str, n_tiles) + +for tile_i=0, n_tiles-1 do begin + antenna[tile_i].rotation = rotation_angles[tile_i] +endfor + +RETURN, antenna +END \ No newline at end of file diff --git a/instrument_config/lofar_hamaker_jones_highband.h5 b/instrument_config/lofar_hamaker_jones_highband.h5 new file mode 100644 index 00000000..f4acd527 Binary files /dev/null and b/instrument_config/lofar_hamaker_jones_highband.h5 differ diff --git a/instrument_config/lofar_hamaker_jones_lowband.h5 b/instrument_config/lofar_hamaker_jones_lowband.h5 new file mode 100644 index 00000000..fd2d75ca Binary files /dev/null and b/instrument_config/lofar_hamaker_jones_lowband.h5 differ diff --git a/instrument_config/lofar_hba_rotation_angles.txt b/instrument_config/lofar_hba_rotation_angles.txt new file mode 100644 index 00000000..5079fce1 --- /dev/null +++ b/instrument_config/lofar_hba_rotation_angles.txt @@ -0,0 +1,70 @@ +0.0 +-45.35586983062188 +-175.04674252215793 +-176.35881474209867 +177.20405545573016 +171.61675151907065 +-159.97238614926368 +-169.0348033148621 +177.13759477388825 +170.95110979511819 +147.87117903780128 +160.79800623878972 +165.96375653207352 +152.0515655107047 +136.35319195377969 +123.69006752597979 +-168.02673220153764 +-157.7063708403059 +132.8622502744911 +124.91813353674416 +-134.61199839498644 +-137.18751181806363 +53.022403647327494 +59.93141717813756 +120.13135624164656 +121.63562804850005 +165.4177137465887 +167.9809975454555 +-142.66394875324377 +-148.5901626376359 +-126.6070748126075 +-122.66886053401196 +-94.34050330614359 +-85.306440092157 +144.24611274556327 +148.23147273706675 +107.75635034065408 +109.2434308549002 +104.09778563492826 +101.82530644714639 +7.594643368591445 +-5.9245819819824925 +-49.436474485657364 +-49.739482012062844 +-96.95295746817392 +-97.22179966818898 +-166.02459081766824 +-167.0121326225314 +7.1614975956203795 +48.47735901884017 +72.3193540785628 +94.1896369496067 +-3.2046792040757697 +-30.722554531175557 +74.20116402361204 +28.71465916844453 +16.013605993880642 +7.337315257448099 +-73.77094609300647 +-70.74692280329062 +-45.746449271186194 +-68.6849556861943 +-133.42664260156545 +-153.7544359047596 +-93.43209043516859 +-151.4448395046497 +175.8965969359885 +-165.95256406786658 +151.18915808995314 +-72.07717844679397