From d068afeebbf8accca1f7d8081feeee98b71f373c Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Thu, 23 Mar 2023 16:44:04 -0700 Subject: [PATCH 01/16] non whitespace changes --- .../beam_modeling/fhd_struct_init_antenna.pro | 60 +++------ fhd_core/beam_modeling/import_az_el_beam.pro | 123 ++++++++++++++++++ .../beam_modeling/pyuvdata_beam_import.pro | 35 +++++ .../polarization/fhd_struct_init_jones.pro | 4 +- instrument_config/general_beam_setup_gain.pro | 40 ++++++ instrument_config/general_beam_setup_init.pro | 44 +++++++ 6 files changed, 265 insertions(+), 41 deletions(-) create mode 100644 fhd_core/beam_modeling/import_az_el_beam.pro create mode 100644 fhd_core/beam_modeling/pyuvdata_beam_import.pro create mode 100644 instrument_config/general_beam_setup_gain.pro create mode 100644 instrument_config/general_beam_setup_init.pro diff --git a/fhd_core/beam_modeling/fhd_struct_init_antenna.pro b/fhd_core/beam_modeling/fhd_struct_init_antenna.pro index a2f0b3d5a..e9f242c19 100644 --- a/fhd_core/beam_modeling/fhd_struct_init_antenna.pro +++ b/fhd_core/beam_modeling/fhd_struct_init_antenna.pro @@ -1,11 +1,13 @@ FUNCTION fhd_struct_init_antenna,obs,beam_model_version=beam_model_version,$ psf_resolution=psf_resolution,psf_intermediate_res=psf_intermediate_res,$ psf_image_resolution=psf_image_resolution,timing=timing,$ - psf_dim=psf_dim,psf_max_dim=psf_max_dim,beam_offset_time=beam_offset_time,debug_dim=debug_dim,$ + psf_dim=psf_dim,psf_max_dim=psf_max_dim,beam_offset_time=beam_offset_time,$ inst_tile_ptr=inst_tile_ptr,ra_arr=ra_arr,dec_arr=dec_arr,fractional_size=fractional_size,$ kernel_window=kernel_window,beam_per_baseline=beam_per_baseline,$ - beam_gaussian_decomp=beam_gaussian_decomp,beam_gauss_param_transfer=beam_gauss_param_transfer,$ - conserve_memory=conserve_memory,_Extra=extra + beam_gaussian_decomp=beam_gaussian_decomp,conserve_memory=conserve_memory,$ + import_pyuvdata_beam_filepath=import_pyuvdata_beam_filepath,$ + use_psf_resolution=use_psf_resolution,_Extra=extra + t0=Systime(1) IF keyword_set(conserve_memory) then begin @@ -13,43 +15,14 @@ IF keyword_set(conserve_memory) then begin ENDIF IF N_Elements(beam_model_version) EQ 0 THEN beam_model_version=1 +IF N_Elements(use_psf_resolution) GT 0 THEN psf_resolution=use_psf_resolution instrument=obs.instrument -tile_gain_fn=instrument+'_beam_setup_gain' ;mwa_beam_setup_gain -tile_init_fn=instrument+'_beam_setup_init' ;mwa_beam_setup_init - -;If other phases of the mwa are used, use the proper gain and init functions -IF STRMID(instrument,0,3) EQ 'mwa' THEN BEGIN - tile_gain_fn='mwa_beam_setup_gain' - tile_init_fn='mwa_beam_setup_init' -ENDIF if keyword_set(beam_gaussian_decomp) and keyword_set(kernel_window) then begin print, 'Gaussian decomposition cannot be used with modified kernel windows. Window not applied.' kernel_window=0 endif -;Default the parameter transfer if not set -if keyword_set(beam_gauss_param_transfer) then begin - ;Default to instrumental beam if not set - if (beam_gauss_param_transfer EQ 1) then beam_gauss_param_transfer = 'decomp' - - ;Set transfer to the instrumental beam or gaussian beam. Currently only available for the MWA - if (beam_gauss_param_transfer EQ 'decomp') or (beam_gauss_param_transfer EQ 'gauss') then begin - if instrument EQ 'mwa' then begin - pointing_num = mwa_get_pointing_number(obs,string=1) - beam_gauss_param_transfer = filepath(instrument + '_decomp_params_pointing' + pointing_num + '.sav',$ - root=rootdir('FHD'),subdir='instrument_config') - endif else begin - message, 'Gaussian decomposition parameter defaults not currently set for non-MWA instruments' - endelse - endif - - ;Match the psf size to the parameters from the file - psf_transfer = getvar_savefile(beam_gauss_param_transfer,'psf') - psf_resolution = psf_transfer.resolution - psf_dim = psf_transfer.dim -endif - n_tiles=obs.n_tile n_freq=obs.n_freq n_pol=obs.n_pol @@ -61,8 +34,6 @@ zenra=obs.zenra zendec=obs.zendec obsx=obs.obsx obsy=obs.obsy -;phasera=obs.phasera -;phasedec=obs.phasedec Jdate=obs.Jd0 IF Keyword_Set(beam_offset_time) THEN Jdate_use=Jdate+beam_offset_time/24./3600. ELSE Jdate_use=Jdate frequency_array=(*obs.baseline_info).freq @@ -103,6 +74,14 @@ antenna_str={n_pol:n_ant_pol,antenna_type:instrument,names:ant_names,model_versi 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 keyword_set(import_pyuvdata_beam_filepath) then begin + tile_gain_fn='general_beam_setup_gain' + tile_init_fn='general_beam_setup_init' +endif else begin + tile_gain_fn=instrument+'_beam_setup_gain' ;mwa_beam_setup_gain + tile_init_fn=instrument+'_beam_setup_init' ;mwa_beam_setup_init +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 antenna=Call_function(tile_init_fn[0],obs,antenna_str,_Extra=extra) ;mwa_beam_setup_init @@ -118,7 +97,6 @@ IF ~Keyword_Set(psf_dim) THEN $ psf_dim=Ceil((Max(antenna.size_meters)*2.*Max(frequency_array)/speed_light)/kbinsize) psf_dim=Ceil(psf_dim/2.)*2. ;dimension MUST be even ;reset psf_dim if cetain conditions met. -if keyword_set(debug_dim) then psf_dim=28. if keyword_set(kernel_window) then psf_dim=18. IF Keyword_Set(psf_max_dim) THEN BEGIN @@ -142,9 +120,10 @@ undefine, xvals_celestial, yvals_celestial valid_i=where(Finite(ra_arr),n_valid) ra_use=ra_arr[valid_i] dec_use=dec_arr[valid_i] + if ~keyword_set(beam_per_baseline) then undefine, dec_arr, ra_arr -; Split the apply astrometry step into mutiple loops if the image dimension is large to reduce +; Split the apply astrometry step into mutiple loops if the image dimension is large to reduce ; memory footprint if keyword_set(conserve_memory) then begin ; calculate bytes required @@ -189,6 +168,9 @@ Eq2Hor,ra_use,dec_use,Jdate_use,alt_arr1,az_arr1,lat=obs.lat,lon=obs.lon,alt=obs za_arr=fltarr(psf_image_dim,psf_image_dim)+90. & za_arr[valid_i]=90.-alt_arr1 az_arr=fltarr(psf_image_dim,psf_image_dim) & az_arr[valid_i]=az_arr1 undefine, ra_use, dec_use, alt_arr1, az_arr1 +;Get pixels which fall within the horizon +horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) +antenna.pix_use=ptr_new(pix_use) if keyword_set(kernel_window) then begin print, 'Applying a modified gridding kernel. Beam is no longer instrumental. Do not use for calibration.' @@ -220,11 +202,11 @@ antenna.pix_use=ptr_new(pix_use) if N_elements(instrument) GT 1 then begin for inst_i=0, N_elements(instrument)-1 do begin antenna_temp = pointer_copy(antenna) - antenna_temp=Call_function(tile_gain_fn[inst_i],obs,antenna_temp,za_arr=za_arr,az_arr=az_arr,psf_image_dim=psf_image_dim,Jdate_use=Jdate_use,_Extra=extra) ;mwa_beam_setup_gain + antenna_temp=Call_function(tile_gain_fn[inst_i],obs,antenna_temp,za_arr=za_arr,az_arr=az_arr,psf_image_dim=psf_image_dim,Jdate_use=Jdate_use,import_pyuvdata_beam_filepath=import_pyuvdata_beam_filepath,_Extra=extra) ;mwa_beam_setup_gain antenna[*inst_tile_ptr[inst_i]] = pointer_copy(antenna_temp[*inst_tile_ptr[inst_i]]) antenna[*inst_tile_ptr[inst_i]].antenna_type = instrument[inst_i] ;if more than one instrument, assign the correct antenna type for each subset for metadata purposes 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 +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,import_pyuvdata_beam_filepath=import_pyuvdata_beam_filepath,_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) diff --git a/fhd_core/beam_modeling/import_az_el_beam.pro b/fhd_core/beam_modeling/import_az_el_beam.pro new file mode 100644 index 000000000..4826d0655 --- /dev/null +++ b/fhd_core/beam_modeling/import_az_el_beam.pro @@ -0,0 +1,123 @@ +FUNCTION import_az_el_beam, obs, antenna, file_path_J_matrix,$ + za_arr=za_arr,az_arr=az_arr,psf_image_dim=psf_image_dim +; Read in an e-field antenna beam stored in evenly-spaced azimuth-elevation (and frequency) + +n_ant_pol=Max(antenna.n_pol) +nfreq_bin=Max(antenna.nfreq_bin) +n_tile=obs.n_tile +xvals_instrument=za_arr*Sin(az_arr*!DtoR) +yvals_instrument=za_arr*Cos(az_arr*!DtoR) +freq_center=antenna[0].freq ;all need to be identical, so just use the first +icomp=Complex(0,1) + +fits_info,file_path_J_matrix,/silent,n_ext=n_ext +;Cols of J matrix: theta phi real(Jxt(t,p)) imag(Jxt(t,p)) real(Jxp(t,p)) imag(Jxp(t,p)) real(Jyt(t,p)) imag(Jyt(t,p)) real(Jyp(t,p)) imag(Jyp(t,p))) +;Where theta is the zenith angle, phi is angle measured clockwise from +east direction looking down +; Jxt is the Jones mapping unit vec in theta (t) direction to the x (east-west) dipole etc +n_ext+=1 ;n_ext starts counting AFTER the 0th extension, which it considers to be the main data unit, but we use that one too +freq_arr_Jmat=Fltarr(n_ext) + +;FOR ext_i=0,n_ext-1 DO BEGIN +; Jmat1=mrdfits(file_path_J_matrix,ext_i,header,status=status,/silent) +; IF ext_i EQ 0 THEN BEGIN +; n_ang=Float(sxpar(header,'NAXIS2')) +; Jmat_arr=Dcomplexarr(n_ext,n_ant_pol,n_ant_pol,n_ang) +; theta_arr=Dblarr(n_ext,n_ang) +; phi_arr=Dblarr(n_ext,n_ang) +; ENDIF +; theta_arr[ext_i,*]=Jmat1[0,*] ;zenith angle in degrees +; phi_arr[ext_i,*]=Jmat1[1,*] ;azimuth angle in degrees, clockwise from East +; FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO BEGIN +; Jmat_arr[ext_i,p_i,p_j,*]=Jmat1[2+p_i*2+p_j*4,*]+icomp*Jmat1[2+p_i*2+p_j*4+1,*] +; if keyword_set(debug_flip) then Jmat_arr[ext_i,p_j,p_i,*]=Jmat1[2+p_i*2+p_j*4,*]+icomp*Jmat1[2+p_i*2+p_j*4+1,*] +; ENDFOR +; freq_arr_Jmat[ext_i]=Float(sxpar(header,'FREQ')) ;in Hz +;ENDFOR +;theta_arr=median(theta_arr,dimension=1) ; all actually the same across freq, so reduce dimension +;phi_arr=median(phi_arr,dimension=1) ; all actually the same across freq, so reduce dimension + +;Begin hacks for pyuvdata beam +; CTYPE1: Azimuth in degrees, from 0 and increment of 1. Don't know origin: guessing clockwise from East +; CTYPE2: Zenith angle in degrees, from 0 and increment of 1 +; CTYPE3: Frequency in Hz, from 50,000,000 and increment of 1,000,000 +; CTYPE4: "Feed index", X then Y +; CTYPE5: Spectral window number. There is only one. +; CTYPE6: Basis vector index. guessing these are theta then phi +; CTYPE7: Complex, real then imaginary component. +Jmat0=mrdfits(file_path_J_matrix,0,header,status=status,/silent) +n_za_ang = sxpar(header,'naxis2') +n_az_ang = sxpar(header,'naxis1') +n_freq = sxpar(header, 'naxis3') +n_feed = sxpar(header, 'naxis4') +n_basis = sxpar(header, 'naxis6') +n_ang = n_za_ang*n_az_ang + +theta_arr0 = findgen(n_za_ang)*sxpar(header,'cdelt2') + sxpar(header, 'crval2') ;in degrees +phi_arr0 = findgen(n_az_ang)*sxpar(header,'cdelt1') + sxpar(header, 'crval1') ;in degrees +freq_arr_Jmat = findgen(sxpar(header,'naxis3'))*sxpar(header,'cdelt3') + sxpar(header, 'crval3') ; in Hz + +az_ind_arr = Reform(Fix(meshgrid(n_za_ang, n_az_ang, 2)), n_ang) +za_ind_arr = Reform(Fix(meshgrid(n_za_ang, n_az_ang, 1)), n_ang) +theta_arr = theta_arr0[za_ind_arr] +phi_arr = phi_arr0[az_ind_arr] + +Jmat_arr=Dcomplexarr(n_freq,n_basis,n_feed,n_ang) +FOR freq_i=0,n_freq-1 DO BEGIN + FOR f_i=0,n_feed-1 DO BEGIN + FOR b_i=0,n_basis-1 DO BEGIN + Jmat1 = Reform(Jmat0[*, *, freq_i, f_i, 0, b_i, 0] + icomp*Jmat0[*, *, freq_i, f_i, 0, b_i, 1]) + Jmat1 = Transpose(Jmat1) + Jmat_arr[freq_i,f_i,b_i,*] = Reform(Jmat1, n_ang) + ENDFOR + ENDFOR +ENDFOR + +phi_arr=270.-phi_arr ;change azimuth convention + + +Jmat_interp=Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin) +FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO Jmat_interp[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(n_ang)) +FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR a_i=0L,n_ang-1 DO BEGIN + Jmat_single_ang=Interpol(Jmat_arr[*,p_i,p_j,a_i],freq_arr_Jmat,freq_center);*norm_factor + FOR freq_i=0L,nfreq_bin-1 DO (*Jmat_interp[p_i,p_j,freq_i])[a_i]=Jmat_single_ang[freq_i] +ENDFOR + +xv_model=theta_arr*Sin(phi_arr*!DtoR) +yv_model=theta_arr*Cos(phi_arr*!DtoR) + +zenith_i=where(theta_arr EQ 0,n_zenith) + +horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) +horizon_mask=fltarr(psf_image_dim,psf_image_dim)+1 +IF n_horizon_test GT 0 THEN horizon_mask[horizon_test]=0 +Jones_matrix=Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin) +FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO $ + Jones_matrix[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(psf_image_dim,psf_image_dim)) + +interp_res=obs.degpix +angle_slice_i0=Uniq(phi_arr) +n_ang_slice=N_Elements(angle_slice_i0) +n_zen_slice=angle_slice_i0[0]+1 +az_ang_in=phi_arr[angle_slice_i0] +zen_ang_in=theta_arr[0:angle_slice_i0[0]] +zen_ang_out=Findgen(Ceil(90./interp_res)+1)*interp_res +az_ang_out=Findgen(Ceil(360./interp_res)+1)*interp_res+Round(Min(az_ang_in)) +n_zen_ang=N_Elements(zen_ang_out) +n_az_ang=N_Elements(az_ang_out) + +xv_model=zen_ang_out*Sin(az_ang_out*!DtoR) +yv_model=zen_ang_out*Cos(az_ang_out*!DtoR) + +zen_ang_inst=Sqrt(xvals_instrument[pix_use]^2+yvals_instrument[pix_use]^2) +az_ang_inst=Atan(yvals_instrument[pix_use],xvals_instrument[pix_use])*!Radeg+180. + +FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO BEGIN + Jmat_use=Reform(*Jmat_interp[p_i,p_j,freq_i],n_zen_slice,n_ang_slice) + Expand,Jmat_use,n_zen_ang,n_az_ang,Jmat_single + (*Jones_matrix[p_i,p_j,freq_i])[pix_use]=Interpolate(Jmat_single,zen_ang_inst/interp_res,az_ang_inst/interp_res,cubic=-0.5) + debug_point=1 +ENDFOR +debug_point=1 + +RETURN, Jones_matrix +END diff --git a/fhd_core/beam_modeling/pyuvdata_beam_import.pro b/fhd_core/beam_modeling/pyuvdata_beam_import.pro new file mode 100644 index 000000000..30dd7e1dd --- /dev/null +++ b/fhd_core/beam_modeling/pyuvdata_beam_import.pro @@ -0,0 +1,35 @@ +FUNCTION pyuvdata_beam_import, obs, antenna_str, pyuvdata_filepath,$ + za_arr=za_arr, az_arr=az_arr, psf_image_dim=psf_image_dim, pix_use=pix_use + +lun = fxposit(pyuvdata_filepath, 0,/readonly) +MRD_HREAD, lun, primary_header, /silent + +coordsys = strtrim(sxpar(primary_header,'coordsys'),2) +CASE coordsys OF + "az_za": BEGIN + Jones_matrix = import_az_el_beam(obs, antenna_str, pyuvdata_filepath,$ + za_arr=za_arr, az_arr=az_arr, psf_image_dim=psf_image_dim) + END + "orthoslant_zenith": BEGIN + message,"COORDSYS type 'orthoslant_zenith' not yet supported!"" + END + "healpix": BEGIN + message,"COORDSYS type 'healpix' not yet supported!" + END + ELSE: message,"Unsupported COORDSYS found in pyuvdata file!" +ENDCASE + +for pol1=0,(size(Jones_matrix))[1]-1 do begin + for pol2=0,(size(Jones_matrix))[2]-1 do begin + if (size(Jones_matrix))[0] gt 2 then begin + for ind=0,(size(Jones_matrix))[3]-1 do begin + *Jones_matrix[pol1, pol2, ind] = (*Jones_matrix[pol1, pol2, ind])[pix_use] + endfor + endif else begin + *Jones_matrix[pol1, pol2] = (*Jones_matrix[pol1, pol2])[pix_use] + endelse + endfor +endfor + +RETURN, Jones_matrix +END diff --git a/fhd_core/polarization/fhd_struct_init_jones.pro b/fhd_core/polarization/fhd_struct_init_jones.pro index 30f0bac08..2f9121b55 100644 --- a/fhd_core/polarization/fhd_struct_init_jones.pro +++ b/fhd_core/polarization/fhd_struct_init_jones.pro @@ -51,8 +51,8 @@ apply_astrometry, obs, ra_arr=obs.obsra, dec_arr=obs.lat, x_arr=x_t, y_arr=y_t, apply_astrometry, obs, dec_arr=lat, x_arr=x_t, y_arr=y_t, /xy2ad, /ignore_refraction obs_temp = fhd_struct_update_obs(obs, beam_nfreq_avg=obs.n_freq) ; Use only one average Jones matrix, not one per frequency -antenna=fhd_struct_init_antenna(obs_temp,beam_model_version=beam_model_version,psf_resolution=1.,$ - psf_dim=obs.dimension,psf_intermediate_res=1.,psf_image_resolution=1.,timing=t_ant) +antenna=fhd_struct_init_antenna(obs_temp,beam_model_version=beam_model_version,use_psf_resolution=1.,$ + psf_dim=obs.dimension,psf_intermediate_res=1.,psf_image_resolution=1.,timing=t_ant,_Extra=extra) Jones=rotate_jones_matrix(obs, antenna[0]) ; Calculate the normalization diff --git a/instrument_config/general_beam_setup_gain.pro b/instrument_config/general_beam_setup_gain.pro new file mode 100644 index 000000000..d0464ee31 --- /dev/null +++ b/instrument_config/general_beam_setup_gain.pro @@ -0,0 +1,40 @@ +FUNCTION general_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,Jdate_use=Jdate_use,$ + import_pyuvdata_beam_filepath=import_pyuvdata_beam_filepath,_Extra=extra + + if ~keyword_set(import_pyuvdata_beam_filepath) then import_pyuvdata_beam_filepath=extra.import_pyuvdata_beam_filepath + + n_ant_pol=Max(antenna.n_pol) + nfreq_bin=Max(antenna.nfreq_bin) + pix_use = *antenna[0].pix_use + 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_interp=za_arr*Sin(az_arr*!DtoR)/obs.degpix+obs.dimension/2. + yvals_interp=za_arr*Cos(az_arr*!DtoR)/obs.degpix+obs.elements/2. + freq_center=antenna[0].freq ;all need to be identical, so just use the first + pyuvdata_filepath=import_pyuvdata_beam_filepath + + ;calculate group identifications (used to set pointers to identical models) + FOR pol_i=0,n_ant_pol-1 DO BEGIN + gi=0 + n_ungrouped=n_tile + ungrouped_i=where(antenna.group_id[pol_i] EQ -1,n_ungrouped) + WHILE n_ungrouped GT 0 DO BEGIN + ref_i=ungrouped_i[0] + antenna[ref_i].group_id[pol_i]=gi + FOR ug_i=1L,n_ungrouped-1 DO IF Total(*antenna[ungrouped_i[ug_i]].gain[pol_i] - *antenna[ref_i].gain[pol_i]) EQ 0 THEN antenna[ungrouped_i[ug_i]].group_id[pol_i]=gi + ungrouped_i=where(antenna.group_id[pol_i] EQ -1,n_ungrouped) + gi+=1 + ENDWHILE + ENDFOR + + ;get the instrumental pol Jones matrix + print,"Reading in: " + pyuvdata_filepath + Jones_matrix = pyuvdata_beam_import(obs, antenna, pyuvdata_filepath,$ + za_arr=za_arr, az_arr=az_arr, psf_image_dim=psf_image_dim, pix_use=pix_use) + + antenna.jones=Jones_matrix + + RETURN,antenna +END diff --git a/instrument_config/general_beam_setup_init.pro b/instrument_config/general_beam_setup_init.pro new file mode 100644 index 000000000..d740c445d --- /dev/null +++ b/instrument_config/general_beam_setup_init.pro @@ -0,0 +1,44 @@ +FUNCTION general_beam_setup_init,obs,antenna_str,_Extra=extra + + ;copied from hera_beam_setup_init + ; + ; + ;polarization 0: x, 1: y + + n_ant_pol=antenna_str.n_pol + n_tiles=obs.n_tile + n_dipoles=1 + nfreq_bin=antenna_str.nfreq_bin + + IF N_Elements(antenna_size) EQ 0 THEN antenna_size=10 ;meters A GUESS + antenna_height=1.5 ;meters A GUESS + velocity_factor=0.673 ;use MWA number (ignored anyway) + + xc_arr=Fltarr(n_dipoles) + yc_arr=Fltarr(n_dipoles) + zc_arr=Fltarr(n_dipoles) + + 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) + + delay_settings=Ptr_new(Fltarr(n_dipoles)) + + freq_center=antenna_str.freq + FOR i=0L,N_Elements(antenna_str.coupling)-1 DO antenna_str.coupling[i]=Ptr_new(Complex(Identity(n_dipoles))) + + base_gain=fltarr(n_dipoles)+1. + gain_arr=Ptrarr(n_ant_pol) + FOR pol_i=0,n_ant_pol-1 DO gain_arr[pol_i]=Ptr_new(Rebin(reform(base_gain,1,n_dipoles),nfreq_bin,n_dipoles,/sample)) + + antenna_str.n_ant_elements=n_dipoles + antenna_str.size_meters=antenna_size + antenna_str.coords=antenna_coords + antenna_str.height=antenna_height + antenna_str.delays=delay_settings + antenna=replicate(antenna_str,n_tiles) + FOR t_i=0L,n_tiles-1 DO antenna[t_i].gain=Pointer_copy(gain_arr) + + RETURN, antenna +END \ No newline at end of file From ce58b771b0d7c50a716ac55b5f5b5fcd992869d8 Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Thu, 23 Mar 2023 17:32:18 -0700 Subject: [PATCH 02/16] update dictionary --- dictionary.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/dictionary.md b/dictionary.md index 188733771..004a95261 100644 --- a/dictionary.md +++ b/dictionary.md @@ -22,7 +22,7 @@ This is a work in progress; please add keywords as you find them in alphabetical -*Dependency*: `model_visibilities` must be set
-*Default*: 0.05, or 0.01 if `allow_sidelobe_sources` set
-**beam_model_version**: a number that indicates the tile beam model calculation. This is dependent on the instrument, and specific calculations are carried out in `_beam_setup_gain.pro`. For the MWA, there are currently three options: 0) !Q, 1) a Hertzian dipole as prescribed by Cheng 1992 and Balanis 1989, 2) the average embedded element model from Sutinjo 2015. For PAPER, there are currently two options: 1) !Q, 2) !Q. For HERA, there are currently two options: 2) 2016 version by Dave Deboer saved as cross and co-pol, 0) an earlier version !Q saved as X and Y pols.
+**beam_model_version**: a number that indicates the tile beam model calculation. This is dependent on the instrument, and specific calculations are carried out in `_beam_setup_gain.pro`. For the MWA, there are currently three options: 0) !Q, 1) a Hertzian dipole as prescribed by Cheng 1992 and Balanis 1989, 2) the average embedded element model from Sutinjo 2015. For PAPER, there are currently two options: 1) !Q, 2) !Q. For HERA, there are currently two options: 2) 2016 version by Dave Deboer saved as cross and co-pol, 0) an earlier version !Q saved as X and Y pols. This keyword is ignored if import_pyuvdata_beam_filepath is set.
-*MWA range*: 0, 1 (or anything else captured in the `else` statement), 2
-*PAPER range*: 1 (or anything else captured in the `else` statement), 2
-*HERA range*: 2 (or anything else captured in the `else` statement)
@@ -49,6 +49,9 @@ This is a work in progress; please add keywords as you find them in alphabetical -*Default*: 1
-*eor_wrapper_defaults*: 1
+**import_pyuvdata_beam_filepath**: path to a .beamfits file that contains a beam model formatted according to the standards of [pyuvdata](https://github.com/RadioAstronomySoftwareGroup/pyuvdata). If set, this beam model will be used instead of an internal beam model specified with beam_model_version.
+ -*Default*: not set
+ **inst_tile_ptr**: a pointer array to designate which tile indices belong to which instrument. The order of the pointer array is assumed to match the order of instruments specified in the keyword string array `instrument`. Only used if there is more than one instrument supplied. Tiles are numbered from 0.
-*Example*: `inst_tile_ptr = PTRARR(2,/allocate)`
`*inst_tile_ptr[0] = [0,1,2,3]`
From f460406a7b64c077d634ef61dde9496fb7921303 Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Mon, 3 Apr 2023 13:56:19 -0700 Subject: [PATCH 03/16] remove obsolete code --- fhd_core/beam_modeling/import_az_el_beam.pro | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/fhd_core/beam_modeling/import_az_el_beam.pro b/fhd_core/beam_modeling/import_az_el_beam.pro index 4826d0655..63d7de6bf 100644 --- a/fhd_core/beam_modeling/import_az_el_beam.pro +++ b/fhd_core/beam_modeling/import_az_el_beam.pro @@ -17,25 +17,6 @@ fits_info,file_path_J_matrix,/silent,n_ext=n_ext n_ext+=1 ;n_ext starts counting AFTER the 0th extension, which it considers to be the main data unit, but we use that one too freq_arr_Jmat=Fltarr(n_ext) -;FOR ext_i=0,n_ext-1 DO BEGIN -; Jmat1=mrdfits(file_path_J_matrix,ext_i,header,status=status,/silent) -; IF ext_i EQ 0 THEN BEGIN -; n_ang=Float(sxpar(header,'NAXIS2')) -; Jmat_arr=Dcomplexarr(n_ext,n_ant_pol,n_ant_pol,n_ang) -; theta_arr=Dblarr(n_ext,n_ang) -; phi_arr=Dblarr(n_ext,n_ang) -; ENDIF -; theta_arr[ext_i,*]=Jmat1[0,*] ;zenith angle in degrees -; phi_arr[ext_i,*]=Jmat1[1,*] ;azimuth angle in degrees, clockwise from East -; FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO BEGIN -; Jmat_arr[ext_i,p_i,p_j,*]=Jmat1[2+p_i*2+p_j*4,*]+icomp*Jmat1[2+p_i*2+p_j*4+1,*] -; if keyword_set(debug_flip) then Jmat_arr[ext_i,p_j,p_i,*]=Jmat1[2+p_i*2+p_j*4,*]+icomp*Jmat1[2+p_i*2+p_j*4+1,*] -; ENDFOR -; freq_arr_Jmat[ext_i]=Float(sxpar(header,'FREQ')) ;in Hz -;ENDFOR -;theta_arr=median(theta_arr,dimension=1) ; all actually the same across freq, so reduce dimension -;phi_arr=median(phi_arr,dimension=1) ; all actually the same across freq, so reduce dimension - ;Begin hacks for pyuvdata beam ; CTYPE1: Azimuth in degrees, from 0 and increment of 1. Don't know origin: guessing clockwise from East ; CTYPE2: Zenith angle in degrees, from 0 and increment of 1 From d9daba789a8e1b156865d114cb5b72cb7599de77 Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Mon, 3 Apr 2023 14:13:01 -0700 Subject: [PATCH 04/16] add comments and remove unnecessary code --- fhd_core/beam_modeling/import_az_el_beam.pro | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/fhd_core/beam_modeling/import_az_el_beam.pro b/fhd_core/beam_modeling/import_az_el_beam.pro index 63d7de6bf..10187b7f3 100644 --- a/fhd_core/beam_modeling/import_az_el_beam.pro +++ b/fhd_core/beam_modeling/import_az_el_beam.pro @@ -55,7 +55,7 @@ ENDFOR phi_arr=270.-phi_arr ;change azimuth convention - +; Interpolate in frequency Jmat_interp=Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin) FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO Jmat_interp[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(n_ang)) FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR a_i=0L,n_ang-1 DO BEGIN @@ -63,14 +63,11 @@ FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR a_i=0L,n_ang-1 DO BEGIN FOR freq_i=0L,nfreq_bin-1 DO (*Jmat_interp[p_i,p_j,freq_i])[a_i]=Jmat_single_ang[freq_i] ENDFOR -xv_model=theta_arr*Sin(phi_arr*!DtoR) -yv_model=theta_arr*Cos(phi_arr*!DtoR) - zenith_i=where(theta_arr EQ 0,n_zenith) -horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) -horizon_mask=fltarr(psf_image_dim,psf_image_dim)+1 -IF n_horizon_test GT 0 THEN horizon_mask[horizon_test]=0 +;horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) +;horizon_mask=fltarr(psf_image_dim,psf_image_dim)+1 +;IF n_horizon_test GT 0 THEN horizon_mask[horizon_test]=0 Jones_matrix=Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin) FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO $ Jones_matrix[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(psf_image_dim,psf_image_dim)) @@ -80,18 +77,17 @@ angle_slice_i0=Uniq(phi_arr) n_ang_slice=N_Elements(angle_slice_i0) n_zen_slice=angle_slice_i0[0]+1 az_ang_in=phi_arr[angle_slice_i0] -zen_ang_in=theta_arr[0:angle_slice_i0[0]] +;zen_ang_in=theta_arr[0:angle_slice_i0[0]] zen_ang_out=Findgen(Ceil(90./interp_res)+1)*interp_res az_ang_out=Findgen(Ceil(360./interp_res)+1)*interp_res+Round(Min(az_ang_in)) n_zen_ang=N_Elements(zen_ang_out) n_az_ang=N_Elements(az_ang_out) -xv_model=zen_ang_out*Sin(az_ang_out*!DtoR) -yv_model=zen_ang_out*Cos(az_ang_out*!DtoR) - +; Convert from x, y to azimuth, zenith angle coordinates zen_ang_inst=Sqrt(xvals_instrument[pix_use]^2+yvals_instrument[pix_use]^2) az_ang_inst=Atan(yvals_instrument[pix_use],xvals_instrument[pix_use])*!Radeg+180. +; Interpolate in azimulth and zenith angle FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO BEGIN Jmat_use=Reform(*Jmat_interp[p_i,p_j,freq_i],n_zen_slice,n_ang_slice) Expand,Jmat_use,n_zen_ang,n_az_ang,Jmat_single From d7b575aceb99a22037a3598d930a7bf72214bedb Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Mon, 3 Apr 2023 14:20:58 -0700 Subject: [PATCH 05/16] modify comments --- fhd_core/beam_modeling/import_az_el_beam.pro | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/fhd_core/beam_modeling/import_az_el_beam.pro b/fhd_core/beam_modeling/import_az_el_beam.pro index 10187b7f3..11c9a934f 100644 --- a/fhd_core/beam_modeling/import_az_el_beam.pro +++ b/fhd_core/beam_modeling/import_az_el_beam.pro @@ -17,14 +17,8 @@ fits_info,file_path_J_matrix,/silent,n_ext=n_ext n_ext+=1 ;n_ext starts counting AFTER the 0th extension, which it considers to be the main data unit, but we use that one too freq_arr_Jmat=Fltarr(n_ext) -;Begin hacks for pyuvdata beam -; CTYPE1: Azimuth in degrees, from 0 and increment of 1. Don't know origin: guessing clockwise from East -; CTYPE2: Zenith angle in degrees, from 0 and increment of 1 -; CTYPE3: Frequency in Hz, from 50,000,000 and increment of 1,000,000 -; CTYPE4: "Feed index", X then Y -; CTYPE5: Spectral window number. There is only one. -; CTYPE6: Basis vector index. guessing these are theta then phi -; CTYPE7: Complex, real then imaginary component. +; Read beam data from fits file +; Format is set by pyuvdata conventions, see https://github.com/RadioAstronomySoftwareGroup/pyuvdata for details Jmat0=mrdfits(file_path_J_matrix,0,header,status=status,/silent) n_za_ang = sxpar(header,'naxis2') n_az_ang = sxpar(header,'naxis1') @@ -65,7 +59,7 @@ ENDFOR zenith_i=where(theta_arr EQ 0,n_zenith) -;horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) +horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) ;horizon_mask=fltarr(psf_image_dim,psf_image_dim)+1 ;IF n_horizon_test GT 0 THEN horizon_mask[horizon_test]=0 Jones_matrix=Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin) From 993eb460d2c83cf7dbfbd74bf6f0c76603cae4c5 Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Mon, 3 Apr 2023 14:25:47 -0700 Subject: [PATCH 06/16] coordsys commenting --- fhd_core/beam_modeling/pyuvdata_beam_import.pro | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/fhd_core/beam_modeling/pyuvdata_beam_import.pro b/fhd_core/beam_modeling/pyuvdata_beam_import.pro index 30dd7e1dd..112992df4 100644 --- a/fhd_core/beam_modeling/pyuvdata_beam_import.pro +++ b/fhd_core/beam_modeling/pyuvdata_beam_import.pro @@ -6,14 +6,14 @@ MRD_HREAD, lun, primary_header, /silent coordsys = strtrim(sxpar(primary_header,'coordsys'),2) CASE coordsys OF - "az_za": BEGIN + "az_za": BEGIN ; Beam is in azimuth/zenith angle coordinates Jones_matrix = import_az_el_beam(obs, antenna_str, pyuvdata_filepath,$ za_arr=za_arr, az_arr=az_arr, psf_image_dim=psf_image_dim) END - "orthoslant_zenith": BEGIN + "orthoslant_zenith": BEGIN ; Beam is in orthoslant coordinates message,"COORDSYS type 'orthoslant_zenith' not yet supported!"" END - "healpix": BEGIN + "healpix": BEGIN ; Beam is in Healpix coordinates message,"COORDSYS type 'healpix' not yet supported!" END ELSE: message,"Unsupported COORDSYS found in pyuvdata file!" From f4c722969fd4b4bcefdd30ba6dd71be787e145ba Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Mon, 3 Apr 2023 14:38:40 -0700 Subject: [PATCH 07/16] remove more unneeded code --- fhd_core/beam_modeling/import_az_el_beam.pro | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/fhd_core/beam_modeling/import_az_el_beam.pro b/fhd_core/beam_modeling/import_az_el_beam.pro index 11c9a934f..67788de1b 100644 --- a/fhd_core/beam_modeling/import_az_el_beam.pro +++ b/fhd_core/beam_modeling/import_az_el_beam.pro @@ -59,9 +59,7 @@ ENDFOR zenith_i=where(theta_arr EQ 0,n_zenith) -horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) -;horizon_mask=fltarr(psf_image_dim,psf_image_dim)+1 -;IF n_horizon_test GT 0 THEN horizon_mask[horizon_test]=0 +horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) Jones_matrix=Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin) FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO $ Jones_matrix[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(psf_image_dim,psf_image_dim)) From 9f8609753112d5e0f906a0ca30cc44a4fbf3ebc8 Mon Sep 17 00:00:00 2001 From: Nichole Barry Date: Wed, 3 Aug 2022 23:08:03 +1000 Subject: [PATCH 08/16] param transfer defaults --- .../beam_modeling/fhd_struct_init_antenna.pro | 26 +++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/fhd_core/beam_modeling/fhd_struct_init_antenna.pro b/fhd_core/beam_modeling/fhd_struct_init_antenna.pro index e9f242c19..f5e73a318 100644 --- a/fhd_core/beam_modeling/fhd_struct_init_antenna.pro +++ b/fhd_core/beam_modeling/fhd_struct_init_antenna.pro @@ -4,8 +4,8 @@ FUNCTION fhd_struct_init_antenna,obs,beam_model_version=beam_model_version,$ psf_dim=psf_dim,psf_max_dim=psf_max_dim,beam_offset_time=beam_offset_time,$ inst_tile_ptr=inst_tile_ptr,ra_arr=ra_arr,dec_arr=dec_arr,fractional_size=fractional_size,$ kernel_window=kernel_window,beam_per_baseline=beam_per_baseline,$ - beam_gaussian_decomp=beam_gaussian_decomp,conserve_memory=conserve_memory,$ - import_pyuvdata_beam_filepath=import_pyuvdata_beam_filepath,$ + beam_gaussian_decomp=beam_gaussian_decomp,beam_gauss_param_transfer=beam_gauss_param_transfer,$ + conserve_memory=conserve_memory,import_pyuvdata_beam_filepath=import_pyuvdata_beam_filepath,$ use_psf_resolution=use_psf_resolution,_Extra=extra t0=Systime(1) @@ -23,6 +23,28 @@ if keyword_set(beam_gaussian_decomp) and keyword_set(kernel_window) then begin kernel_window=0 endif +;Default the parameter transfer if not set +if keyword_set(beam_gauss_param_transfer) then begin + ;Default to instrumental beam if not set + if (beam_gauss_param_transfer EQ 1) then beam_gauss_param_transfer = 'decomp' + + ;Set transfer to the instrumental beam or gaussian beam. Currently only available for the MWA + if (beam_gauss_param_transfer EQ 'decomp') or (beam_gauss_param_transfer EQ 'gauss') then begin + if instrument EQ 'mwa' then begin + pointing_num = mwa_get_pointing_number(obs,string=1) + beam_gauss_param_transfer = filepath(instrument + '_decomp_params_pointing' + pointing_num + '.sav',$ + root=rootdir('FHD'),subdir='instrument_config') + endif else begin + message, 'Gaussian decomposition parameter defaults not currently set for non-MWA instruments' + endelse + endif + + ;Match the psf size to the parameters from the file + psf_transfer = getvar_savefile(beam_gauss_param_transfer,'psf') + psf_resolution = psf_transfer.resolution + psf_dim = psf_transfer.dim +endif + n_tiles=obs.n_tile n_freq=obs.n_freq n_pol=obs.n_pol From cc717f8e030f68620153036a0b1b7f7b724bd487 Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Wed, 17 May 2023 14:11:47 -0700 Subject: [PATCH 09/16] remove extra keywords --- fhd_core/polarization/fhd_struct_init_jones.pro | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fhd_core/polarization/fhd_struct_init_jones.pro b/fhd_core/polarization/fhd_struct_init_jones.pro index 2f9121b55..05dd9f208 100644 --- a/fhd_core/polarization/fhd_struct_init_jones.pro +++ b/fhd_core/polarization/fhd_struct_init_jones.pro @@ -52,7 +52,7 @@ apply_astrometry, obs, dec_arr=lat, x_arr=x_t, y_arr=y_t, /xy2ad, /ignore_refrac obs_temp = fhd_struct_update_obs(obs, beam_nfreq_avg=obs.n_freq) ; Use only one average Jones matrix, not one per frequency antenna=fhd_struct_init_antenna(obs_temp,beam_model_version=beam_model_version,use_psf_resolution=1.,$ - psf_dim=obs.dimension,psf_intermediate_res=1.,psf_image_resolution=1.,timing=t_ant,_Extra=extra) + psf_dim=obs.dimension,psf_intermediate_res=1.,psf_image_resolution=1.,timing=t_ant,import_pyuvdata_beam_filepath=import_pyuvdata_beam_filepath) Jones=rotate_jones_matrix(obs, antenna[0]) ; Calculate the normalization From 5a9d452197df39dce97b1581ed74ddf6500d8a46 Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Wed, 17 May 2023 17:25:44 -0700 Subject: [PATCH 10/16] get keyword from extra structure --- fhd_core/polarization/fhd_struct_init_jones.pro | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fhd_core/polarization/fhd_struct_init_jones.pro b/fhd_core/polarization/fhd_struct_init_jones.pro index 05dd9f208..737bebb15 100644 --- a/fhd_core/polarization/fhd_struct_init_jones.pro +++ b/fhd_core/polarization/fhd_struct_init_jones.pro @@ -52,7 +52,7 @@ apply_astrometry, obs, dec_arr=lat, x_arr=x_t, y_arr=y_t, /xy2ad, /ignore_refrac obs_temp = fhd_struct_update_obs(obs, beam_nfreq_avg=obs.n_freq) ; Use only one average Jones matrix, not one per frequency antenna=fhd_struct_init_antenna(obs_temp,beam_model_version=beam_model_version,use_psf_resolution=1.,$ - psf_dim=obs.dimension,psf_intermediate_res=1.,psf_image_resolution=1.,timing=t_ant,import_pyuvdata_beam_filepath=import_pyuvdata_beam_filepath) + psf_dim=obs.dimension,psf_intermediate_res=1.,psf_image_resolution=1.,timing=t_ant,import_pyuvdata_beam_filepath=extra.import_pyuvdata_beam_filepath) Jones=rotate_jones_matrix(obs, antenna[0]) ; Calculate the normalization From 8a94394e828d9633e2552fbca6b3bc5cf51a6506 Mon Sep 17 00:00:00 2001 From: Nichole Barry Date: Thu, 25 May 2023 14:39:36 +1000 Subject: [PATCH 11/16] Reduce mem by only using selected pix --- fhd_core/beam_modeling/import_az_el_beam.pro | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/fhd_core/beam_modeling/import_az_el_beam.pro b/fhd_core/beam_modeling/import_az_el_beam.pro index 67788de1b..4026f0ef2 100644 --- a/fhd_core/beam_modeling/import_az_el_beam.pro +++ b/fhd_core/beam_modeling/import_az_el_beam.pro @@ -62,7 +62,7 @@ zenith_i=where(theta_arr EQ 0,n_zenith) horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) Jones_matrix=Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin) FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO $ - Jones_matrix[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(psf_image_dim,psf_image_dim)) + Jones_matrix[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(n_pix)) interp_res=obs.degpix angle_slice_i0=Uniq(phi_arr) @@ -83,10 +83,8 @@ az_ang_inst=Atan(yvals_instrument[pix_use],xvals_instrument[pix_use])*!Radeg+180 FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO BEGIN Jmat_use=Reform(*Jmat_interp[p_i,p_j,freq_i],n_zen_slice,n_ang_slice) Expand,Jmat_use,n_zen_ang,n_az_ang,Jmat_single - (*Jones_matrix[p_i,p_j,freq_i])[pix_use]=Interpolate(Jmat_single,zen_ang_inst/interp_res,az_ang_inst/interp_res,cubic=-0.5) - debug_point=1 + (*Jones_matrix[p_i,p_j,freq_i])=Interpolate(Jmat_single,zen_ang_inst/interp_res,az_ang_inst/interp_res,cubic=-0.5) ENDFOR -debug_point=1 RETURN, Jones_matrix END From 87c9db9530f9530e6615260652af6157b57f7be4 Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Mon, 26 Jun 2023 14:13:05 -0700 Subject: [PATCH 12/16] undo reduce mem changes --- fhd_core/beam_modeling/import_az_el_beam.pro | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/fhd_core/beam_modeling/import_az_el_beam.pro b/fhd_core/beam_modeling/import_az_el_beam.pro index 4026f0ef2..67788de1b 100644 --- a/fhd_core/beam_modeling/import_az_el_beam.pro +++ b/fhd_core/beam_modeling/import_az_el_beam.pro @@ -62,7 +62,7 @@ zenith_i=where(theta_arr EQ 0,n_zenith) horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) Jones_matrix=Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin) FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO $ - Jones_matrix[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(n_pix)) + Jones_matrix[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(psf_image_dim,psf_image_dim)) interp_res=obs.degpix angle_slice_i0=Uniq(phi_arr) @@ -83,8 +83,10 @@ az_ang_inst=Atan(yvals_instrument[pix_use],xvals_instrument[pix_use])*!Radeg+180 FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO BEGIN Jmat_use=Reform(*Jmat_interp[p_i,p_j,freq_i],n_zen_slice,n_ang_slice) Expand,Jmat_use,n_zen_ang,n_az_ang,Jmat_single - (*Jones_matrix[p_i,p_j,freq_i])=Interpolate(Jmat_single,zen_ang_inst/interp_res,az_ang_inst/interp_res,cubic=-0.5) + (*Jones_matrix[p_i,p_j,freq_i])[pix_use]=Interpolate(Jmat_single,zen_ang_inst/interp_res,az_ang_inst/interp_res,cubic=-0.5) + debug_point=1 ENDFOR +debug_point=1 RETURN, Jones_matrix END From 89fa16d8a3d5ab80d35e5b6abf692917d6aabb53 Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Tue, 27 Jun 2023 15:07:55 -0700 Subject: [PATCH 13/16] transpose Jones matrix --- fhd_core/beam_modeling/import_az_el_beam.pro | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fhd_core/beam_modeling/import_az_el_beam.pro b/fhd_core/beam_modeling/import_az_el_beam.pro index 67788de1b..308d41f06 100644 --- a/fhd_core/beam_modeling/import_az_el_beam.pro +++ b/fhd_core/beam_modeling/import_az_el_beam.pro @@ -42,7 +42,7 @@ FOR freq_i=0,n_freq-1 DO BEGIN FOR b_i=0,n_basis-1 DO BEGIN Jmat1 = Reform(Jmat0[*, *, freq_i, f_i, 0, b_i, 0] + icomp*Jmat0[*, *, freq_i, f_i, 0, b_i, 1]) Jmat1 = Transpose(Jmat1) - Jmat_arr[freq_i,f_i,b_i,*] = Reform(Jmat1, n_ang) + Jmat_arr[freq_i,b_i,f_i,*] = Reform(Jmat1, n_ang) ENDFOR ENDFOR ENDFOR From ae2c5b5164c44be8d691bceb07d0da3b514c0c50 Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Mon, 31 Jul 2023 16:05:51 -0700 Subject: [PATCH 14/16] flip Jones elements --- fhd_core/beam_modeling/import_az_el_beam.pro | 1 + 1 file changed, 1 insertion(+) diff --git a/fhd_core/beam_modeling/import_az_el_beam.pro b/fhd_core/beam_modeling/import_az_el_beam.pro index 308d41f06..3da42425e 100644 --- a/fhd_core/beam_modeling/import_az_el_beam.pro +++ b/fhd_core/beam_modeling/import_az_el_beam.pro @@ -43,6 +43,7 @@ FOR freq_i=0,n_freq-1 DO BEGIN Jmat1 = Reform(Jmat0[*, *, freq_i, f_i, 0, b_i, 0] + icomp*Jmat0[*, *, freq_i, f_i, 0, b_i, 1]) Jmat1 = Transpose(Jmat1) Jmat_arr[freq_i,b_i,f_i,*] = Reform(Jmat1, n_ang) + if b_i eq 1 then Jmat_arr[freq_i,b_i,f_i,*] *= -1 ;flip to align with MWA beam conventions ENDFOR ENDFOR ENDFOR From 6afaef3b3a77d7f98ac717d9b8d1f883f244fb1e Mon Sep 17 00:00:00 2001 From: Nichole Barry Date: Mon, 28 Aug 2023 18:59:00 +1000 Subject: [PATCH 15/16] add mem reduction, fix further pix selection bug --- fhd_core/beam_modeling/import_az_el_beam.pro | 6 ++---- fhd_core/beam_modeling/pyuvdata_beam_import.pro | 12 ------------ 2 files changed, 2 insertions(+), 16 deletions(-) diff --git a/fhd_core/beam_modeling/import_az_el_beam.pro b/fhd_core/beam_modeling/import_az_el_beam.pro index 3da42425e..39c109a30 100644 --- a/fhd_core/beam_modeling/import_az_el_beam.pro +++ b/fhd_core/beam_modeling/import_az_el_beam.pro @@ -63,7 +63,7 @@ zenith_i=where(theta_arr EQ 0,n_zenith) horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix) Jones_matrix=Ptrarr(n_ant_pol,n_ant_pol,nfreq_bin) FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO $ - Jones_matrix[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(psf_image_dim,psf_image_dim)) + Jones_matrix[p_i,p_j,freq_i]=Ptr_new(Dcomplexarr(n_pix)) interp_res=obs.degpix angle_slice_i0=Uniq(phi_arr) @@ -84,10 +84,8 @@ az_ang_inst=Atan(yvals_instrument[pix_use],xvals_instrument[pix_use])*!Radeg+180 FOR p_i=0,n_ant_pol-1 DO FOR p_j=0,n_ant_pol-1 DO FOR freq_i=0L,nfreq_bin-1 DO BEGIN Jmat_use=Reform(*Jmat_interp[p_i,p_j,freq_i],n_zen_slice,n_ang_slice) Expand,Jmat_use,n_zen_ang,n_az_ang,Jmat_single - (*Jones_matrix[p_i,p_j,freq_i])[pix_use]=Interpolate(Jmat_single,zen_ang_inst/interp_res,az_ang_inst/interp_res,cubic=-0.5) - debug_point=1 + (*Jones_matrix[p_i,p_j,freq_i])=Interpolate(Jmat_single,zen_ang_inst/interp_res,az_ang_inst/interp_res,cubic=-0.5) ENDFOR -debug_point=1 RETURN, Jones_matrix END diff --git a/fhd_core/beam_modeling/pyuvdata_beam_import.pro b/fhd_core/beam_modeling/pyuvdata_beam_import.pro index 112992df4..052fbac36 100644 --- a/fhd_core/beam_modeling/pyuvdata_beam_import.pro +++ b/fhd_core/beam_modeling/pyuvdata_beam_import.pro @@ -19,17 +19,5 @@ CASE coordsys OF ELSE: message,"Unsupported COORDSYS found in pyuvdata file!" ENDCASE -for pol1=0,(size(Jones_matrix))[1]-1 do begin - for pol2=0,(size(Jones_matrix))[2]-1 do begin - if (size(Jones_matrix))[0] gt 2 then begin - for ind=0,(size(Jones_matrix))[3]-1 do begin - *Jones_matrix[pol1, pol2, ind] = (*Jones_matrix[pol1, pol2, ind])[pix_use] - endfor - endif else begin - *Jones_matrix[pol1, pol2] = (*Jones_matrix[pol1, pol2])[pix_use] - endelse - endfor -endfor - RETURN, Jones_matrix END From 07bc4f83921210aa1894fee63b26045d3e22f63c Mon Sep 17 00:00:00 2001 From: Ruby Byrne Date: Wed, 24 Jul 2024 00:00:47 +0000 Subject: [PATCH 16/16] update skyh5 heading to pyradiosky standard --- fhd_core/setup_metadata/load_skyh5_source_catalog.pro | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fhd_core/setup_metadata/load_skyh5_source_catalog.pro b/fhd_core/setup_metadata/load_skyh5_source_catalog.pro index a2bec2c76..b23447125 100644 --- a/fhd_core/setup_metadata/load_skyh5_source_catalog.pro +++ b/fhd_core/setup_metadata/load_skyh5_source_catalog.pro @@ -5,8 +5,8 @@ FUNCTION load_skyh5_source_catalog, catalog_path, $ if ~keyword_set(skyh5_name_stokes) then skyh5_name_stokes='/Data/stokes' if ~keyword_set(skyh5_name_nsources) then skyh5_name_nsources='/Header/Ncomponents' - if ~keyword_set(skyh5_name_ra) then skyh5_name_ra= '/Header/lon' - if ~keyword_set(skyh5_name_dec) then skyh5_name_dec='/Header/lat' + if ~keyword_set(skyh5_name_ra) then skyh5_name_ra= '/Header/skycoord/ra' + if ~keyword_set(skyh5_name_dec) then skyh5_name_dec='/Header/skycoord/dec' if ~keyword_set(skyh5_name_freqs) then skyh5_name_freqs='/Header/reference_frequency' if ~keyword_set(skyh5_name_alphas) then skyh5_name_alphas='/Header/spectral_index'