Read pyuvdata-formatted beamfits beam models

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 <br />
   -*Default*: 0.05, or 0.01 if `allow_sidelobe_sources` set <br />
 
-**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 `<instrument>_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. <br />
+**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 `<instrument>_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. <br />
   -*MWA range*: 0, 1 (or anything else captured in the `else` statement), 2 <br />
   -*PAPER range*: 1 (or anything else captured in the `else` statement), 2 <br />
   -*HERA range*: 2 (or anything else captured in the `else` statement) <br />
@@ -49,6 +49,9 @@ This is a work in progress; please add keywords as you find them in alphabetical
   -*Default*: 1 <br />
   -*eor_wrapper_defaults*: 1 <br />
 
+**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. <br />
+  -*Default*: not set <br />
+
 **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. <br />
   -*Example*: `inst_tile_ptr = PTRARR(2,/allocate)`<br />
   `*inst_tile_ptr[0] = [0,1,2,3]`<br />

fhd_core/beam_modeling/fhd_struct_init_antenna.pro

@@ -1,27 +1,22 @@
 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
+    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
   IF conserve_memory GT 1E6 THEN mem_thresh=conserve_memory ELSE mem_thresh=1E8 ;in bytes
 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.'
@@ -61,8 +56,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 +96,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 +119,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 +142,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 +190,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 +224,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)

fhd_core/beam_modeling/import_az_el_beam.pro

@@ -0,0 +1,91 @@
+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)
+
+; 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')
+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,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
+
+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
+    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
+
+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))
+
+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)
+
+; 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
+    (*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
+
+RETURN, Jones_matrix
+END

fhd_core/beam_modeling/pyuvdata_beam_import.pro

@@ -0,0 +1,23 @@
+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  ; 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  ; Beam is in orthoslant coordinates
+        message,"COORDSYS type 'orthoslant_zenith' not yet supported!""
+    END
+    "healpix": BEGIN  ; Beam is in Healpix coordinates
+        message,"COORDSYS type 'healpix' not yet supported!"
+    END
+    ELSE: message,"Unsupported COORDSYS found in pyuvdata file!"
+ENDCASE
+
+RETURN, Jones_matrix
+END

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,import_pyuvdata_beam_filepath=extra.import_pyuvdata_beam_filepath)
 Jones=rotate_jones_matrix(obs, antenna[0])
 
 ; Calculate the normalization

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'
 

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

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