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]`
diff --git a/fhd_core/beam_modeling/fhd_struct_init_antenna.pro b/fhd_core/beam_modeling/fhd_struct_init_antenna.pro index a2f0b3d5a..f5e73a318 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 + 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,15 +15,8 @@ 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.' @@ -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) 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..39c109a30 --- /dev/null +++ b/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 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..052fbac36 --- /dev/null +++ b/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 diff --git a/fhd_core/polarization/fhd_struct_init_jones.pro b/fhd_core/polarization/fhd_struct_init_jones.pro index 30f0bac08..737bebb15 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,import_pyuvdata_beam_filepath=extra.import_pyuvdata_beam_filepath) Jones=rotate_jones_matrix(obs, antenna[0]) ; Calculate the normalization 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' 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