EveryBeam MWA integration tests

Goes without saying, but the test here rely on WODEN being compiled with Everybeam support. First import some code, and up some observational parameters.

import os
from subprocess import call
from astropy.io import fits
import numpy as np
from astropy.table import Column, Table
from astropy.time import Time
from astropy.coordinates import SkyCoord
from astropy.coordinates import EarthLocation
from astropy import units as u
import matplotlib.pyplot as plt
import numpy.testing as npt
from astropy.constants import c
from astropy.wcs import WCS
from wodenpy.primary_beam.use_everybeam import run_everybeam, run_everybeam_over_threads
import erfa
import mwa_hyperbeam
from wodenpy.use_libwoden.create_woden_struct_classes import Woden_Struct_Classes
from ctypes import c_double
from astropy.time import TimeDelta
from wodenpy.array_layout.precession import RTS_Precess_LST_Lat_to_J2000

import sys
# sys.path.append('../../scripts/')
# from run_woden import main as run_woden

##A bunch of test and plotting code lives here so we can use it in multile notebooks
from eb_testing_code import create_WCS, plot_jones_on_sky, plot_everybeam_on_sky, make_sky_models, read_uvfits, convert_inst_to_stokes, test_stokes_recovery, getFDF, findpeaks, test_RM_recovery, make_RM_skymodel

C = c.to('m/s').value

# ra0 = 0.0
# dec0 = -26.7
MWA_LAT=-26.703319405555554
MWA_LONG=116.67081523611111

##this date has highlighted issues with parallactic rotation when it's not
##working correctly for me, so use this one in testing
date = "2013-05-16T10:48:48"


##Assume that the OSKAR telescope is near the MWA??
mwa_location = EarthLocation(lat=MWA_LAT*u.deg, 
                            lon=MWA_LONG*u.deg,
                            height=377.827)
observing_time = Time(date, scale='utc', location=mwa_location)

##Grab the LST
LST = observing_time.sidereal_time('mean')
LST_deg = LST.value*15
print(f"LST: {LST_deg} deg")#, RA: {ra0}")

print(observing_time.mjd*86400)
LST: 153.18482172499105 deg
4875418128.0

What does the beam look like on the sky?

First up, let’s plot the hyperbeam MWA_FEE beam model on the sky to see what we’re expecting. WODEN likes the beam to be in IAU ordering, where the first (often called x) polarisation is aligned north-south. We also want this the have parallactic rotation, so switch that one. We’ll use hyperbeam first, which we already use and have tested extensively within WODEN.

Note that we precess the array and LST back to J2000; this is how WODEN works internally to make everything happen in the J2000 frame. I’m doing it here explicitly to check the beam works in J2000.

##Set up a grid of RA, Dec points

ra0 = LST_deg
dec0 = MWA_LAT

mjd_current = observing_time.mjd


lst_J2000, latitude_J2000 = RTS_Precess_LST_Lat_to_J2000(np.radians(LST_deg),
                                                         np.radians(MWA_LAT),
                                                         mjd_current)


nside=64
radec_reso=120/nside

header, wcs = create_WCS(ra0, dec0, nside, radec_reso)
x_mesh, y_mesh = np.meshgrid(np.arange(nside), np.arange(nside))
ras, decs = wcs.all_pix2world(x_mesh, y_mesh, 0)

ras = np.radians(ras.flatten())
decs = np.radians(decs.flatten())

##Then use erfa to convert these values into azs, els
# has = np.radians(LST_deg) - ras
has = lst_J2000 - ras

##use this erfa function to convert to azimuth and elevation
##were using degrees, but erfa uses rads, so convert here
az_grid, els = erfa.hd2ae(has, decs, latitude_J2000)
za_grid = np.pi/2 - els

##mwa hyperbeam
beam = mwa_hyperbeam.FEEBeam(hdf5_file=os.environ["MWA_FEE_HDF5"])

freq = 180e+6

parallactic_rotation = False
IAU_order = True

jones = beam.calc_jones_array(az_grid, za_grid, freq, [0]*16, [1]*16,
                              norm_to_zenith=True,
                              latitude_rad=latitude_J2000, iau_order=IAU_order)

all_gx_hyp = jones[:,0]
all_Dx_hyp = jones[:,1]
all_Dy_hyp = jones[:,2]
all_gy_hyp = jones[:,3]

all_gx_hyp.shape = (nside, nside)
all_Dx_hyp.shape = (nside, nside)
all_Dy_hyp.shape = (nside, nside)
all_gy_hyp.shape = (nside, nside)

plot_jones_on_sky(all_gx_hyp, all_Dx_hyp, all_Dy_hyp, all_gy_hyp, wcs, title="hyperbeam with parallactic rotation")
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
../_images/612eb9fb0db96a4c5b8d5e1592dc62cf255d9ce7e97f9616a7b129f5c3bdaa45.png

Do a sanity check that we have polarisations in the order that we want. North dipoles are most sensitive towards the East-West, so if we subtract Y from X, we should see a positive signal towards the East-West, and negative towards the North-South.

diff = np.abs(all_gx_hyp) - np.abs(all_gy_hyp)
im = plt.imshow(diff, origin='lower')
plt.colorbar(im)
plt.xlabel("RA")
plt.ylabel("Dec")
plt.show()
../_images/ce9949197e41bdb3a90edece91cbf02921cdd100eb69a9c50cd50f7ac429e21f.png

Aight, cool, we have a sensible beam patter benchmark to compare against.

EveryBeam MWA

Let’s now do the same with everybeam MWA and see how it compares. We’ll use the wodenpy.use_everybeam.run_everybeam function to do this, which is what the wodenpy sky model reading code uses. However, it’s a fairly slow function, so we’ll run it in parallel over a number of cores. To do this, we’ll use the helper function wodenpy.use_everybeam.run_everybeam_over_threads. This uses the same technique that run_woden.py uses to parallelise the sky model/beam calculations.

Unlike other EveryBeam telescopes, we directly hack into the EveryBeam function that takes az,za as inputs. All other telescopes we use don’t have a direct az,za interface (that I’ve been able to find). But by inputting az,za directly, we can explicitly control precession effects ourselves. Inputting RA,Dec makes EveryBeam use casacore to do the precession. Note of course we are inputting ra,dec into run_everybeam_over_threads, but that function does the same erfa.ha2ae as shown in the cell above to get az,za before passing to EveryBeam.

##Example MWA MS from the EveryBeam package
ms_path="MWA-single-timeslot.ms"
coeff_path=os.environ["MWA_FEE_HDF5"]
j2000_latitudes = [latitude_J2000]
j2000_lsts = [lst_J2000]
freqs = [freq]
parallactic_rotate=True
iau_order=True
##zenith pointing
mwa_dipole_delays = np.zeros(16)

##How many cores to split the work over
num_threads=8


all_jones = run_everybeam_over_threads(num_threads, ras, decs, freqs,
                            mwa_coeff_path=coeff_path,
                            mwa_dipole_delays=mwa_dipole_delays,
                            j2000_latitudes=j2000_latitudes,
                            j2000_lsts=j2000_lsts,
                            iau_order=iau_order,
                            parallactic_rotate=parallactic_rotate)

beam_ind, time_ind, freq_ind = 0, 0, 0

all_gx_eb = all_jones[beam_ind, time_ind, freq_ind, :, 0, 0]
all_Dx_eb = all_jones[beam_ind, time_ind, freq_ind, :, 0, 1]
all_Dy_eb = all_jones[beam_ind, time_ind, freq_ind, :, 1, 0]
all_gy_eb = all_jones[beam_ind, time_ind, freq_ind, :, 1, 1]

all_gx_eb.shape = (nside, nside)
all_Dx_eb.shape = (nside, nside)
all_Dy_eb.shape = (nside, nside)
all_gy_eb.shape = (nside, nside)

diff_gx = all_gx_eb - all_gx_hyp
diff_Dx = all_Dx_eb - all_Dx_hyp
diff_Dy = all_Dy_eb - all_Dy_hyp
diff_gy = all_gy_eb - all_gy_hyp

diff_gx[int(nside//2), int(nside//2)] = np.nan
diff_Dx[int(nside//2), int(nside//2)] = np.nan
diff_Dy[int(nside//2), int(nside//2)] = np.nan
diff_gy[int(nside//2), int(nside//2)] = np.nan

plot_jones_on_sky(all_gx_eb, all_Dx_eb, all_Dy_eb, all_gy_eb, wcs, 'Everybeam with parallactic rotation')
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
../_images/953c4bc791ba68879ad3d9eac4315a5e1932171717812e39597ee80c82d61332.png
diff = np.abs(all_gx_eb) - np.abs(all_gy_eb)
im = plt.imshow(diff, origin='lower')
plt.colorbar(im)
plt.xlabel("RA")
plt.ylabel("Dec")
plt.show()
../_images/ce9949197e41bdb3a90edece91cbf02921cdd100eb69a9c50cd50f7ac429e21f.png

Alright let’s see how the hyperbeam and everybeam MWA beams compare!

diff_gx = all_gx_eb - all_gx_hyp
diff_Dx = all_Dx_eb - all_Dx_hyp
diff_Dy = all_Dy_eb - all_Dy_hyp
diff_gy = all_gy_eb - all_gy_hyp

plot_jones_on_sky(diff_gx, diff_Dx, diff_Dy, diff_gy, wcs, 'Difference between everybeam and hyperbeam')
../_images/0aeb3f98ea2c8c603a25a7dd0d433677aa710065eac4257a378a3e6a47a422e3.png

Lovely. Things are looking extremely similar, which is great.

Stokes recovery

Now let’s check that we can recover the Stokes parameters from the beam. We’ll just test a single point source at zenith, where the gains should essentially be 1.0 and leakage zero. We’ll run all Stokes of I, Q, U, V through the full simulator, recombine the linear polarisations, and check we get back what we put in.

First, make our sky models. We’ll make 4 Stokes I models, where each has either I, Q, U, V = 1, and all other Stokes are zero. To make that happen for Q/U, make the component a list type so explicitly set one or the other to zero.

make_sky_models(ra0, dec0)

Keep the simulations tiny by making a 3 ant array

np.random.seed(234987)

##make a random array layout. Source is at phase centre, so the visibilities
##should all be full real and just be set by the flux in the sky
num_antennas = 3
east = np.random.uniform(-1000, 1000, num_antennas)
north = np.random.uniform(-1000, 1000, num_antennas)
height = np.zeros(num_antennas)

array = np.empty((num_antennas, 3))
array[:,0] = east
array[:,1] = north
array[:,2] = height

array_name = "eg_array.txt"
np.savetxt(array_name, array)

Run the 4 different Stokes sky models through hyperbeam first as a point of comparison. We’re avoiding using metafits here to be as apples to apples with EveryBeam as possible.

##parameters for the simulation; setting the date sets an LST of about 0 deg
##to match the phase centre

freq_reso = 1e+6
low_freq = 180e+6
high_freq = 210e+6
num_freq_chans = int((high_freq - low_freq) / freq_reso)

primary_beam = "MWA_FEE"

print("DATE", date)

for pol in ['I', 'Q', 'U', 'V']:

    uvfits_name = f"stokes{pol}_{primary_beam}"
    cat_name = f'{pol}_source.fits'

    cmd = 'run_woden.py '

    cmd += f'--ra0={ra0} '
    cmd += f'--dec0={dec0} '
    cmd += f'--array_layout={array_name} '
    cmd += f'--date={date} '
    cmd += f'--output_uvfits_prepend={uvfits_name} '
    cmd += f'--primary_beam={primary_beam} '
    cmd += f'--freq_res={freq_reso} '
    cmd += f'--band_nums=1 '
    cmd += f'--num_time_steps=1 '
    cmd += f'--IAU_order '
    cmd += f'--cat_filename={cat_name} '
    cmd += f'--lowest_channel_freq={low_freq} '
    cmd += f'--num_freq_channels={num_freq_chans} '
    cmd += f'--time_res=2 '
    cmd += f'--MWA_FEE_delays=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] '
    cmd += f'--num_threads=1 '
    
    call(cmd, shell=True)
      
DATE 2013-05-16T10:48:48
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:10:56 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:10:56 - INFO - Input arguments after parsing:
                             	Phase centre: 153.18482, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 180.000 MHz
                             	Channel frequency resolution: 1000.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 30
                             	Start date: 2013-05-16T10:48:48
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 3 antenna positions from array layout file: eg_array.txt
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:10:56 - INFO - Obs epoch initial LST was 153.1889998083 deg
2025-10-23 10:10:56 - INFO - Setting initial J2000 LST to 153.0310851527 deg
2025-10-23 10:10:56 - INFO - Setting initial mjd to 56428.4505671296
2025-10-23 10:10:56 - INFO - After precession initial latitude of the array is -26.6349592375 deg
2025-10-23 10:10:56 - INFO - Precessing array layout to J2000
2025-10-23 10:10:56 - INFO - Using MWA primary beam via hyperbeam with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:10:56 - INFO - Doing the initial mapping of sky model
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
2025-10-23 10:10:57 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:10:57 - INFO - Have read in 1 components
2025-10-23 10:10:57 - INFO - After cropping there are 1 components
2025-10-23 10:10:57 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:10:57 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:10:57 - INFO - Reading set 0 sky models
2025-10-23 10:10:57 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:10:57 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:10:57 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-23 10:10:58 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:10:58 - INFO - Full run took 0:00:01.839038
2025-10-23 10:10:58 - INFO - WODEN is done. Closing the log. S'later
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:10:59 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:10:59 - INFO - Input arguments after parsing:
                             	Phase centre: 153.18482, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 180.000 MHz
                             	Channel frequency resolution: 1000.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 30
                             	Start date: 2013-05-16T10:48:48
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 3 antenna positions from array layout file: eg_array.txt
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:11:00 - INFO - Obs epoch initial LST was 153.1889998083 deg
2025-10-23 10:11:00 - INFO - Setting initial J2000 LST to 153.0310851527 deg
2025-10-23 10:11:00 - INFO - Setting initial mjd to 56428.4505671296
2025-10-23 10:11:00 - INFO - After precession initial latitude of the array is -26.6349592375 deg
2025-10-23 10:11:00 - INFO - Precessing array layout to J2000
2025-10-23 10:11:00 - INFO - Using MWA primary beam via hyperbeam with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:11:00 - INFO - Doing the initial mapping of sky model
2025-10-23 10:11:00 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:11:00 - INFO - Have read in 1 components
2025-10-23 10:11:00 - INFO - After cropping there are 1 components
2025-10-23 10:11:00 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:11:00 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:11:00 - INFO - Reading set 0 sky models
2025-10-23 10:11:00 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:11:00 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:11:01 - INFO - Sky set 0 has returned from the GPU after 0.2 seconds
2025-10-23 10:11:01 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:11:01 - INFO - Full run took 0:00:01.804909
2025-10-23 10:11:01 - INFO - WODEN is done. Closing the log. S'later
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:11:02 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:11:02 - INFO - Input arguments after parsing:
                             	Phase centre: 153.18482, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 180.000 MHz
                             	Channel frequency resolution: 1000.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 30
                             	Start date: 2013-05-16T10:48:48
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 3 antenna positions from array layout file: eg_array.txt
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:11:03 - INFO - Obs epoch initial LST was 153.1889998083 deg
2025-10-23 10:11:03 - INFO - Setting initial J2000 LST to 153.0310851527 deg
2025-10-23 10:11:03 - INFO - Setting initial mjd to 56428.4505671296
2025-10-23 10:11:03 - INFO - After precession initial latitude of the array is -26.6349592375 deg
2025-10-23 10:11:03 - INFO - Precessing array layout to J2000
2025-10-23 10:11:03 - INFO - Using MWA primary beam via hyperbeam with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:11:03 - INFO - Doing the initial mapping of sky model
2025-10-23 10:11:03 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:11:03 - INFO - Have read in 1 components
2025-10-23 10:11:03 - INFO - After cropping there are 1 components
2025-10-23 10:11:03 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:11:03 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:11:03 - INFO - Reading set 0 sky models
2025-10-23 10:11:03 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:11:04 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:11:04 - INFO - Sky set 0 has returned from the GPU after 0.2 seconds
2025-10-23 10:11:04 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:11:04 - INFO - Full run took 0:00:01.877820
2025-10-23 10:11:04 - INFO - WODEN is done. Closing the log. S'later
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:11:06 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:11:06 - INFO - Input arguments after parsing:
                             	Phase centre: 153.18482, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 180.000 MHz
                             	Channel frequency resolution: 1000.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 30
                             	Start date: 2013-05-16T10:48:48
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 3 antenna positions from array layout file: eg_array.txt
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:11:07 - INFO - Obs epoch initial LST was 153.1889998083 deg
2025-10-23 10:11:07 - INFO - Setting initial J2000 LST to 153.0310851527 deg
2025-10-23 10:11:07 - INFO - Setting initial mjd to 56428.4505671296
2025-10-23 10:11:07 - INFO - After precession initial latitude of the array is -26.6349592375 deg
2025-10-23 10:11:07 - INFO - Precessing array layout to J2000
2025-10-23 10:11:07 - INFO - Using MWA primary beam via hyperbeam with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:11:07 - INFO - Doing the initial mapping of sky model
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
2025-10-23 10:11:07 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:11:07 - INFO - Have read in 1 components
2025-10-23 10:11:07 - INFO - After cropping there are 1 components
2025-10-23 10:11:07 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:11:07 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:11:07 - INFO - Reading set 0 sky models
2025-10-23 10:11:07 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:11:07 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:11:08 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-23 10:11:08 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:11:08 - INFO - Full run took 0:00:02.147754
2025-10-23 10:11:08 - INFO - WODEN is done. Closing the log. S'later

Tests are defined in eb_testing_code.py, which just asserts that whichever Stokes param we set to one comes out as one, and everything else is zero. Make the absolute tolerance of the test something to toggle, because real beam has leakage, and we don’t have the beam at exactly the phase centre for hyperbeam. The Stokes recovery is only really legit when the XX and YY beams are the same, so it’s good at phase centre and zenith.

Do the actual tests for MWA FEE beam

for pol in ['I', 'Q', 'U', 'V']:
    test_stokes_recovery(pol, 'MWA_FEE', atol=5e-3)
    
    # #Uncomment to print out example values
    # uvfits_name = f"stokes{pol}_MWA_FEE"
    # XX, XY, YX, YY = read_uvfits(f'{uvfits_name}_band01.uvfits')
    # ##pick a random baseline to plot, they should all be the sam
    # baseline = 0

    # recover_I, recover_Q, recover_U, recover_V = convert_inst_to_stokes(XX[baseline], XY[baseline], YX[baseline], YY[baseline])
    # print(f"{pol}, {recover_I[0].real:.2f}, {recover_Q[0].real:.2f}, {recover_U[0].real:.2f}, {recover_V[0].real:.2f}")
Testing Stokes I
Stokes I passed
Testing Stokes Q
Stokes Q passed
Testing Stokes U
Stokes U passed
Testing Stokes V
Stokes V passed

Right, try and do the same thing with EveryBeam

##parameters for the simulation; setting the date sets an LST of about 0 deg
##to match the phase centre

freq_reso = 1e+6
low_freq = 180e+6
high_freq = 210e+6
num_freq_chans = int((high_freq - low_freq) / freq_reso)

primary_beam = "everybeam_MWA"

for pol in ['I', 'Q', 'U', 'V']:

    uvfits_name = f"stokes{pol}_{primary_beam}"
    cat_name = f'{pol}_source.fits'
    
    cmd = 'run_woden.py '

    cmd += f'--ra0={ra0} '
    cmd += f'--dec0={dec0} '
    cmd += f'--array_layout={array_name} '
    cmd += f'--date={date} '
    cmd += f'--output_uvfits_prepend={uvfits_name} '
    cmd += f'--primary_beam={primary_beam} '
    cmd += f'--freq_res={freq_reso} '
    cmd += f'--band_nums=1 '
    cmd += f'--num_time_steps=1 '
    cmd += f'--IAU_order '
    cmd += f'--cat_filename={cat_name} '
    cmd += f'--lowest_channel_freq={low_freq} '
    cmd += f'--num_freq_channels={num_freq_chans} '
    cmd += f'--time_res=2 '
    cmd += f'--MWA_FEE_delays=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] '
    cmd += f'--num_threads=1 '
    
    call(cmd, shell=True)
    # print(cmd)
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:26:38 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:26:38 - INFO - Input arguments after parsing:
                             	Phase centre: 153.18482, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 180.000 MHz
                             	Channel frequency resolution: 1000.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 30
                             	Start date: 2013-05-16T10:48:48
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 3 antenna positions from array layout file: eg_array.txt
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:26:38 - INFO - Obs epoch initial LST was 153.1889998083 deg
2025-10-23 10:26:38 - INFO - Setting initial J2000 LST to 153.0310851527 deg
2025-10-23 10:26:38 - INFO - Setting initial mjd to 56428.4505671296
2025-10-23 10:26:38 - INFO - After precession initial latitude of the array is -26.6349592375 deg
2025-10-23 10:26:38 - INFO - Precessing array layout to J2000
2025-10-23 10:26:38 - INFO - Using MWA primary beam via EveryBeam MWA with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:26:38 - INFO - Doing the initial mapping of sky model
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
2025-10-23 10:26:38 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:26:38 - INFO - Have read in 1 components
2025-10-23 10:26:38 - INFO - After cropping there are 1 components
2025-10-23 10:26:38 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:26:39 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:26:39 - INFO - Reading set 0 sky models
2025-10-23 10:26:39 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:26:39 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:26:39 - INFO - libwoden: Loading EveryBeam MWA tile beam
2025-10-23 10:26:39 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-23 10:26:39 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:26:39 - INFO - Full run took 0:00:01.910199
2025-10-23 10:26:39 - INFO - WODEN is done. Closing the log. S'later
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:26:41 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:26:41 - INFO - Input arguments after parsing:
                             	Phase centre: 153.18482, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 180.000 MHz
                             	Channel frequency resolution: 1000.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 30
                             	Start date: 2013-05-16T10:48:48
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 3 antenna positions from array layout file: eg_array.txt
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:26:42 - INFO - Obs epoch initial LST was 153.1889998083 deg
2025-10-23 10:26:42 - INFO - Setting initial J2000 LST to 153.0310851527 deg
2025-10-23 10:26:42 - INFO - Setting initial mjd to 56428.4505671296
2025-10-23 10:26:42 - INFO - After precession initial latitude of the array is -26.6349592375 deg
2025-10-23 10:26:42 - INFO - Precessing array layout to J2000
2025-10-23 10:26:42 - INFO - Using MWA primary beam via EveryBeam MWA with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:26:42 - INFO - Doing the initial mapping of sky model
2025-10-23 10:26:42 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:26:42 - INFO - Have read in 1 components
2025-10-23 10:26:42 - INFO - After cropping there are 1 components
2025-10-23 10:26:42 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:26:42 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:26:42 - INFO - Reading set 0 sky models
2025-10-23 10:26:42 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:26:42 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:26:42 - INFO - libwoden: Loading EveryBeam MWA tile beam
2025-10-23 10:26:43 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-23 10:26:43 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:26:43 - INFO - Full run took 0:00:01.967022
2025-10-23 10:26:43 - INFO - WODEN is done. Closing the log. S'later
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:26:44 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:26:44 - INFO - Input arguments after parsing:
                             	Phase centre: 153.18482, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 180.000 MHz
                             	Channel frequency resolution: 1000.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 30
                             	Start date: 2013-05-16T10:48:48
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 3 antenna positions from array layout file: eg_array.txt
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:26:45 - INFO - Obs epoch initial LST was 153.1889998083 deg
2025-10-23 10:26:45 - INFO - Setting initial J2000 LST to 153.0310851527 deg
2025-10-23 10:26:45 - INFO - Setting initial mjd to 56428.4505671296
2025-10-23 10:26:45 - INFO - After precession initial latitude of the array is -26.6349592375 deg
2025-10-23 10:26:45 - INFO - Precessing array layout to J2000
2025-10-23 10:26:45 - INFO - Using MWA primary beam via EveryBeam MWA with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:26:45 - INFO - Doing the initial mapping of sky model
2025-10-23 10:26:45 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:26:45 - INFO - Have read in 1 components
2025-10-23 10:26:45 - INFO - After cropping there are 1 components
2025-10-23 10:26:45 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:26:45 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:26:45 - INFO - Reading set 0 sky models
2025-10-23 10:26:45 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:26:46 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:26:46 - INFO - libwoden: Loading EveryBeam MWA tile beam
2025-10-23 10:26:46 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-23 10:26:46 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:26:46 - INFO - Full run took 0:00:01.919606
2025-10-23 10:26:46 - INFO - WODEN is done. Closing the log. S'later
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:26:48 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:26:48 - INFO - Input arguments after parsing:
                             	Phase centre: 153.18482, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 180.000 MHz
                             	Channel frequency resolution: 1000.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 30
                             	Start date: 2013-05-16T10:48:48
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 3 antenna positions from array layout file: eg_array.txt
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:26:48 - INFO - Obs epoch initial LST was 153.1889998083 deg
2025-10-23 10:26:48 - INFO - Setting initial J2000 LST to 153.0310851527 deg
2025-10-23 10:26:48 - INFO - Setting initial mjd to 56428.4505671296
2025-10-23 10:26:48 - INFO - After precession initial latitude of the array is -26.6349592375 deg
2025-10-23 10:26:48 - INFO - Precessing array layout to J2000
2025-10-23 10:26:48 - INFO - Using MWA primary beam via EveryBeam MWA with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:26:48 - INFO - Doing the initial mapping of sky model
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
2025-10-23 10:26:48 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:26:48 - INFO - Have read in 1 components
2025-10-23 10:26:48 - INFO - After cropping there are 1 components
2025-10-23 10:26:48 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:26:49 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:26:49 - INFO - Reading set 0 sky models
2025-10-23 10:26:49 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:26:49 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:26:49 - INFO - libwoden: Loading EveryBeam MWA tile beam
2025-10-23 10:26:49 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-23 10:26:49 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:26:49 - INFO - Full run took 0:00:01.869466
2025-10-23 10:26:49 - INFO - WODEN is done. Closing the log. S'later
for pol in ['I', 'Q', 'U', 'V']:
    test_stokes_recovery(pol, 'everybeam_MWA', atol=5e-3)
    
    # #Uncomment to print out example values
    # uvfits_name = f"stokes{pol}_everybeam_MWA"
    # XX, XY, YX, YY = read_uvfits(f'{uvfits_name}_band01.uvfits')
    # ##pick a random baseline to plot, they should all be the sam
    # baseline = 0

    # recover_I, recover_Q, recover_U, recover_V = convert_inst_to_stokes(XX[baseline], XY[baseline], YX[baseline], YY[baseline])
    # print(f"{pol}, {recover_I[0].real:.2f}, {recover_Q[0].real:.2f}, {recover_U[0].real:.2f}, {recover_V[0].real:.2f}")
Testing Stokes I
Stokes I passed
Testing Stokes Q
Stokes Q passed
Testing Stokes U
Stokes U passed
Testing Stokes V
Stokes V passed

Now test an RM recovery

Check we recover the correct RM and sign for a linearly polarised source. First of all, make the sky model.

phi_RM, pol_frac = make_RM_skymodel(ra0, dec0)

Run that through WODEN using hyperbeam to check the test is working.

freq_reso = 0.1e+6
low_freq = 180e+6
high_freq = 210e+6
num_freq_chans = int((high_freq - low_freq) / freq_reso)


primary_beam = "MWA_FEE"

uvfits_name = f"rm_source_{primary_beam}"
cat_name = 'RM_source.fits'

cmd = 'run_woden.py '

cmd += f'--ra0={ra0} '
cmd += f'--dec0={dec0} '
cmd += f'--array_layout={array_name} '
cmd += f'--date={date} '
cmd += f'--output_uvfits_prepend={uvfits_name} '
cmd += f'--primary_beam={primary_beam} '
cmd += f'--freq_res={freq_reso} '
cmd += f'--band_nums=1 '
cmd += f'--num_time_steps=1 '
cmd += f'--IAU_order '
cmd += f'--cat_filename={cat_name} '
cmd += f'--lowest_channel_freq={low_freq} '
cmd += f'--num_freq_channels={num_freq_chans} '
cmd += f'--time_res=2 '
cmd += f'--MWA_FEE_delays=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] '
cmd += f'--num_threads=1 '


call(cmd, shell=True)    
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:27:02 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:27:02 - INFO - Input arguments after parsing:
                             	Phase centre: 153.18482, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 180.000 MHz
                             	Channel frequency resolution: 100.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 300
                             	Start date: 2013-05-16T10:48:48
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 3 antenna positions from array layout file: eg_array.txt
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:27:03 - INFO - Obs epoch initial LST was 153.1889998083 deg
2025-10-23 10:27:03 - INFO - Setting initial J2000 LST to 153.0310851527 deg
2025-10-23 10:27:03 - INFO - Setting initial mjd to 56428.4505671296
2025-10-23 10:27:03 - INFO - After precession initial latitude of the array is -26.6349592375 deg
2025-10-23 10:27:03 - INFO - Precessing array layout to J2000
2025-10-23 10:27:03 - INFO - Using MWA primary beam via hyperbeam with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:27:03 - INFO - Doing the initial mapping of sky model
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
2025-10-23 10:27:03 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:27:03 - INFO - Have read in 1 components
2025-10-23 10:27:03 - INFO - After cropping there are 1 components
2025-10-23 10:27:03 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:27:03 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:27:03 - INFO - Reading set 0 sky models
2025-10-23 10:27:03 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:27:04 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:27:04 - INFO - Sky set 0 has returned from the GPU after 0.2 seconds
2025-10-23 10:27:04 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:27:04 - INFO - Full run took 0:00:01.850418
2025-10-23 10:27:04 - INFO - WODEN is done. Closing the log. S'later
0

Have a look and check our visis make sense

XX, XY, YX, YY = read_uvfits('rm_source_MWA_FEE_band01.uvfits')
baseline = 0    
recover_I, recover_Q, recover_U, recover_V = convert_inst_to_stokes(XX[baseline], XY[baseline], YX[baseline], YY[baseline])

    
freqs = np.arange(low_freq, high_freq, freq_reso)

fig, axs = plt.subplots(1, 1, figsize=(8, 6))

axs.plot(freqs / 1e+6, recover_I.real, 'k-', label='Recovered I')
axs.plot(freqs / 1e+6, recover_Q.real, 'r', label='Recovered Q')
axs.plot(freqs / 1e+6, recover_U.real, 'b', label='Recovered U')
axs.plot(freqs / 1e+6, recover_V.real, 'g', label='Recovered V')

axs.legend()

axs.set_xlabel('Frequency (MHz)')
axs.set_ylabel('Flux (Jy)')

plt.show()

# fig, axs = plt.subplots(1, 1, figsize=(8, 6))
# axs.plot(freqs / 1e+6, recover_I.real, 'k-', label='Polarisation Fraction')
../_images/dffbec841399e84c0a5fa43da1d0fc87f24d3223bb52ff7a7fc7f03354e60d23.png

Now define a test, using a bunch of Emil’s RM synthesis code.

test_RM_recovery(uvfits_name, phi_RM, pol_frac, freqs)
../_images/9f3ea772e5cfa968beffc5448aa6b0d131864d4e0ab69bbe4154405e06bd4f6c.png
Recovered RM: 50.0 Expected RM: 50
Recovered Pol. Fraction: 0.49999502 Expected Pol Fraction 0.5

Loverly, works a treat. Now let’s do that with EveryBeam.

freq_reso = 0.1e+6
low_freq = 180e+6
high_freq = 210e+6
num_freq_chans = int((high_freq - low_freq) / freq_reso)

primary_beam = "everybeam_MWA"
uvfits_name = f"rm_source_{primary_beam}"
cat_name = 'RM_source.fits'

cmd = 'run_woden.py '

cmd += f'--ra0={ra0} '
cmd += f'--dec0={dec0} '
cmd += f'--array_layout={array_name} '
cmd += f'--date={date} '
cmd += f'--output_uvfits_prepend={uvfits_name} '
cmd += f'--primary_beam={primary_beam} '
cmd += f'--freq_res={freq_reso} '
cmd += f'--band_nums=1 '
cmd += f'--num_time_steps=1 '
cmd += f'--IAU_order '
cmd += f'--cat_filename={cat_name} '
cmd += f'--lowest_channel_freq={low_freq} '
cmd += f'--num_freq_channels={num_freq_chans} '
cmd += f'--MWA_FEE_delays=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] '
cmd += f'--time_res=2 '
cmd += f'--num_threads=1 '

call(cmd, shell=True)
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:27:18 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:27:18 - INFO - Input arguments after parsing:
                             	Phase centre: 153.18482, -26.70332 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 180.000 MHz
                             	Channel frequency resolution: 100.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 300
                             	Start date: 2013-05-16T10:48:48
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 3 antenna positions from array layout file: eg_array.txt
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:27:18 - INFO - Obs epoch initial LST was 153.1889998083 deg
2025-10-23 10:27:18 - INFO - Setting initial J2000 LST to 153.0310851527 deg
2025-10-23 10:27:18 - INFO - Setting initial mjd to 56428.4505671296
2025-10-23 10:27:18 - INFO - After precession initial latitude of the array is -26.6349592375 deg
2025-10-23 10:27:18 - INFO - Precessing array layout to J2000
2025-10-23 10:27:18 - INFO - Using MWA primary beam via EveryBeam MWA with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:27:18 - INFO - Doing the initial mapping of sky model
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
2025-10-23 10:27:18 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:27:18 - INFO - Have read in 1 components
2025-10-23 10:27:18 - INFO - After cropping there are 1 components
2025-10-23 10:27:18 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:27:19 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:27:19 - INFO - Reading set 0 sky models
2025-10-23 10:27:19 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:27:19 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:27:19 - INFO - libwoden: Loading EveryBeam MWA tile beam
2025-10-23 10:27:19 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-23 10:27:20 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:27:20 - INFO - Full run took 0:00:01.973844
2025-10-23 10:27:20 - INFO - WODEN is done. Closing the log. S'later
0
XX, XY, YX, YY = read_uvfits('rm_source_everybeam_MWA_band01.uvfits')
baseline = 0    
recover_I, recover_Q, recover_U, recover_V = convert_inst_to_stokes(XX[baseline], XY[baseline], YX[baseline], YY[baseline])

    
freqs = np.arange(low_freq, high_freq, freq_reso)

fig, axs = plt.subplots(1, 1, figsize=(8, 6))

axs.plot(freqs / 1e+6, recover_I.real, 'k-', label='Recovered I')
axs.plot(freqs / 1e+6, recover_Q.real, 'r', label='Recovered Q')
axs.plot(freqs / 1e+6, recover_U.real, 'b', label='Recovered U')
axs.plot(freqs / 1e+6, recover_V.real, 'g', label='Recovered V')

axs.legend()

axs.set_xlabel('Frequency (MHz)')
axs.set_ylabel('Flux (Jy)')

plt.show()
../_images/cb2f30fa0fe914217fc0e0b5aabd003f2792227715256609142068356b8d0974.png

Looks good!

test_RM_recovery(uvfits_name, phi_RM, pol_frac, freqs)
../_images/9f3ea772e5cfa968beffc5448aa6b0d131864d4e0ab69bbe4154405e06bd4f6c.png
Recovered RM: 50.0 Expected RM: 50
Recovered Pol. Fraction: 0.49999502 Expected Pol Fraction 0.5

Yup works woot

Try off zenith pointings

Now we’ll test whether we can do off-zenith pointings. We’ll also test that we can run multiple time steps and frequencies.

ra0 = 60.0
dec0 = -30.0
date = "2015-09-10T19:44:15"
observing_time = Time(date, scale='utc', location=mwa_location)

##Setup a grid of RA/Dec on the sky
nside=50
radec_reso = 120/nside

header, wcs = create_WCS(ra0, dec0, nside, radec_reso)
x_mesh, y_mesh = np.meshgrid(np.arange(nside), np.arange(nside))
ras, decs = wcs.all_pix2world(x_mesh, y_mesh, 0)

ras = np.radians(ras.flatten())
decs = np.radians(decs.flatten())

num_comps = len(ras)

##Do a single freqs and time and beam
all_freqs = np.array([100e+6, 200e+6])
all_times = np.array([observing_time, observing_time + TimeDelta(2*3600.0, format='sec')])

lsts = []
latitudes = []

for obs_time in all_times:
    
    LST = obs_time.sidereal_time('mean')
    lst_current = np.radians(LST.value*15)
    
    lst_J2000, latitude_J2000 = RTS_Precess_LST_Lat_to_J2000(
                                lst_current,
                                np.radians(MWA_LAT),
                                obs_time.mjd)
    
    lsts.append(lst_J2000)
    latitudes.append(latitude_J2000)

coeff_path= os.environ["MWA_FEE_HDF5"]
delays = [0,2,4,6,0,2,4,6,0,2,4,6,0,2,4,6]


num_threads = 8

all_jones = run_everybeam_over_threads(num_threads, ras, decs, all_freqs,
                            mwa_coeff_path=coeff_path,
                            mwa_dipole_delays=delays,
                            j2000_latitudes=latitudes,
                            j2000_lsts=lsts,
                            iau_order=True,
                            parallactic_rotate=True)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
beam_ind = 0

for time_ind in range(len(all_times)):
    for freq_ind in range(len(all_freqs)):
        this_gx = all_jones[beam_ind, time_ind, freq_ind, :, 0, 0]
        this_Dx = all_jones[beam_ind, time_ind, freq_ind, :, 0, 1]
        this_Dy = all_jones[beam_ind, time_ind, freq_ind, :, 1, 0]
        this_gy = all_jones[beam_ind, time_ind, freq_ind, :, 1, 1]

        this_gx.shape = (nside, nside)
        this_Dx.shape = (nside, nside)
        this_Dy.shape = (nside, nside)
        this_gy.shape = (nside, nside)

        plot_jones_on_sky(this_gx, this_Dx, this_Dy, this_gy, wcs, title="everybeam with parallactic rotation")
../_images/1e86551a934557e946115b58e0a3dc57e02d3c90ef0368cbbeaad4366c4dc629.png ../_images/96a763372538db66038f218fa9565d1942ee18b775ccdf7b78bd4d6517b45fbf.png ../_images/9c128b20ace2305e2ffdd2a317b1bd288fb95da109a3b34f180ecd3b6f11938a.png ../_images/8f09dab7f85be30a9e2e171e177f85cc8e9d1c6e9bb62a810ff90ee0f3538ae1.png

Double check what it should look like using hyperdrive, where we can explicitly put the delays that we want in.

##use this erfa function to convert to azimuth and elevation
##were using degrees, but erfa uses rads, so convert here

for this_time in all_times:
    for freq in all_freqs:

        has = np.radians(this_time.sidereal_time('mean').value*15.0) - ras

        az_grid, els = erfa.hd2ae(has, decs, np.radians(MWA_LAT))

        za_grid = np.pi/2 - els

        beam = mwa_hyperbeam.FEEBeam(hdf5_file=os.environ["MWA_FEE_HDF5"])

        parallactic_rotation = True
        delays = [0,2,4,6,0,2,4,6,0,2,4,6,0,2,4,6]

        jones = beam.calc_jones_array(az_grid, za_grid, freq, delays, [1]*16, True, np.radians(MWA_LAT), parallactic_rotation)

        all_gx = jones[:,0]
        all_Dx = jones[:,1]
        all_Dy = jones[:,2]
        all_gy = jones[:,3]

        all_gx.shape = (nside, nside)
        all_Dx.shape = (nside, nside)
        all_Dy.shape = (nside, nside)
        all_gy.shape = (nside, nside)

        plot_jones_on_sky(all_gx, all_Dx, all_Dy, all_gy, wcs, title="hyperbeam with parallactic rotation")
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
../_images/0b1b77dd38c56f431cde162c1ce624e7fc78b3502726111b13c8dbc3abf26cb9.png
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
../_images/b17b8064cd20f2b4532ead904c341705f7f322414ccc78697773ee826b0dc02b.png
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
../_images/e005ef8851748c58499ad466c36f84e1c9a1847596703f983c406329afb34163.png
/home/jack-line/software/WODEN_dev/woden_dev/lib/python3.12/site-packages/erfa/core.py:18088: RuntimeWarning: invalid value encountered in hd2ae
  az, el = ufunc.hd2ae(ha, dec, phi)
../_images/ab00c6b017866d841a12eaa031b504bad6c3c0a48e726030e4fbdad4a8b8c4be.png

Interpolated beam

Can we use the interpolated Daniel Ung .h5 file for highband? Let’s try it.

from astropy.table import Column, Table
import numpy as np
import os

sing_ra = 15.0
sing_dec = -40.0

ras = [sing_ra]
decs = [sing_dec]

num_comps = len(ras)

c_ids = Column(data=np.array(['source']*num_comps), name='UNQ_SOURCE_ID', dtype='|S20')
c_names = Column(data=np.array([f'source_C{i:04d}' for i in range(num_comps)]), name='NAME', dtype='|S20')

##Component position
c_ras = Column(data=ras, name='RA')
c_decs = Column(data=decs, name='DEC')

##This says we have a point source
c_comp_types = Column(data=np.array(['P']*num_comps, dtype='|S1'), name="COMP_TYPE", dtype='|S1')
##This says we have a Stokes I power-law SED
c_mod_types = Column(data=np.array(['pl']*num_comps, dtype='|S3'), name="MOD_TYPE", dtype='|S3')

##Set everything as a flat spectrum 1 Jy stokes I source. That way we can image it and see the beam pattern
c_stokes_I_ref = Column(data=np.ones(num_comps), name='NORM_COMP_PL')
c_stokes_I_SI = Column(data=np.zeros(num_comps), name='ALPHA_PL')


cols = [c_ids, c_names, c_ras, c_decs, c_comp_types, c_mod_types, c_stokes_I_ref, c_stokes_I_SI]

main_table = Table(cols)

single_cat = 'single_point.fits'

main_table.write(single_cat, format='fits', overwrite=True)


ra0=0.0
dec0=-26.7
date="2024-07-21T20:13:00"
primary_beam="everybeam_MWA"
num_times=1
num_freq_chans=390
low_freq=167e+6
ms_path="MWA-single-timeslot.ms"
uvfits_name="interp_run_woden"
freq_reso=80e+3

args = []

interp_h5 = os.environ['MWA_FEE_HDF5_INTERP']

cmd = 'run_woden.py '

cmd += f'--ra0={sing_ra} '
cmd += f'--dec0={-sing_dec} '
cmd += f'--metafits_filename=../../examples/metafits/1125949472_metafits.fits '
cmd += f'--output_uvfits_prepend={uvfits_name} '
cmd += f'--primary_beam={primary_beam} '
cmd += f'--freq_res={freq_reso} '
cmd += f'--band_nums=1 '
cmd += f'--num_time_steps={num_times} '
cmd += f'--IAU_order '
cmd += f'--cat_filename={single_cat} '
cmd += f'--lowest_channel_freq={low_freq} '
cmd += f'--num_freq_channels={num_freq_chans} '
cmd += f'--time_res=2 '
cmd += f'--num_threads=1 '
cmd += f'--hdf5_beam_path={interp_h5} '

call(cmd, shell=True)
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:33:10 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:33:10 - INFO - Input arguments after parsing:
                             	Phase centre: 15.00000, 40.00000 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 167.000 MHz
                             	Channel frequency resolution: 80.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 390
                             	Start date: 2015-09-10T19:44:15
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 128 antenna positions from metafits file: ../../examples/metafits/1125949472_metafits.fits
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:33:11 - INFO - Obs epoch initial LST was 42.2620215148 deg
2025-10-23 10:33:11 - INFO - Setting initial J2000 LST to 42.0911874758 deg
2025-10-23 10:33:11 - INFO - Setting initial mjd to 57275.8224074073
2025-10-23 10:33:11 - INFO - After precession initial latitude of the array is -26.7665041369 deg
2025-10-23 10:33:11 - INFO - Precessing array layout to J2000
2025-10-23 10:33:11 - INFO - Using MWA primary beam via EveryBeam MWA with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/MWA_embedded_element_pattern_rev2_interp_167_197MHz.h5
                             	delays: [0 2 4 6 0 2 4 6 0 2 4 6 0 2 4 6]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:33:11 - INFO - Doing the initial mapping of sky model
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
2025-10-23 10:33:11 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:33:11 - INFO - Have read in 1 components
2025-10-23 10:33:11 - INFO - After cropping there are 1 components
2025-10-23 10:33:11 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:33:12 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:33:12 - INFO - Reading set 0 sky models
2025-10-23 10:33:12 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:33:13 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:33:14 - INFO - libwoden: Loading EveryBeam MWA tile beam
2025-10-23 10:33:46 - INFO - Sky set 0 has returned from the GPU after 32.2 seconds
2025-10-23 10:33:48 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:33:50 - INFO - Full run took 0:00:39.770120
2025-10-23 10:33:50 - INFO - WODEN is done. Closing the log. S'later
0
coarse_h5 = os.environ['MWA_FEE_HDF5']

uvfits_name = 'coarse_run_woden'

cmd = 'run_woden.py '

cmd += f'--ra0={sing_ra} '
cmd += f'--dec0={-sing_dec} '
cmd += f'--metafits_filename=../../examples/metafits/1125949472_metafits.fits '
cmd += f'--output_uvfits_prepend={uvfits_name} '
cmd += f'--primary_beam={primary_beam} '
cmd += f'--freq_res={freq_reso} '
cmd += f'--band_nums=1 '
cmd += f'--num_time_steps={num_times} '
cmd += f'--IAU_order '
cmd += f'--station_id=0 '
cmd += f'--cat_filename={single_cat} '
cmd += f'--lowest_channel_freq={low_freq} '
cmd += f'--num_freq_channels={num_freq_chans} '
cmd += f'--time_res=2 '
cmd += f'--num_threads=1 '
cmd += f'--hdf5_beam_path={coarse_h5} '

call(cmd, shell=True)
/home/jack-line/software/WODEN_dev/woden_dev/bin/run_woden.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('wodenpy==2.6.0')
2025-10-23 10:34:41 - INFO - 
                                              )  (              )  
                                  (  (     ( /(  )\ )        ( /(  
                                  )\))(   ')\())(()/(   (    )\()) 
                                 ((_)()\ )((_)\  /(_))  )\  ((_)\  
                                 _(())\_)() ((_)(_))_  ((_)  _((_) 
                                 \ \((_)/ // _ \ |   \ | __|| \| | 
                                  \ \/\/ /| (_) || |) || _| | .` | 
                                   \_/\_/  \___/ |___/ |___||_|\_| 
                                   
                                 You are using wodenpy version/git hash: c8b7005
                                 
2025-10-23 10:34:41 - INFO - Input arguments after parsing:
                             	Phase centre: 15.00000, 40.00000 deg
                             	Array central latitude: -26.703 deg
                             	Array central longitude: 116.671 deg
                             	Array central height: 377.827 m
                             	Lowest channel frequency: 167.000 MHz
                             	Channel frequency resolution: 80.000 kHz
                             	Coarse band bandwidth: 1.280 MHz
                             	Num channels per coarse band: 390
                             	Start date: 2015-09-10T19:44:15
                             	Time resolution: 2.0 (s)
                             	Num time steps: 1
                             	Have read 128 antenna positions from metafits file: ../../examples/metafits/1125949472_metafits.fits
                             	Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-23 10:34:43 - INFO - Obs epoch initial LST was 42.2620215148 deg
2025-10-23 10:34:43 - INFO - Setting initial J2000 LST to 42.0911874758 deg
2025-10-23 10:34:43 - INFO - Setting initial mjd to 57275.8224074073
2025-10-23 10:34:43 - INFO - After precession initial latitude of the array is -26.7665041369 deg
2025-10-23 10:34:43 - INFO - Precessing array layout to J2000
2025-10-23 10:34:43 - INFO - Using MWA primary beam via EveryBeam MWA with the following parameters:
                             	hdf5 file: /home/jack-line/software/mwa_beam_files/mwa_full_embedded_element_pattern.h5
                             	delays: [0 2 4 6 0 2 4 6 0 2 4 6 0 2 4 6]
                             	setting all dipole amplitudes to 1.0
                             	will not flag any dipoles
2025-10-23 10:34:43 - INFO - Doing the initial mapping of sky model
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
2025-10-23 10:34:43 - INFO - Sky model mapping took 0.0 mins
2025-10-23 10:34:43 - INFO - Have read in 1 components
2025-10-23 10:34:43 - INFO - After cropping there are 1 components
2025-10-23 10:34:43 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-23 10:34:43 - INFO - Running in serial on GPU.
                             Will read sky model using 1 threads
                             There are 1 sets of sky models to run
2025-10-23 10:34:43 - INFO - Reading set 0 sky models
2025-10-23 10:34:43 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-23 10:34:44 - INFO - Sending Sky set 0 to the GPU
2025-10-23 10:34:45 - INFO - libwoden: Loading EveryBeam MWA tile beam
2025-10-23 10:34:46 - INFO - Sky set 0 has returned from the GPU after 1.8 seconds
2025-10-23 10:34:49 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-23 10:34:50 - INFO - Full run took 0:00:08.903036
2025-10-23 10:34:50 - INFO - WODEN is done. Closing the log. S'later
0
xx_interp, xy_interp, yx_interp, yy_interp = read_uvfits('interp_run_woden_band01.uvfits')
xx_coarse, xy_coarse, yx_coarse, yy_coarse = read_uvfits('coarse_run_woden_band01.uvfits')

freqs = np.arange(low_freq, low_freq + num_freq_chans*freq_reso, freq_reso)

xx_abs = np.abs(xx_coarse[0])
plt.plot(freqs/1e+6, xx_abs, label='Coarse h5 file')

xx_abs = np.abs(xx_interp[0])
plt.plot(freqs/1e+6, xx_abs, label='Interpolated h5 file')
plt.xlabel('Frequency (MHz)')
plt.ylabel('XX Amplitude')
plt.legend()
plt.show()
../_images/b1be143895803afd3039f8896594c8a591fefd3e14044cc86010b4f402254ecd.png

Yup, seems like we can indeed; the interpolated beam is slower however. I’m calling EveryBeam the exact same number of times, but EveryBeam is smart enough to cache repeated values, so the coarse h5 file runs faster.