EveryBeam OSKAR “MWA” integration tests¶
As WODEN has been heavily tested with the MWA beam, it makes sense to start by trying the replicate an MWA beam through everybeam OSKAR. We can do that by creating a measurement set with the MWA array layout and the tile configuration, and then just use the OSKAR dipole model to do the beam forming. This let’s us compare results with known tested hyperbeam MWA_FEE. Should be good enough to check we have the right sign conventions, and things look sensible on the sky.
To see how the OSKAR-MWA measurement set was created, see the WODEN/test_installation/everybeam/create_OSKAR-MWA_ms/make_oskar-mwa_beam.ipynb notebook.
Goes without saying, but the test here rely on everybeam being installed.
import sys
##This is to get hold of EveryBeam on my machine
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
os.environ["MWA_BEAM_FILE"] = "/home/jline/software/MWA_beam_files/mwa_full_embedded_element_pattern.h5"
C = c.to('m/s').value
# ra0 = 0.0
# dec0 = -26.7
MWA_LAT=-26.703319405555554
MWA_LONG=116.67081523611111
##pick a time/date that sticks our phase centre overhead
date = "2024-07-21T20:13:00"
##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
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
##Set up a grid of RA, Dec points
from os import environ
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=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())
##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(environ["MWA_FEE_HDF5"])
freq = 160e+6
parallactic_rotation = True
IAU_order = True
jones = beam.calc_jones_array(az_grid, za_grid, freq, [0]*16, [1]*16,
IAU_order, latitude_J2000, parallactic_rotation)
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)
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()
Next, let’s do the equivalent for the everybeam OSKAR “MWA” beam model and see how it compares. We’ll ask EveryBeam to do the parallactic rotation for us
##I made the measurement set using OSKAR with the MWA layout
ms_path="create_OSKAR-MWA_ms/OSKAR-MWA-layout.ms"
coeff_path=""
station_id = 0
station_ids = [0]
j2000_latitudes = [latitude_J2000]
j2000_lsts = [lst_J2000]
times = [observing_time]
freqs = [freq]
current_latitude = np.radians(MWA_LAT)
current_longitude = np.radians(MWA_LONG)
apply_beam_norms=True
parallactic_rotate=True
iau_order=True
element_response = "skala40_wave"
##How many cores to split the work over
num_threads=8
# telescope = load_OSKAR_telescope(ms_path)
all_jones = run_everybeam_over_threads(num_threads, ms_path,
coeff_path,
ras, decs,
np.radians(ra0), np.radians(dec0),
j2000_latitudes, j2000_lsts,
times, freqs, station_ids,
apply_beam_norms=apply_beam_norms,
iau_order=iau_order,
parallactic_rotate=parallactic_rotate,
element_response_model=element_response)
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)
plot_jones_on_sky(all_gx_eb, all_Dx_eb, all_Dy_eb, all_gy_eb, wcs, 'Everybeam OSKAR-MWA eb-rotate')
Thread 2 processing coords 626 to 939Thread 3 processing coords 939 to 1252Thread 0 processing coords 0 to 313Thread 1 processing coords 313 to 626Thread 4 processing coords 1252 to 1565Thread 5 processing coords 1565 to 1878Thread 7 processing coords 2191 to 2504
Thread 6 processing coords 1878 to 2191
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:00 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
Thread 6 finished
Thread 7 finished
Thread 2 finished
Thread 3 finished
Thread 4 finished
Thread 0 finished
Thread 1 finished
Thread 5 finished
OK. So it’s not identical, but the absolute power is very similar. The phases are different, but this is using whatever the skala40_wave dipole model is. Some kind of logarithmic dipole, so a different element to the bowtie dipoles used in the MWA FEE model. So we’ll say that we’re using things correctly. Just double check that we have the expected polarisation ordering.
Under IAU conventions, the X polarisation is the North-South Dipole. The NS dipole is most sensitive to the East-West, so if we subtract the Y polarisation from the X polarisation, we should see a positive signal towards the East-West, and negative towards the North-South.
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()
Stokes IQUV recovery¶
Now test that we can recover Stokes IQUV from a point source.
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. This just checks that all the beam is creating XX,XY,YX,YY visibilities that come out with the sign conventions we expect.
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)
Right, run the simulations.
##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_OSKAR"
for pol in ['I', 'Q', 'U', 'V']:
# for pol in ['I']:
uvfits_name = f"stokes{pol}_{primary_beam}"
cat_name = f'{pol}_source.fits'
##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
command += f'--time_res=2 --num_time_steps=1 --IAU_order '
command += f'--station_id=0 ' ##default behaviour is to use a different primary beam for each station. All stations are the same here so pick first one.
command += f'--beam_ms_path=create_OSKAR-MWA_ms/OSKAR-MWA-layout.ms '
command += '--num_threads=1'
call(command, shell=True)
# print(command)
Note I’ve deleted the outputs of the cell above, because EveryBeam spits out a bunch of error messages with the OSKAR beam warning about frequencies mismatching, and not having a central direction set. The outputs seems fine. It’s a TODO to silence these warnings.
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.
for pol in ['I', 'Q', 'U', 'V']:
# for pol in ['I']:
test_stokes_recovery(pol, 'everybeam_OSKAR', atol=1e-1)
##Uncomment to print out example values
uvfits_name = f"stokes{pol}_everybeam_OSKAR"
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
I, 0.93, 0.00, -0.00, -0.00
Testing Stokes Q
Stokes Q passed
Q, 0.00, 0.93, -0.00, -0.00
Testing Stokes U
Stokes U passed
U, -0.00, 0.00, 0.93, 0.00
Testing Stokes V
Stokes V passed
V, -0.00, 0.00, -0.00, 0.93
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
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_OSKAR"
uvfits_name = f"rm_source_{primary_beam}"
cat_name = 'RM_source.fits'
##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
command += f'--time_res=2 --num_time_steps=1 --IAU_order '
command += f' --station_id=0 '
command += f'--beam_ms_path=create_OSKAR-MWA_ms/OSKAR-MWA-layout.ms '
command += '--num_threads=1'
call(command, shell=True)
Note I’ve deleted the outputs of the cell above, because EveryBeam spits out a bunch of error messages with the OSKAR beam warning about frequencies mismatching, and not having a central direction set. The outputs seems fine. It’s a TODO to silence these warnings.
XX, XY, YX, YY = read_uvfits('rm_source_everybeam_OSKAR_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()
Righto, because we are normalising to beam centre, and also simulating the point source at beam centre, we get something lovely and perfect.
test_RM_recovery(uvfits_name, phi_RM, pol_frac, freqs, atol=1e-2)
Recovered RM: 50.0 Expected RM: 50
Recovered Pol. Fraction: 0.50002867 Expected Pol Fraction 0.5
Works loverly.
Test antenna locations¶
Now check to see if the antenna locations are read in correctly from a measurement set. We’ll run WODEN using the metafits file that was used to make the OSKAR MWA measurement set. Then compare that to running the with the same settings but using the OSKAR MWA measurement set directly.
freq_reso = 0.1e+6
low_freq = 180e+6
num_freq_chans = 1
primary_beam = "MWA_FEE"
uvfits_name = f"ant_locs_{primary_beam}"
cat_name = 'RM_source.fits'
##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
command += f'--metafits_filename=../../examples/metafits/1126115208_metafits.fits '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
command += f'--time_res=2 --num_time_steps=1 --IAU_order '
command += '--MWA_FEE_delays=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] '
command += '--num_threads=1'
call(command, 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.5.0')
2025-03-06 09:00:34 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-06 09:00:34 - INFO - Input arguments after parsing:
Phase centre: 0.00455, -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: 1
Start date: 2024-07-21T20:13:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 128 antenna positions from metafits file: ../../examples/metafits/1126115208_metafits.fits
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-03-06 09:00:34 - INFO - Obs epoch initial LST was 0.0087300922 deg
2025-03-06 09:00:34 - INFO - Setting initial J2000 LST to 359.6931723607 deg
2025-03-06 09:00:34 - INFO - Setting initial mjd to 60512.8423726853
2025-03-06 09:00:34 - INFO - After precession initial latitude of the array is -26.8398178265 deg
2025-03-06 09:00:34 - INFO - Precessing array layout to J2000
2025-03-06 09:00:34 - 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-03-06 09:00:34 - 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-03-06 09:00:34 - INFO - Sky model mapping took 0.0 mins
2025-03-06 09:00:34 - INFO - Have read in 1 components
2025-03-06 09:00:34 - INFO - After cropping there are 1 components
2025-03-06 09:00:34 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-06 09:00:35 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-03-06 09:00:35 - INFO - Reading set 0 sky models
2025-03-06 09:00:35 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 09:00:35 - INFO - Sending Sky set 0 to the GPU
2025-03-06 09:00:35 - INFO - Sky set 0 has returned from the GPU after 0.1 seconds
2025-03-06 09:00:35 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-03-06 09:00:35 - INFO - Full run took 0:00:01.517604
2025-03-06 09:00:35 - INFO - WODEN is done. Closing the log. S'later
0
primary_beam = "everybeam_OSKAR"
uvfits_name = f"ant_locs_{primary_beam}"
cat_name = 'RM_source.fits'
##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
command += f'--time_res=2 --num_time_steps=1 --IAU_order '
command += f' --station_id=0 '
command += f'--beam_ms_path=create_OSKAR-MWA_ms/OSKAR-MWA-layout.ms '
command += '--num_threads=1'
call(command, 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.5.0')
Successful readonly open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-layout.ms: 22 columns, 8128 rows
Successful readonly open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-layout.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-layout.ms::FIELD: 9 columns, 1 rows
Successful readonly open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-layout.ms/ANTENNA: 8 columns, 128 rows
2025-03-06 09:00:36 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-06 09:00:36 - INFO - Input arguments after parsing:
Phase centre: 0.00455, -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: 0.040 MHz
Num channels per coarse band: 1
Start date: 2024-07-21T20:13:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 128 antenna positions from measurement set: create_OSKAR-MWA_ms/OSKAR-MWA-layout.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-03-06 09:00:37 - INFO - Obs epoch initial LST was 0.0087300922 deg
2025-03-06 09:00:37 - INFO - Setting initial J2000 LST to 359.6931723607 deg
2025-03-06 09:00:37 - INFO - Setting initial mjd to 60512.8423726853
2025-03-06 09:00:37 - INFO - After precession initial latitude of the array is -26.8398178265 deg
2025-03-06 09:00:37 - INFO - Precessing array layout to J2000
2025-03-06 09:00:37 - INFO - Will run with EveryBeam OSKAR primary beam, based on this measurement set:
create_OSKAR-MWA_ms/OSKAR-MWA-layout.ms
Primary beam is pointed at RA,Dec = 0.0, -29.999999999999996 deg
2025-03-06 09:00:37 - 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-03-06 09:00:37 - INFO - Sky model mapping took 0.0 mins
2025-03-06 09:00:37 - INFO - Have read in 1 components
2025-03-06 09:00:37 - INFO - After cropping there are 1 components
2025-03-06 09:00:37 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-06 09:00:37 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-03-06 09:00:37 - INFO - Reading set 0 sky models
2025-03-06 09:00:37 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 09:00:37 - INFO - Finshed sky set 0 reading thread num 0 in 0.0 seconds
2025-03-06 09:00:38 - INFO - Sending Sky set 0 to the GPU
POINTING LOFAR beam 0.0 -30.0
Could not load dataset for frequency 180 MHz, using the nearest neighbor with frequency 160 MHz instead
2025-03-06 09:00:38 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-03-06 09:00:38 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-03-06 09:00:38 - INFO - Full run took 0:00:01.751037
2025-03-06 09:00:38 - INFO - WODEN is done. Closing the log. S'later
2025-03-05 23:00:39 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:39 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:39 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
0
Read in the XYZ positions from the resultant uvfits files, as well as the uvw coordinates.
def read_xyz_uvfits(uvfits_name):
with fits.open(uvfits_name) as hdu:
xyz = hdu[1].data['STABXYZ']
uu = hdu[0].data['UU']*C
vv = hdu[0].data['VV']*C
ww = hdu[0].data['WW']*C
return xyz, uu, vv, ww
xyz_meta, uu_meta, vv_meta, ww_meta = read_xyz_uvfits('ant_locs_MWA_FEE_band01.uvfits')
xyz_ms, uu_ms, vv_ms, ww_ms = read_xyz_uvfits('ant_locs_everybeam_OSKAR_band01.uvfits')
First up compare the XYZ coords
def plot(xyz, ax1, ax2, marker, label):
ax1.plot(xyz[:,0], xyz[:,1], marker, mfc='none', label=label)
ax2.plot(xyz[:,0], xyz[:,2], marker, mfc='none', label=label)
fig, axs = plt.subplots(2, 2, figsize=(12, 12), layout='constrained')
plot(xyz_meta, axs[0,0], axs[1,0], 'o', 'metafits')
plot(xyz_ms, axs[0,0], axs[1,0], 'x', 'ms')
plot(xyz_meta, axs[0,1], axs[1,1], 'o', 'metafits')
plot(xyz_ms, axs[0,1], axs[1,1], 'x', 'ms')
axs[0,1].set_xlim(425, 500)
axs[0,1].set_ylim(-200, 0)
axs[1,1].set_xlim(430, 480)
axs[1,1].set_ylim(0, 100)
for ax in axs.flatten():
ax.legend()
axs[0,0].set_xlabel('X (m)')
axs[0,0].set_xlabel('Y (m)')
axs[0,1].set_xlabel('X (m)')
axs[0,1].set_xlabel('Y (m)')
axs[1,0].set_xlabel('X (m)')
axs[1,0].set_xlabel('Z (m)')
axs[1,1].set_xlabel('X (m)')
axs[1,1].set_xlabel('Z (m)')
plt.show()
np.allclose(xyz_meta, xyz_ms, atol=1e-3, rtol=1e-6)
True
Now compare the uv coords.
fig, axs = plt.subplots(1, 2, figsize=(12, 6), layout='constrained')
for ax in axs:
ax.plot(uu_meta, vv_meta, 'C0o', mfc='none', label='metafits')
# ax.plot(-uu_meta, -vv_meta, 'C0o', mfc='none')
ax.plot(uu_ms, vv_ms, 'C1x', mfc='none', label='ms')
# ax.plot(-uu_ms, -vv_ms, 'C1x', mfc='none')
ax.set_xlabel('UU (m)')
ax.set_ylabel('VV (m)')
ax.legend()
axs[1].set_xlim(-100, 100)
axs[1].set_ylim(-100, 100)
plt.show()
print(f"u coords max, mean diff {np.max(np.abs(uu_meta - uu_ms)):.2e}, {np.mean(np.abs(uu_meta - uu_ms)):.2e} (metres)")
print(f"v coords max, mean diff {np.max(np.abs(vv_meta - vv_ms)):.2e}, {np.mean(np.abs(vv_meta - vv_ms)):.2e} (metres)")
print(f"w coords max, mean diff {np.max(np.abs(ww_meta - ww_ms)):.2e}, {np.mean(np.abs(ww_meta - ww_ms)):.2e} (metres)")
u coords max, mean diff 5.96e-08, 7.33e-12 (metres)
v coords max, mean diff 0.00e+00, 0.00e+00 (metres)
w coords max, mean diff 4.77e-07, 2.28e-10 (metres)
OK, so this is interesting. In test_MWA.ipynb (in the same directory as this notebook), we see that the uv coords are the same, but the XYZ coords are offset. That was after reading in antenna locations from a hyperbeam measurement set. OSKAR doesn’t do precession, and I write out the current XYZ to the uvfits. Maybe hyperbeam writes out the precessed XYZ? I don’t think this is a problem, as the uv coords are the same, and all software I know of / care about reads the uvw not the XYZ. But it’s worth noting.
Let’s see if we can make flagged dipole beams¶
Again, I’ve used OSKAR to make a measurement set via create_OSKAR-MWA_ms/make_oskar-mwa_beam.ipynb. This time I’ve made the first 16 tiles have a flagged dipole each; this should give us visible beam shape deformities. First plot them on the sky to check we are getting the right thing. Then we’ll look inside the visibilities for a single off-zenith point source as different beams means the amplitudes of the visis should vary. We’ll also do two time steps, and two frequencies, to check we’re calling the right beam for the right time/frequency. Note this takes about 6 minutes on my deskop, even spread across 8 threads.
##Setup a grid of RA/Dec on the sky
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())
num_comps = len(ras)
all_freqs = np.array([160e+6, 230e+6])
all_times = np.array([observing_time, observing_time + TimeDelta(3*3600.0, format='sec')])
num_beams = 1
station_ids = np.arange(16)
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=""
ms_path="create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms"
apply_beam_norms = True
iau_order = True
parallactic_rotate = True
element_response = "skala40_wave"
num_threads = 8
all_jones = run_everybeam_over_threads(num_threads, ms_path,
coeff_path,
ras, decs,
np.radians(ra0), np.radians(dec0),
latitudes, lsts,
all_times, all_freqs, station_ids,
apply_beam_norms=apply_beam_norms,
iau_order=iau_order,
parallactic_rotate=parallactic_rotate,
element_response_model=element_response)
Thread 1 processing coords 512 to 1024Thread 0 processing coords 0 to 512Thread 3 processing coords 1536 to 2048Thread 2 processing coords 1024 to 1536
Thread 4 processing coords 2048 to 2560Thread 5 processing coords 2560 to 3072
Thread 6 processing coords 3072 to 3584Thread 7 processing coords 3584 to 4096
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:00:41 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
Thread 0 finished
Thread 1 finished
Thread 2 finished
Thread 5 finished
Thread 4 finished
Thread 7 finished
Thread 3 finished
Thread 6 finished
def plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs):
fig, axs = plt.subplots(4, 4, figsize=(12, 12), layout='constrained',
subplot_kw={'projection': wcs})
for col in range(4):
for row in range(4):
beam_ind = col*4 + row
gx = all_jones[beam_ind, time_ind, freq_ind, :, 0,0]
gx.shape = (nside, nside)
im = axs[row, col].imshow(np.abs(gx), origin='lower', vmax=1, vmin=0)
axs[row, col].grid(alpha=0.5)
lon = axs[row, col].coords[0]
lat = axs[row, col].coords[1]
lon.set_ticks_visible(False)
lat.set_ticks_visible(False)
lon.set_axislabel(' ')
lat.set_axislabel(' ')
if col > 0:
lat.set_ticklabel_visible(False)
if row < 3:
lon.set_ticklabel_visible(False)
if col == 3:
plt.colorbar(im, ax=axs[row, col])
fig.suptitle(f"Time: {all_times[time_ind].isot}, Freq: {all_freqs[freq_ind]/1e+6} MHz")
plt.show()
time_ind = 0
freq_ind = 0
plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs)
Great success, we can see our beams are deformed in different ways for the first 16 dipoles. Let’s try a different frequency:
time_ind = 0
freq_ind = 1
plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs)
Now we’ll try 3 hours later in time:
time_ind = 1
freq_ind = 1
plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs)
OK, this looks really weird with a second beam maximum not being at phase centre. But we can get the same thing using hyperbeam; it’s just a phased array MWA thing.
##Then use erfa to convert these values into azs, els
has = lsts[1] - 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, np.radians(MWA_LAT))
za_grid = np.pi/2 - els
beam = mwa_hyperbeam.FEEBeam(environ['MWA_FEE_HDF5'])
freq = 230e+6
delays = [4, 8, 12, 16]*4
jones = beam.calc_jones_array(az_grid, za_grid, freq, delays, [1]*16, True, np.radians(MWA_LAT), 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)
all_gx = jones[:,0]
all_gx = np.abs(all_gx).reshape(nside, nside)
mask = np.ones((nside, nside))
mask[els.reshape(nside, nside) < 0] = np.nan
fig, axs = plt.subplots(1, 1, figsize=(8, 6), subplot_kw={'projection': wcs})
axs.grid()
axs.imshow(all_gx*mask, origin='lower')
plt.show()
Now let’s check the visibilities. First of all, we’ll check all the visis are the same when we use the same primary beam for all tiles. Then we’ll play with the flagged tiles to see what happens
##make some more single point source sky models, but this time for an offzenith source
make_sky_models(ra0 + 15, dec0)
freq_reso = 40e+3
low_freq = 160e+6
num_freq_chans = 2
primary_beam = "everybeam_OSKAR"
uvfits_name = f"noflag_{primary_beam}"
cat_name = 'I_source.fits'
##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
command += f'--time_res=2 --num_time_steps=1 --IAU_order '
command += f'--station_id=20 '
command += f'--beam_ms_path=create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms '
command += '--num_threads=1'
call(command, 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.5.0')
Successful readonly open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms: 22 columns, 8128 rows
Successful readonly open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms::FIELD: 9 columns, 1 rows
Successful readonly open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms/ANTENNA: 8 columns, 128 rows
2025-03-06 09:06:49 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-06 09:06:49 - INFO - Input arguments after parsing:
Phase centre: 0.00455, -26.70332 deg
Array central latitude: -26.703 deg
Array central longitude: 116.671 deg
Array central height: 377.827 m
Lowest channel frequency: 160.000 MHz
Channel frequency resolution: 40.000 kHz
Coarse band bandwidth: 0.040 MHz
Num channels per coarse band: 2
Start date: 2024-07-21T20:13:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 128 antenna positions from measurement set: create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-03-06 09:06:49 - INFO - Obs epoch initial LST was 0.0087300922 deg
2025-03-06 09:06:49 - INFO - Setting initial J2000 LST to 359.6931723607 deg
2025-03-06 09:06:49 - INFO - Setting initial mjd to 60512.8423726853
2025-03-06 09:06:49 - INFO - After precession initial latitude of the array is -26.8398178265 deg
2025-03-06 09:06:49 - INFO - Precessing array layout to J2000
2025-03-06 09:06:49 - INFO - Will run with EveryBeam OSKAR primary beam, based on this measurement set:
create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms
Primary beam is pointed at RA,Dec = 0.0, -29.999999999999996 deg
2025-03-06 09:06:49 - 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-03-06 09:06:49 - INFO - Sky model mapping took 0.0 mins
2025-03-06 09:06:49 - INFO - Have read in 1 components
2025-03-06 09:06:49 - INFO - After cropping there are 1 components
2025-03-06 09:06:49 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-06 09:06:50 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-03-06 09:06:50 - INFO - Reading set 0 sky models
2025-03-06 09:06:50 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 09:06:50 - INFO - Finshed sky set 0 reading thread num 0 in 0.0 seconds
2025-03-06 09:06:50 - INFO - Sending Sky set 0 to the GPU
POINTING LOFAR beam 0.0 -30.0
2025-03-06 09:06:51 - INFO - Sky set 0 has returned from the GPU after 0.4 seconds
2025-03-06 09:06:51 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-03-06 09:06:51 - INFO - Full run took 0:00:01.996194
2025-03-06 09:06:51 - INFO - WODEN is done. Closing the log. S'later
2025-03-05 23:06:51 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:06:51 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:06:51 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
0
Read in the visibilities and plot the amplitudes of the first freq
XX, XY, YX, YY = read_uvfits('noflag_everybeam_OSKAR_band01.uvfits')
##take the XX and YY for the first freq and do some histograms
fig, axs = plt.subplots(1, 1, figsize=(12, 6))
axs.hist(np.abs(XX[:, 0]), bins=100, histtype='step', color='C0', label='XX')
axs.hist(np.abs(YY[:, 0]), bins=100, histtype='step', color='C1', label='YY')
xx_noflag_osk = np.mean(np.abs(XX[:, 0]))
yy_noflag_osk = np.mean(np.abs(YY[:, 0]))
axs.set_xlabel('Visi amplitude (Jy)')
axs.set_ylabel('Count')
axs.set_yscale('log')
axs.legend()
plt.show()
As expected, the XX and YY have different amplitudes, but all visibilities are the same. Now let’s check the flagged tiles. First we’ll estimate how many visibilities should be different, given only the first 16 tiles have different beams.
unchanged = 0
num_ants = 128
num_baseline = 0
first_unchanged = np.nan
for i in range(num_ants - 1):
for j in range(i+1, num_ants):
if i > 15 and j > 15:
if np.isnan(first_unchanged):
first_unchanged = num_baseline
unchanged += 1
num_baseline += 1
print(f"Percentage of expected unchanged baselines: {unchanged} / {num_baseline} = {unchanged/num_baseline*100:.2f}%")
print(f"First unchanged baseline: {first_unchanged}")
Percentage of expected unchanged baselines: 6216 / 8128 = 76.48%
First unchanged baseline: 1912
##parameters for the simulation; setting the date sets an LST of about 0 deg
##to match the phase centre
freq_reso = 40e+3
low_freq = 160e+6
num_freq_chans = 2
primary_beam = "everybeam_OSKAR"
uvfits_name = f"flagged_{primary_beam}"
cat_name = 'I_source.fits'
print(ra0, dec0)
##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
command += f'--time_res=2 --num_time_steps=1 --IAU_order '
command += f'--beam_ms_path=create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms '
command += '--num_threads=1'
call(command, shell=True)
0.0045520088542088475 -26.703319405555554
/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.5.0')
Successful readonly open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms: 22 columns, 8128 rows
Successful readonly open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms::FIELD: 9 columns, 1 rows
Successful readonly open of default-locked table create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms/ANTENNA: 8 columns, 128 rows
2025-03-06 09:06:52 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-06 09:06:52 - INFO - Input arguments after parsing:
Phase centre: 0.00455, -26.70332 deg
Array central latitude: -26.703 deg
Array central longitude: 116.671 deg
Array central height: 377.827 m
Lowest channel frequency: 160.000 MHz
Channel frequency resolution: 40.000 kHz
Coarse band bandwidth: 0.040 MHz
Num channels per coarse band: 2
Start date: 2024-07-21T20:13:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 128 antenna positions from measurement set: create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-03-06 09:06:53 - INFO - Obs epoch initial LST was 0.0087300922 deg
2025-03-06 09:06:53 - INFO - Setting initial J2000 LST to 359.6931723607 deg
2025-03-06 09:06:53 - INFO - Setting initial mjd to 60512.8423726853
2025-03-06 09:06:53 - INFO - After precession initial latitude of the array is -26.8398178265 deg
2025-03-06 09:06:53 - INFO - Precessing array layout to J2000
2025-03-06 09:06:53 - INFO - Will run with EveryBeam OSKAR primary beam, based on this measurement set:
create_OSKAR-MWA_ms/OSKAR-MWA-flags-layout.ms
Primary beam is pointed at RA,Dec = 0.0, -29.999999999999996 deg
2025-03-06 09:06:53 - 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-03-06 09:06:53 - INFO - Sky model mapping took 0.0 mins
2025-03-06 09:06:53 - INFO - Have read in 1 components
2025-03-06 09:06:53 - INFO - After cropping there are 1 components
2025-03-06 09:06:53 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-06 09:06:54 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-03-06 09:06:54 - INFO - Reading set 0 sky models
2025-03-06 09:06:54 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 09:06:54 - INFO - Finshed sky set 0 reading thread num 0 in 0.0 seconds
2025-03-06 09:06:54 - INFO - Sending Sky set 0 to the GPU
POINTING LOFAR beam 0.0 -30.0
2025-03-05 23:06:55 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292) Leap second table TAI_UTC seems out-of-date.
2025-03-05 23:06:55 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ Until the table is updated (see the CASA documentation or your system admin),
2025-03-05 23:06:55 SEVERE MeasTable::dUTC(Double) (file ./measures/Measures/MeasTable.cc, line 4292)+ times and coordinates derived from UTC could be wrong by 1s or more.
2025-03-06 09:06:55 - INFO - Sky set 0 has returned from the GPU after 0.9 seconds
2025-03-06 09:06:55 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-03-06 09:06:55 - INFO - Full run took 0:00:02.582723
2025-03-06 09:06:55 - INFO - WODEN is done. Closing the log. S'later
0
Note I’ve deleted the outputs of the cell above, because EveryBeam spits out a bunch of error messages with the OSKAR beam warning about frequencies mismatching, and not having a central direction set. The outputs seems fine. It’s a TODO to silence these warnings.
XX, XY, YX, YY = read_uvfits('flagged_everybeam_OSKAR_band01.uvfits')
##take the XX and YY for the first freq and do some histograms
fig, axs = plt.subplots(1, 1, figsize=(12, 6))
axs.hist(np.abs(XX[:, 0]), bins=60, histtype='step', color='C0', label='XX')
axs.hist(np.abs(YY[:, 0]), bins=60, histtype='step', color='C1', label='YY')
axs.axhline(unchanged, color='k', linestyle='-', alpha=0.5, label='Number expected to be unchanged')
axs.axvline(xx_noflag_osk, color='C0', linestyle='--', label='XX no flag value')
axs.axvline(yy_noflag_osk, color='C1', linestyle='--', label='YY no flag value')
axs.set_xlabel('Visi amplitude (Jy)')
axs.set_ylabel('Count')
axs.set_yscale('log')
axs.legend()
plt.show()
OK, this makes sense. I only put dipole flags on the first 16 tiles, so a large number of the visibilities should have the same value (which are the big central peaks of each distribution). The XX and YY are different amplitudes, as I put the off-zenith. This means the XX and YY beams should have different sensitivities. They both have the same dipole flags however, so the distributions should have a similar shape as they do.
Maybe counter-intuitive, but with the change in beam shape per tiles, some beams can become relatively more sensitive than others. As I’m normalising each beam individually towards the phase centre, this means we can high a higher ampltiude than the full 16 dipole beam. This behaviour probably needs to be changed to have some kind of overall normalisation, but this is a good first step. Doesn’t seem to be a problem in this example however.
Let’s compare that to a hyperbeam model with the same flags
##parameters for the simulation; setting the date sets an LST of about 0 deg
##to match the phase centre
freq_reso = 40e+3
low_freq = 160e+6
num_freq_chans = 2
primary_beam = "MWA_FEE"
uvfits_name = f"noflags_{primary_beam}"
cat_name = 'I_source.fits'
##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
command += f'--time_res=2 --num_time_steps=1 --IAU_order '
command += f'--metafits=create_OSKAR-MWA_ms/1088285600_flags.metafits '
command += f'--num_threads=1'
call(command, 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.5.0')
2025-03-06 09:06:57 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-06 09:06:57 - INFO - Input arguments after parsing:
Phase centre: 0.00455, -26.70332 deg
Array central latitude: -26.703 deg
Array central longitude: 116.671 deg
Array central height: 377.827 m
Lowest channel frequency: 160.000 MHz
Channel frequency resolution: 40.000 kHz
Coarse band bandwidth: 1.280 MHz
Num channels per coarse band: 2
Start date: 2024-07-21T20:13:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 128 antenna positions from metafits file: create_OSKAR-MWA_ms/1088285600_flags.metafits
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-03-06 09:06:58 - INFO - Obs epoch initial LST was 0.0087300922 deg
2025-03-06 09:06:58 - INFO - Setting initial J2000 LST to 359.6931723607 deg
2025-03-06 09:06:58 - INFO - Setting initial mjd to 60512.8423726853
2025-03-06 09:06:58 - INFO - After precession initial latitude of the array is -26.8398178265 deg
2025-03-06 09:06:58 - INFO - Precessing array layout to J2000
2025-03-06 09:06:58 - 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-03-06 09:06:58 - 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-03-06 09:06:58 - INFO - Sky model mapping took 0.0 mins
2025-03-06 09:06:58 - INFO - Have read in 1 components
2025-03-06 09:06:58 - INFO - After cropping there are 1 components
2025-03-06 09:06:58 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-06 09:06:58 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-03-06 09:06:58 - INFO - Reading set 0 sky models
2025-03-06 09:06:58 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 09:06:58 - INFO - Sending Sky set 0 to the GPU
2025-03-06 09:06:58 - INFO - Sky set 0 has returned from the GPU after 0.1 seconds
2025-03-06 09:06:58 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-03-06 09:06:58 - INFO - Full run took 0:00:01.607278
2025-03-06 09:06:58 - INFO - WODEN is done. Closing the log. S'later
0
XX, XY, YX, YY = read_uvfits('noflags_MWA_FEE_band01.uvfits')
##take the XX and YY for the first freq and do some histograms
fig, axs = plt.subplots(1, 1, figsize=(12, 6))
axs.hist(np.abs(XX[:, 0]), bins=100, histtype='step', color='C0', label='XX')
axs.hist(np.abs(YY[:, 0]), bins=100, histtype='step', color='C1', label='YY')
xx_noflag = np.mean(np.abs(XX[:, 0]))
yy_noflag = np.mean(np.abs(YY[:, 0]))
axs.set_xlabel('Visi amplitude (Jy)')
axs.set_ylabel('Count')
axs.set_yscale('log')
axs.legend()
plt.show()
##parameters for the simulation; setting the date sets an LST of about 0 deg
##to match the phase centre
freq_reso = 70e+6
low_freq = 160e+6
num_freq_chans = 2
primary_beam = "MWA_FEE"
uvfits_name = f"flagged_{primary_beam}"
cat_name = 'I_source.fits'
##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam={primary_beam} '
command += f'--lowest_channel_freq={low_freq} --freq_res={freq_reso} '
command += f'--num_freq_channels={num_freq_chans} --band_nums=1 '
command += f'--time_res=2 --num_time_steps=1 --IAU_order '
command += f'--metafits_filename=create_OSKAR-MWA_ms/1088285600_flags.metafits '
command += f'--use_MWA_dipflags '
command += f'--num_threads=1'
call(command, 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.5.0')
Num tiles with dipole flags: 16
2025-03-06 09:07:00 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-06 09:07:00 - INFO - Input arguments after parsing:
Phase centre: 0.00455, -26.70332 deg
Array central latitude: -26.703 deg
Array central longitude: 116.671 deg
Array central height: 377.827 m
Lowest channel frequency: 160.000 MHz
Channel frequency resolution: 70000.000 kHz
Coarse band bandwidth: 1.280 MHz
Num channels per coarse band: 2
Start date: 2024-07-21T20:13:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 128 antenna positions from metafits file: create_OSKAR-MWA_ms/1088285600_flags.metafits
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-03-06 09:07:01 - INFO - Obs epoch initial LST was 0.0087300922 deg
2025-03-06 09:07:01 - INFO - Setting initial J2000 LST to 359.6931723607 deg
2025-03-06 09:07:01 - INFO - Setting initial mjd to 60512.8423726853
2025-03-06 09:07:01 - INFO - After precession initial latitude of the array is -26.8398178265 deg
2025-03-06 09:07:01 - INFO - Precessing array layout to J2000
2025-03-06 09:07:01 - 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]
will use dipole amplitudes from given metafits file
will use dipole flags from given metafits file
metafits file: create_OSKAR-MWA_ms/1088285600_flags.metafits
2025-03-06 09:07:01 - 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-03-06 09:07:01 - INFO - Sky model mapping took 0.0 mins
2025-03-06 09:07:01 - INFO - Have read in 1 components
2025-03-06 09:07:01 - INFO - After cropping there are 1 components
2025-03-06 09:07:01 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-06 09:07:01 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-03-06 09:07:01 - INFO - Reading set 0 sky models
2025-03-06 09:07:01 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-06 09:07:02 - INFO - Sending Sky set 0 to the GPU
2025-03-06 09:07:02 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-03-06 09:07:02 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-03-06 09:07:02 - INFO - Full run took 0:00:01.880184
2025-03-06 09:07:02 - INFO - WODEN is done. Closing the log. S'later
0
XX, XY, YX, YY = read_uvfits('flagged_MWA_FEE_band01.uvfits')
##take the XX and YY for the first freq and do some histograms
fig, axs = plt.subplots(1, 1, figsize=(12, 6))
axs.hist(np.abs(XX[:, 0]), bins=50, histtype='step', color='C0', label='XX')
axs.hist(np.abs(YY[:, 0]), bins=50, histtype='step', color='C1', label='YY')
axs.axhline(unchanged, color='k', linestyle='-', alpha=0.5, label='Number expected to be unchanged')
axs.axvline(xx_noflag, color='C0', linestyle='--', label='XX no flag value')
axs.axvline(yy_noflag, color='C1', linestyle='--', label='YY no flag value')
axs.set_xlabel('Visi amplitude (Jy)')
axs.set_ylabel('Count')
axs.set_yscale('log')
axs.legend()
plt.show()
So a similar result to the everybeam experiment, with two important differences. Firstly, the hyperbeam model is normalised to the full 16 dipole beam, so any dipole flag causes a drop in sensitivity. Hence all visibilities are either unchanged, or lower than full 16 dipole beam. Secondly, the XX and YY amps have a slightly different shape. Maybe this has something to do with the different models having different mutual coupling effects, as well as the SKA beam is a christmas tree whereas the MWA dipole is a bowtie.