EveryBeam LOFAR LBA integration tests¶
These tests are all setup and explained in test_installation/everybeam/test_NWA.ipynb, so check that out first if you haven’t already. In short, we’ll check we can make images of the LBA primary beam, that the beams interact with Stokes polarisation correctly, and that we can change frequency and time correctly.
Goes without saying, but these tests rely on everybeam being installed.
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
import erfa
##A bunch of test and plotting code lives here so we can use it in multile notebooks
from wodenpy.primary_beam.use_everybeam import run_everybeam, run_everybeam_over_threads
##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
from wodenpy.array_layout.precession import RTS_Precess_LST_Lat_to_J2000
C = c.to('m/s').value
ra0 = 0
dec0 = 53
lofar_lat = 52.905329712
lofar_long = 6.867996528
##pick a time/date that sticks our phase centre overhead
date = "2024-07-21T03:35:00"
##Assume that the OSKAR telescope is near the MWA??
lofar_location = EarthLocation(lat=lofar_lat*u.deg,
lon=lofar_long*u.deg)
observing_time = Time(date, scale='utc', location=lofar_location)
##Grab the LST
LST = observing_time.sidereal_time('mean')
LST_deg = LST.value*15
print(f"LST: {LST_deg} deg, RA: {ra0}")
LST: 0.01862197128531748 deg, RA: 0
OK, so the most efficient way I’ve found of calling EveryBeam for LOFAR does not allow you to directly control the beam pointing. Instead, you have to edit the underlying Measurement Set (MS) to change the pointing. So we’ll do that now.
from wodenpy.primary_beam.use_everybeam import create_filtered_ms
create_filtered_ms("lba.MS", "pointed_LBA.ms", np.radians(ra0), np.radians(dec0))
ms_path = "pointed_LBA.ms"
Successful readonly open of default-locked table lba.MS: 24 columns, 4218 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Image the beam on the sky¶
First up, let’s see if we can just plot a LOFAR station beam to check it looks sensible.
##setup coords
nside = 100
radec_reso = 60 / 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())
freq = 50e+6
station_id = 0
##Stick all our variables in lists as `run_everybeam` is designed to be
## run over multiple times/frequencies
times = [observing_time]
freqs = [freq]
station_ids = [0]
##Everybeam doesn't need a coeff path for this beam
coeff_path=""
##We'll run EveryBeam in parallel to speed things up
num_threads=8
all_jones = run_everybeam_over_threads(num_threads,
ras, decs,
ms_path=ms_path,
times=times, freqs=freqs,
station_ids=station_ids,
beam_ra0=np.radians(ra0),beam_dec0=np.radians(dec0),
apply_beam_norms=True,
iau_order=True,
parallactic_rotate=True)
beam_ind, time_ind, freq_ind = 0, 0, 0
all_gx = all_jones[beam_ind, time_ind, freq_ind, :, 0, 0]
all_Dx = all_jones[beam_ind, time_ind, freq_ind, :, 0, 1]
all_Dy = all_jones[beam_ind, time_ind, freq_ind, :, 1, 0]
all_gy = all_jones[beam_ind, time_ind, freq_ind, :, 1, 1]
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)
This looks decent, gains look way stronger than the leakages because I’ve set up the run_everybeam function to do a normalistion towards the beam phase centre. However, these beams aren’t aligned north-east and south-west as the dipoles are different in a LOFAR station.
Off cardinal dipoles¶
This brings up a need for new functionality in WODEN to handle 45/135 deg dipoles. As noted in Table 4.1 of TMS, the combinations of Stokes IQUV parameters measured by each instrumental polarisation are different as compared to 0/90 deg dipoles. This means we need a different mixing matrix. If we label the 45/135 deg dipoles as p and q, then the instrumental visibilities are:
\[\begin{split}
\begin{bmatrix}
\mathrm{V}_{pp} \\
\mathrm{V}_{pq} \\
\mathrm{V}_{qp} \\
\mathrm{V}_{qq} \\
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & 1 & 0 \\
0 & -1 & 0 & j \\
0 & -1 & 0 & -j \\
1 & 0 & -1 & 0 \\
\end{bmatrix}
\begin{bmatrix}
\mathrm{V}_I \\
\mathrm{V}_Q \\
\mathrm{V}_U \\
\mathrm{V}_V \\
\end{bmatrix}
\end{split}\]
If we then say the beam Jones matrix is:
\[\begin{split}
\mathbf{J} =
\begin{bmatrix}
g_p & D_p \\
D_q & g_q \\
\end{bmatrix}
\end{split}\]
then the measured visibilities between antennas 1 and 2 are:
\[\begin{split}
\begin{bmatrix}
\mathrm{V}_{12\,pp} \\
\mathrm{V}_{12\,pq} \\
\mathrm{V}_{12\,qp} \\
\mathrm{V}_{12\,qq} \\
\end{bmatrix} =
\mathbf{J}_1 \otimes \mathbf{J}_2^*
\begin{bmatrix}
1 & 0 & 1 & 0 \\
0 & -1 & 0 & i \\
0 & -1 & 0 & -i \\
1 & 0 & -1 & 0 \\
\end{bmatrix}
\begin{bmatrix}
\mathrm{V}_{12\,I} \\
\mathrm{V}_{12\,Q} \\
\mathrm{V}_{12\,U} \\
\mathrm{V}_{12\,V}
\end{bmatrix}
\end{split}\]
\[\begin{split}=
\begin{bmatrix}
g_{1p} g_{2p}^* & g_{1p} D_{2p}^* & D_{1p} g_{2p}^* & D_{1p} D_{2p}^* \\
g_{1p} D_{2q}^* & g_{1p} g_{2q}^* & D_{1p} D_{2q}^* & D_{1p} g_{2q}^* \\
D_{1q} g_{2p}^* & D_{1q} D_{2p}^* & g_{1q} g_{2p}^* & g_{1q} D_{2p}^* \\
D_{1q} D_{2q}^* & D_{1q} g_{2q}^* & g_{1q} D_{2q}^* & g_{1q} g_{2q}^*
\end{bmatrix}
\begin{bmatrix}
1 & 0 & 1 & 0 \\
0 & -1 & 0 & i \\
0 & -1 & 0 & -i \\
1 & 0 & -1 & 0 \\
\end{bmatrix}
\begin{bmatrix}
\mathrm{V}_{12\,I} \\
\mathrm{V}_{12\,Q} \\
\mathrm{V}_{12\,U} \\
\mathrm{V}_{12\,V}
\end{bmatrix}
\end{split}\]
\[\begin{split}=
\begin{bmatrix}
g_{1p} g_{2p}^* + D_{1p} D_{2p}^* & -g_{1p} D_{2p}^* - D_{1p} g_{2p}^* & g_{1p} g_{2p}^* - D_{1p} D_{2p}^* & i(g_{1p} D_{2p}^* -D_{1p} g_{2p}^*) \\
g_{1p} D_{2q}^* + D_{1p} g_{2q}^* & -g_{1p} g_{2q}^* - D_{1p} D_{2q}^* & g_{1p} D_{2q}^* - D_{1p} g_{2q}^* & i(g_{1p} g_{2q}^* -D_{1p} D_{2q}^*) \\
D_{1q} g_{2p}^* + g_{1q} D_{2p}^* & -D_{1q} D_{2p}^* - g_{1q} g_{2p}^* & D_{1q} g_{2p}^* - g_{1q} D_{2p}^* & i(D_{1q} D_{2p}^* -g_{1q} g_{2p}^*) \\
D_{1q} D_{2q}^* + g_{1q} g_{2q}^* & -D_{1q} g_{2q}^* - g_{1q} D_{2q}^* & D_{1q} D_{2q}^* - g_{1q} g_{2q}^* & i(D_{1q} g_{2q}^* -g_{1q} D_{2q}^*) \\
\end{bmatrix}
\begin{bmatrix}
\mathrm{V}_{12\,I} \\
\mathrm{V}_{12\,Q} \\
\mathrm{V}_{12\,U} \\
\mathrm{V}_{12\,V}
\end{bmatrix}
\end{split}\]
\[\begin{split}
\begin{align*}
\mathrm{V}_{12\,pp} = (g_{1p} g_{2p}^* + D_{1p} D_{2p}^*)\mathrm{V}_{12\,I} - (g_{1p} D_{2p}^* + D_{1p} g_{2p}^*)\mathrm{V}_{12\,Q} \\+ (g_{1p} g_{2p}^* - D_{1p} D_{2p}^*)\mathrm{V}_{12\,U} + i(g_{1p} D_{2p}^* - D_{1p} g_{2p}^*)\mathrm{V}_{12\,V}
\end{align*}
\end{split}\]
\[\begin{split}
\begin{align*}
\mathrm{V}_{12\,pq} = (g_{1p} D_{2q}^* + D_{1p} g_{2q}^*)\mathrm{V}_{12\,I} - (g_{1p} g_{2q}^* + D_{1p} D_{2q}^*)\mathrm{V}_{12\,Q} \\+ (g_{1p} D_{2q}^* - D_{1p} g_{2q}^*)\mathrm{V}_{12\,U} + i(g_{1p} g_{2q}^* -D_{1p} D_{2q}^*)\mathrm{V}_{12\,V}
\end{align*}
\end{split}\]
\[\begin{split}
\begin{align*}
\mathrm{V}_{12\,qp} = (D_{1q} g_{2p}^* + g_{1q} D_{2p}^*)\mathrm{V}_{12\,I} - (D_{1q} D_{2p}^* + g_{1q} g_{2p}^*)\mathrm{V}_{12\,Q} \\+ (D_{1q} g_{2p}^* - g_{1q} D_{2p}^*)\mathrm{V}_{12\,U} + i(D_{1q} D_{2p}^* -g_{1q} g_{2p}^*)\mathrm{V}_{12\,V}
\end{align*}
\end{split}\]
\[\begin{split}
\begin{align*}
\mathrm{V}_{12\,qq} = (D_{1q} D_{2q}^* + g_{1q} g_{2q}^*)\mathrm{V}_{12\,I} - (D_{1q} g_{2q}^* + g_{1q} D_{2q}^*)\mathrm{V}_{12\,Q} \\+ (D_{1q} D_{2q}^* - g_{1q} g_{2q}^*)\mathrm{V}_{12\,U} + i(D_{1q} g_{2q}^* -g_{1q} D_{2q}^*)\mathrm{V}_{12\,V}
\end{align*}
\end{split}\]
Repeating the whole process but ignoring the beam, we see that
\[
\begin{align*}
\mathrm{V}_{pp} = \mathrm{V}_{I} + \mathrm{V}_{U}
\end{align*}
\]
\[
\begin{align*}
\mathrm{V}_{pq} = -\mathrm{V}_{Q} + i\mathrm{V}_{V}
\end{align*}
\]
\[
\begin{align*}
\mathrm{V}_{qp} = -\mathrm{V}_{Q} - i\mathrm{V}_{V}
\end{align*}
\]
\[
\begin{align*}
\mathrm{V}_{qq} = \mathrm{V}_{I} - \mathrm{V}_{U}
\end{align*}
\]
Rearranging this we see to recover Stokes IQUV visibilities, we obviously need a different set of equations compared to 0/90 deg dipoles. These are
\[
\begin{align*}
\mathrm{V}_{I} = \frac{\mathrm{V}_{pp} + \mathrm{V}_{qq}}{2}
\end{align*}
\]
\[
\begin{align*}
\mathrm{V}_{Q} = -\frac{\mathrm{V}_{pq} + \mathrm{V}_{qp}}{2}
\end{align*}
\]
\[
\begin{align*}
\mathrm{V}_{U} = \frac{\mathrm{V}_{pp} - \mathrm{V}_{qq}}{2}
\end{align*}
\]
\[
\begin{align*}
\mathrm{V}_{V} = \frac{\mathrm{V}_{qp} - \mathrm{V}_{pq}}{2i}
\end{align*}
\]
This maths has been added to WODEN. Loading a LOFAR beam model should trigger the correct mixing matrix to be used, but there is also a --off_cardinal_dipoles flag that can be set to force the use of the off cardinal mixing matrix.
Stokes recovery with off cardinal dipoles¶
Let’s try to recover single point source of either Stokes I, Q, U or V. First of all, we’ll do it without any primary beam effects, to check the mixing matrix. We’ll then add the primary beam effects and see how it affects the recovery.
make_sky_models(ra0, dec0)
freq_reso = 1e+6
low_freq = 100e+6
high_freq = 150e+6
num_freq_chans = int((high_freq - low_freq) / freq_reso)
primary_beam = "none"
for pol in ['I', 'Q', 'U', 'V']:
# for pol in ['I']:
uvfits_name = f"stokes{pol}_{primary_beam}"
cat_name = f'{pol}_source.fits'
cmd = "run_woden.py "
##The command to run WODEN
cmd += f'--ra0={ra0} '
cmd += f'--dec0={dec0} '
cmd += f'--latitude={lofar_lat} '
cmd += f'--longitude={lofar_long} '
cmd += f'--date={date} '
cmd += f'--output_uvfits_prepend={uvfits_name} '
cmd += f'--cat_filename={cat_name} '
cmd += f'--primary_beam={primary_beam} '
cmd += f'--lowest_channel_freq={low_freq} '
cmd += f'--freq_res={freq_reso} '
cmd += f'--num_freq_channels={num_freq_chans} '
cmd += f'--band_nums=1 '
cmd += f'--time_res=2 '
cmd += f'--num_time_steps=1 '
cmd += f'--IAU_order '
cmd += f'--off_cardinal_dipoles '
cmd += '--num_threads=1 '
cmd += f'--beam_ms_path=pointed_LBA.ms '
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')
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful readonly open of default-locked table pointed_LBA.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
2025-10-22 14:33:20 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: c8b7005
2025-10-22 14:33:20 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 53.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 377.827 m
Lowest channel frequency: 100.000 MHz
Channel frequency resolution: 1000.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 50
Start date: 2024-07-21T03:35:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 37 antenna positions from measurement set: pointed_LBA.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-22 14:33:20 - INFO - Obs epoch initial LST was 0.0228000546 deg
2025-10-22 14:33:20 - INFO - Setting initial J2000 LST to 359.7122369008 deg
2025-10-22 14:33:20 - INFO - Setting initial mjd to 60512.1493171296
2025-10-22 14:33:20 - INFO - After precession initial latitude of the array is 52.7688502802 deg
2025-10-22 14:33:20 - INFO - Precessing array layout to J2000
2025-10-22 14:33:20 - INFO - No primary beam selected, no beam attenuation will be applied.
2025-10-22 14:33:20 - 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-22 14:33:20 - INFO - Sky model mapping took 0.0 mins
2025-10-22 14:33:20 - INFO - Have read in 1 components
2025-10-22 14:33:20 - INFO - After cropping there are 1 components
2025-10-22 14:33:20 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-22 14:33:21 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-10-22 14:33:21 - INFO - Reading set 0 sky models
2025-10-22 14:33:21 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:33:21 - INFO - Sending Sky set 0 to the GPU
2025-10-22 14:33:21 - INFO - Sky set 0 has returned from the GPU after 0.1 seconds
2025-10-22 14:33:21 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-22 14:33:21 - INFO - Full run took 0:00:01.827117
2025-10-22 14:33:21 - 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')
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful readonly open of default-locked table pointed_LBA.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
2025-10-22 14:33:23 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: c8b7005
2025-10-22 14:33:23 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 53.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 377.827 m
Lowest channel frequency: 100.000 MHz
Channel frequency resolution: 1000.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 50
Start date: 2024-07-21T03:35:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 37 antenna positions from measurement set: pointed_LBA.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-22 14:33:24 - INFO - Obs epoch initial LST was 0.0228000546 deg
2025-10-22 14:33:24 - INFO - Setting initial J2000 LST to 359.7122369008 deg
2025-10-22 14:33:24 - INFO - Setting initial mjd to 60512.1493171296
2025-10-22 14:33:24 - INFO - After precession initial latitude of the array is 52.7688502802 deg
2025-10-22 14:33:24 - INFO - Precessing array layout to J2000
2025-10-22 14:33:24 - INFO - No primary beam selected, no beam attenuation will be applied.
2025-10-22 14:33:24 - INFO - Doing the initial mapping of sky model
2025-10-22 14:33:24 - INFO - Sky model mapping took 0.0 mins
2025-10-22 14:33:24 - INFO - Have read in 1 components
2025-10-22 14:33:24 - INFO - After cropping there are 1 components
2025-10-22 14:33:24 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-22 14:33:24 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-10-22 14:33:24 - INFO - Reading set 0 sky models
2025-10-22 14:33:24 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:33:24 - INFO - Sending Sky set 0 to the GPU
2025-10-22 14:33:24 - INFO - Sky set 0 has returned from the GPU after 0.1 seconds
2025-10-22 14:33:25 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-22 14:33:25 - INFO - Full run took 0:00:01.946008
2025-10-22 14:33:25 - 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')
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful readonly open of default-locked table pointed_LBA.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
2025-10-22 14:33:26 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: c8b7005
2025-10-22 14:33:26 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 53.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 377.827 m
Lowest channel frequency: 100.000 MHz
Channel frequency resolution: 1000.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 50
Start date: 2024-07-21T03:35:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 37 antenna positions from measurement set: pointed_LBA.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-22 14:33:27 - INFO - Obs epoch initial LST was 0.0228000546 deg
2025-10-22 14:33:27 - INFO - Setting initial J2000 LST to 359.7122369008 deg
2025-10-22 14:33:27 - INFO - Setting initial mjd to 60512.1493171296
2025-10-22 14:33:27 - INFO - After precession initial latitude of the array is 52.7688502802 deg
2025-10-22 14:33:27 - INFO - Precessing array layout to J2000
2025-10-22 14:33:27 - INFO - No primary beam selected, no beam attenuation will be applied.
2025-10-22 14:33:27 - INFO - Doing the initial mapping of sky model
2025-10-22 14:33:27 - INFO - Sky model mapping took 0.0 mins
2025-10-22 14:33:27 - INFO - Have read in 1 components
2025-10-22 14:33:27 - INFO - After cropping there are 1 components
2025-10-22 14:33:27 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-22 14:33:27 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-10-22 14:33:27 - INFO - Reading set 0 sky models
2025-10-22 14:33:27 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:33:28 - INFO - Sending Sky set 0 to the GPU
2025-10-22 14:33:28 - INFO - Sky set 0 has returned from the GPU after 0.1 seconds
2025-10-22 14:33:28 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-22 14:33:28 - INFO - Full run took 0:00:01.878011
2025-10-22 14:33:28 - 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')
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful readonly open of default-locked table pointed_LBA.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
2025-10-22 14:33:29 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: c8b7005
2025-10-22 14:33:29 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 53.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 377.827 m
Lowest channel frequency: 100.000 MHz
Channel frequency resolution: 1000.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 50
Start date: 2024-07-21T03:35:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 37 antenna positions from measurement set: pointed_LBA.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-22 14:33:30 - INFO - Obs epoch initial LST was 0.0228000546 deg
2025-10-22 14:33:30 - INFO - Setting initial J2000 LST to 359.7122369008 deg
2025-10-22 14:33:30 - INFO - Setting initial mjd to 60512.1493171296
2025-10-22 14:33:30 - INFO - After precession initial latitude of the array is 52.7688502802 deg
2025-10-22 14:33:30 - INFO - Precessing array layout to J2000
2025-10-22 14:33:30 - INFO - No primary beam selected, no beam attenuation will be applied.
2025-10-22 14:33:30 - 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-22 14:33:30 - INFO - Sky model mapping took 0.0 mins
2025-10-22 14:33:30 - INFO - Have read in 1 components
2025-10-22 14:33:30 - INFO - After cropping there are 1 components
2025-10-22 14:33:30 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-22 14:33:31 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-10-22 14:33:31 - INFO - Reading set 0 sky models
2025-10-22 14:33:31 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:33:31 - INFO - Sending Sky set 0 to the GPU
2025-10-22 14:33:31 - INFO - Sky set 0 has returned from the GPU after 0.1 seconds
2025-10-22 14:33:31 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-22 14:33:31 - INFO - Full run took 0:00:01.946729
2025-10-22 14:33:31 - INFO - WODEN is done. Closing the log. S'later
for pol in ['I', 'Q', 'U', 'V']:
uvfits_name = f"stokes{pol}_none"
XX, XY, YX, YY = read_uvfits(f'{uvfits_name}_band01.uvfits')
##on_cardinal=False means we're using the maths above to recombine the
##XX, XY, YX, YY into Stokes I, Q, U, V
test_stokes_recovery(pol, 'none', atol=5e-2, on_cardinal=False)
# ##uncomment to print out the exact values recovered
# ##pick a random baseline to print, they should all be the sam
# baseline = 0
# ##Supply on_cardinal=False to use the right maths
# recover_I, recover_Q, recover_U, recover_V = convert_inst_to_stokes(XX[baseline], XY[baseline], YX[baseline], YY[baseline],
# on_cardinal=False)
# print(f"Input {pol}, Recover I {recover_I[0].real:.2f}, Q {recover_Q[0].real:.2f}, U {recover_U[0].real:.2f}, V {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
Sweeeeet, works in the case of no primary beam. How about when we add the LOFAR primary beam?
make_sky_models(ra0, dec0)
freq_reso=1e+6
low_freq=100e+6
high_freq = 150e+6
num_freq_chans = int((high_freq - low_freq) / freq_reso)
primary_beam = "everybeam_LOFAR"
for pol in ['I', 'Q', 'U', 'V']:
uvfits_name = f"stokes{pol}_{primary_beam}"
cat_name = f'{pol}_source.fits'
cmd = "run_woden.py "
##The command to run WODEN
cmd += f'--ra0={ra0} '
cmd += f'--dec0={dec0} '
cmd += f'--latitude={lofar_lat} '
cmd += f'--longitude={lofar_long} '
cmd += f'--date={date} '
cmd += f'--output_uvfits_prepend={uvfits_name} '
cmd += f'--cat_filename={cat_name} '
cmd += f'--primary_beam={primary_beam} '
cmd += f'--lowest_channel_freq={low_freq} '
cmd += f'--freq_res={freq_reso} '
cmd += f'--num_freq_channels={num_freq_chans} '
cmd += f'--band_nums=1 '
cmd += f'--time_res=2 '
cmd += f'--num_time_steps=1 '
cmd += f'--IAU_order '
cmd += f'--off_cardinal_dipoles '
cmd += '--num_threads=1 '
cmd += f'--beam_ms_path=pointed_LBA.ms '
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')
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful readonly open of default-locked table pointed_LBA.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
2025-10-22 14:33:50 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: c8b7005
2025-10-22 14:33:50 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 53.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 0.000 m
Lowest channel frequency: 100.000 MHz
Channel frequency resolution: 1000.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 50
Start date: 2024-07-21T03:35:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 37 antenna positions from measurement set: pointed_LBA.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-22 14:33:51 - INFO - Obs epoch initial LST was 0.0228000546 deg
2025-10-22 14:33:51 - INFO - Setting initial J2000 LST to 359.7122369008 deg
2025-10-22 14:33:51 - INFO - Setting initial mjd to 60512.1493171296
2025-10-22 14:33:51 - INFO - After precession initial latitude of the array is 52.7688502802 deg
2025-10-22 14:33:51 - INFO - Precessing array layout to J2000
2025-10-22 14:33:51 - INFO - Will run with EveryBeam LOFAR primary beam, based on this measurement set:
pointed_LBA.ms
Primary beam is pointed at RA,Dec = 0.0, 53.0 deg
2025-10-22 14:33:51 - 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-22 14:33:51 - INFO - Sky model mapping took 0.0 mins
2025-10-22 14:33:51 - INFO - Have read in 1 components
2025-10-22 14:33:51 - INFO - After cropping there are 1 components
2025-10-22 14:33:51 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-22 14:33:52 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-10-22 14:33:52 - INFO - Reading set 0 sky models
2025-10-22 14:33:52 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:33:52 - INFO - Sending Sky set 0 to the GPU
2025-10-22 14:33:52 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-22 14:33:52 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-22 14:33:52 - INFO - Full run took 0:00:02.015367
2025-10-22 14:33:52 - 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')
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful readonly open of default-locked table pointed_LBA.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
2025-10-22 14:33:54 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: c8b7005
2025-10-22 14:33:54 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 53.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 0.000 m
Lowest channel frequency: 100.000 MHz
Channel frequency resolution: 1000.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 50
Start date: 2024-07-21T03:35:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 37 antenna positions from measurement set: pointed_LBA.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-22 14:33:55 - INFO - Obs epoch initial LST was 0.0228000546 deg
2025-10-22 14:33:55 - INFO - Setting initial J2000 LST to 359.7122369008 deg
2025-10-22 14:33:55 - INFO - Setting initial mjd to 60512.1493171296
2025-10-22 14:33:55 - INFO - After precession initial latitude of the array is 52.7688502802 deg
2025-10-22 14:33:55 - INFO - Precessing array layout to J2000
2025-10-22 14:33:55 - INFO - Will run with EveryBeam LOFAR primary beam, based on this measurement set:
pointed_LBA.ms
Primary beam is pointed at RA,Dec = 0.0, 53.0 deg
2025-10-22 14:33:55 - INFO - Doing the initial mapping of sky model
2025-10-22 14:33:55 - INFO - Sky model mapping took 0.0 mins
2025-10-22 14:33:55 - INFO - Have read in 1 components
2025-10-22 14:33:55 - INFO - After cropping there are 1 components
2025-10-22 14:33:55 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-22 14:33:55 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-10-22 14:33:55 - INFO - Reading set 0 sky models
2025-10-22 14:33:55 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:33:55 - INFO - Sending Sky set 0 to the GPU
2025-10-22 14:33:56 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-22 14:33:56 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-22 14:33:56 - INFO - Full run took 0:00:02.008347
2025-10-22 14:33:56 - 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')
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful readonly open of default-locked table pointed_LBA.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
2025-10-22 14:33:57 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: c8b7005
2025-10-22 14:33:57 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 53.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 0.000 m
Lowest channel frequency: 100.000 MHz
Channel frequency resolution: 1000.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 50
Start date: 2024-07-21T03:35:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 37 antenna positions from measurement set: pointed_LBA.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-22 14:33:58 - INFO - Obs epoch initial LST was 0.0228000546 deg
2025-10-22 14:33:58 - INFO - Setting initial J2000 LST to 359.7122369008 deg
2025-10-22 14:33:58 - INFO - Setting initial mjd to 60512.1493171296
2025-10-22 14:33:58 - INFO - After precession initial latitude of the array is 52.7688502802 deg
2025-10-22 14:33:58 - INFO - Precessing array layout to J2000
2025-10-22 14:33:58 - INFO - Will run with EveryBeam LOFAR primary beam, based on this measurement set:
pointed_LBA.ms
Primary beam is pointed at RA,Dec = 0.0, 53.0 deg
2025-10-22 14:33:58 - INFO - Doing the initial mapping of sky model
2025-10-22 14:33:58 - INFO - Sky model mapping took 0.0 mins
2025-10-22 14:33:58 - INFO - Have read in 1 components
2025-10-22 14:33:58 - INFO - After cropping there are 1 components
2025-10-22 14:33:58 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-22 14:33:58 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-10-22 14:33:58 - INFO - Reading set 0 sky models
2025-10-22 14:33:58 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:33:59 - INFO - Sending Sky set 0 to the GPU
2025-10-22 14:33:59 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-22 14:33:59 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-22 14:33:59 - INFO - Full run took 0:00:02.003016
2025-10-22 14:33:59 - 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')
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful readonly open of default-locked table pointed_LBA.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
2025-10-22 14:34:01 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: c8b7005
2025-10-22 14:34:01 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 53.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 0.000 m
Lowest channel frequency: 100.000 MHz
Channel frequency resolution: 1000.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 50
Start date: 2024-07-21T03:35:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 37 antenna positions from measurement set: pointed_LBA.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-22 14:34:02 - INFO - Obs epoch initial LST was 0.0228000546 deg
2025-10-22 14:34:02 - INFO - Setting initial J2000 LST to 359.7122369008 deg
2025-10-22 14:34:02 - INFO - Setting initial mjd to 60512.1493171296
2025-10-22 14:34:02 - INFO - After precession initial latitude of the array is 52.7688502802 deg
2025-10-22 14:34:02 - INFO - Precessing array layout to J2000
2025-10-22 14:34:02 - INFO - Will run with EveryBeam LOFAR primary beam, based on this measurement set:
pointed_LBA.ms
Primary beam is pointed at RA,Dec = 0.0, 53.0 deg
2025-10-22 14:34:02 - 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-22 14:34:02 - INFO - Sky model mapping took 0.0 mins
2025-10-22 14:34:02 - INFO - Have read in 1 components
2025-10-22 14:34:02 - INFO - After cropping there are 1 components
2025-10-22 14:34:02 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-22 14:34:02 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-10-22 14:34:02 - INFO - Reading set 0 sky models
2025-10-22 14:34:02 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:34:02 - INFO - Sending Sky set 0 to the GPU
2025-10-22 14:34:02 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-22 14:34:03 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-22 14:34:03 - INFO - Full run took 0:00:02.004226
2025-10-22 14:34:03 - INFO - WODEN is done. Closing the log. S'later
for pol in ['I', 'Q', 'U', 'V']:
uvfits_name = f"stokes{pol}_everybeam_LOFAR"
XX, XY, YX, YY = read_uvfits(f'{uvfits_name}_band01.uvfits')
##on_cardinal=False means we're using the maths above to recombine the
##XX, XY, YX, YY into Stokes I, Q, U, V
test_stokes_recovery(pol, 'none', atol=5e-2, on_cardinal=False)
# ##uncomment to print out the exact values recovered
# ##pick a random baseline to print, they should all be the sam
# baseline = 0
# ##Supply on_cardinal=False to use the right maths
# recover_I, recover_Q, recover_U, recover_V = convert_inst_to_stokes(XX[baseline], XY[baseline], YX[baseline], YY[baseline],
# on_cardinal=False)
# print(f"Input {pol}, Recover I {recover_I[0].real:.2f}, Q {recover_Q[0].real:.2f}, U {recover_U[0].real:.2f}, V {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
Booya, works perfectly. Now to test the RM recovery.
RM recovery¶
phi_RM, pol_frac = make_RM_skymodel(ra0, dec0)
freq_reso = 0.1e+6
low_freq = 100e+6
high_freq = 130e+6
num_freq_chans = int((high_freq - low_freq) / freq_reso)
primary_beam = "everybeam_LOFAR"
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'--latitude={lofar_lat} --longitude={lofar_long} '
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 --num_threads=1 '
command += f' --station_id=0 --off_cardinal_dipoles '
command += f'--eb_point_to_phase '
command += f'--beam_ms_path=pointed_LBA.ms'
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.6.0')
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful readonly open of default-locked table pointed_LBA.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful read/write open of default-locked table /home/jack-line/software/WODEN_dev/test_installation/everybeam/pointed_rm_source_everybeam_LOFAR_band01.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
2025-10-22 14:34:13 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: c8b7005
2025-10-22 14:34:13 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 53.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 0.000 m
Lowest channel frequency: 100.000 MHz
Channel frequency resolution: 100.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 300
Start date: 2024-07-21T03:35:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 37 antenna positions from measurement set: pointed_LBA.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-22 14:34:13 - INFO - Obs epoch initial LST was 0.0228000546 deg
2025-10-22 14:34:13 - INFO - Setting initial J2000 LST to 359.7122369008 deg
2025-10-22 14:34:13 - INFO - Setting initial mjd to 60512.1493171296
2025-10-22 14:34:13 - INFO - After precession initial latitude of the array is 52.7688502802 deg
2025-10-22 14:34:13 - INFO - Precessing array layout to J2000
2025-10-22 14:34:13 - INFO - Will run with EveryBeam LOFAR primary beam, based on this measurement set:
pointed_LBA.ms
Created the following minimal MS to point the beam:
/home/jack-line/software/WODEN_dev/test_installation/everybeam/pointed_rm_source_everybeam_LOFAR_band01.ms
Primary beam is pointed at RA,Dec = 0.0, 53.0 deg
2025-10-22 14:34:13 - 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-22 14:34:13 - INFO - Sky model mapping took 0.0 mins
2025-10-22 14:34:13 - INFO - Have read in 1 components
2025-10-22 14:34:13 - INFO - After cropping there are 1 components
2025-10-22 14:34:13 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-22 14:34:14 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-10-22 14:34:14 - INFO - Reading set 0 sky models
2025-10-22 14:34:14 - INFO - From sky set 0 thread num 0 reading 1 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:34:14 - INFO - Sending Sky set 0 to the GPU
2025-10-22 14:34:14 - INFO - Sky set 0 has returned from the GPU after 0.3 seconds
2025-10-22 14:34:15 - INFO - Have completed 1 of 1 components calcs (100.0%)
2025-10-22 14:34:15 - INFO - Deleting /home/jack-line/software/WODEN_dev/test_installation/everybeam/pointed_rm_source_everybeam_LOFAR_band01.ms...
2025-10-22 14:34:15 - INFO - Done
2025-10-22 14:34:15 - INFO - Full run took 0:00:02.404859
2025-10-22 14:34:15 - INFO - WODEN is done. Closing the log. S'later
0
XX, XY, YX, YY = read_uvfits('rm_source_everybeam_LOFAR_band01.uvfits')
baseline = 0
recover_I, recover_Q, recover_U, recover_V = convert_inst_to_stokes(XX[baseline], XY[baseline], YX[baseline], YY[baseline],
on_cardinal=False)
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()
test_RM_recovery(uvfits_name, phi_RM, pol_frac, freqs, atol=0.02)
Recovered RM: 50.0 Expected RM: 50
Recovered Pol. Fraction: 0.5 Expected Pol Fraction 0.5
Looks good.
Test unique primary beam per station¶
run_everybeam (which is called by run_everybeam_over_threads) has been setup to run any given number of stations, times, frequencies and directions. Here we’ll test things are returned in the expecting order by running two times and two frequencies. As waiting around is boring, we’ll also run things across multiple threads; this is a mircocosm of what happens in run_woden.py, as we split the sky model reading and beam calculations across threads.
from astropy.time import TimeDelta
from wodenpy.primary_beam.use_everybeam import get_num_stations
##Setup a grid of RA/Dec on the sky
nside=100
radec_reso = 100/nside
##lock frame to first LST, so we can see if beam is locked to RA/Dec
header, wcs = create_WCS(LST_deg, 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)
num_beams = get_num_stations(ms_path)
##Do a couple of frequencies to check mapping is working
all_freqs = np.array([50e+6, 100e+6])
##Do two times, 4 hours apart to check things move on the sky
##Our obs frame is locked to beam phase, so should still be centred, but it'll
##be elongated due to projection effects/pointing
all_times = np.array([observing_time, observing_time + TimeDelta(4*3600, format='sec')])
station_ids = np.arange(num_beams)
##We'll run EveryBeam in parallel to speed things up
num_threads=8
print(LST_deg, np.radians(ra0), np.radians(dec0))
all_jones = run_everybeam_over_threads(num_threads,
ras, decs,
ms_path=ms_path,
times=all_times, freqs=all_freqs,
station_ids=station_ids,
beam_ra0=np.radians(ra0),beam_dec0=np.radians(dec0),
apply_beam_norms=True,
iau_order=True,
parallactic_rotate=True)
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
0.01862197128531748 0.0 0.9250245035569946
OK, now setup code to plot the the absolute values of the first polarisation for the first 36 stations, per time, per freq, to check our outputs.
def plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs):
fig, axs = plt.subplots(6, 6, figsize=(12, 12), layout='constrained',
subplot_kw={'projection': wcs})
for col in range(6):
for row in range(6):
beam_ind = col*6 + 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 < 5:
lon.set_ticklabel_visible(False)
if col == 5:
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)
Cool, we get a different primary beam shape for each station, they’re all normalised to 1. So far so good.
time_ind = 0
freq_ind = 1
plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs)
Yup, if we switch from 50 to 100MHz the beam gets smaller on the sky, happy days.
time_ind = 1
freq_ind = 0
plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs)
Now we’ve gone from beam at zenith to 4 hours off zenith. We’re still plotting with RA/Dec centred on beam, but we can see projection effects so we are defo off zenith.
time_ind = 1
freq_ind = 1
plot_beam_selection(all_jones, time_ind, freq_ind, nside, wcs)
And again changing freq changes size. Things are looking good.
Check positions of components via image¶
Let’s bung a grid of RA/Dec through the whole simulation code, so multiple frequencies and time steps. This should check we’re calling the above beam code correctly and integrating things into the code proper.
from astropy.table import Column, Table
nside=11
half_width = 2
ras = np.linspace(ra0 - half_width, ra0 + half_width, nside)
decs = np.linspace(dec0 - half_width, dec0 + half_width, nside)
ras, decs = np.meshgrid(ras, decs)
ras = ras.flatten()
decs = decs.flatten()
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)
cat_name = 'check_positions.fits'
main_table.write(cat_name, format='fits', overwrite=True)
args = []
primary_beam="everybeam_LOFAR"
low_freq=50e+6
num_freq_chans=100
freq_reso=1e+4
uvfits_name=f"check_positions_{primary_beam}"
##The command to run WODEN
cmd = "run_woden.py "
cmd += f'--ra0={ra0} '
cmd += f'--dec0={dec0} '
cmd += f'--latitude={lofar_lat} '
cmd += f'--longitude={lofar_long} '
cmd += f'--date={date} '
cmd += f'--output_uvfits_prepend={uvfits_name} '
cmd += f'--cat_filename=check_positions.fits '
cmd += f'--primary_beam={primary_beam} '
cmd += f'--lowest_channel_freq={low_freq} '
cmd += f'--freq_res={freq_reso} '
cmd += f'--num_freq_channels={num_freq_chans} '
cmd += f'--band_nums=1 '
cmd += f'--time_res=10 '
cmd += f'--num_time_steps=10 '
cmd += f'--IAU_order '
cmd += f'--off_cardinal_dipoles '
cmd += f'--eb_point_to_phase '
cmd += '--num_threads=8 '
cmd += f'--station_id=0 '
cmd += '--cpu_mode '
cmd += f'--beam_ms_path=pointed_LBA.ms '
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')
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful readonly open of default-locked table pointed_LBA.ms/SPECTRAL_WINDOW: 14 columns, 1 rows
Successful read/write open of default-locked table pointed_LBA.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms: 24 columns, 703 rows
Successful read/write open of default-locked table /home/jack-line/software/WODEN_dev/test_installation/everybeam/pointed_check_positions_everybeam_LOFAR_band01.ms::FIELD: 10 columns, 1 rows
Successful readonly open of default-locked table pointed_LBA.ms/ANTENNA: 10 columns, 37 rows
2025-10-22 14:41:25 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: c8b7005
2025-10-22 14:41:25 - INFO - Input arguments after parsing:
Phase centre: 0.00000, 53.00000 deg
Array central latitude: 52.905 deg
Array central longitude: 6.868 deg
Array central height: 0.000 m
Lowest channel frequency: 50.000 MHz
Channel frequency resolution: 10.000 kHz
Coarse band bandwidth: 0.195 MHz
Num channels per coarse band: 100
Start date: 2024-07-21T03:35:00
Time resolution: 10.0 (s)
Num time steps: 10
Have read 37 antenna positions from measurement set: pointed_LBA.ms
Will write outputs to: /home/jack-line/software/WODEN_dev/test_installation/everybeam
2025-10-22 14:41:25 - INFO - Obs epoch initial LST was 0.0395123880 deg
2025-10-22 14:41:25 - INFO - Setting initial J2000 LST to 359.7288967971 deg
2025-10-22 14:41:25 - INFO - Setting initial mjd to 60512.1493634259
2025-10-22 14:41:25 - INFO - After precession initial latitude of the array is 52.7688494985 deg
2025-10-22 14:41:25 - INFO - Precessing array layout to J2000
2025-10-22 14:41:25 - INFO - Will run with EveryBeam LOFAR primary beam, based on this measurement set:
pointed_LBA.ms
Created the following minimal MS to point the beam:
/home/jack-line/software/WODEN_dev/test_installation/everybeam/pointed_check_positions_everybeam_LOFAR_band01.ms
Primary beam is pointed at RA,Dec = 0.0, 53.0 deg
2025-10-22 14:41:25 - 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-22 14:41:25 - INFO - Sky model mapping took 0.0 mins
2025-10-22 14:41:25 - INFO - Have read in 121 components
2025-10-22 14:41:25 - INFO - After cropping there are 121 components
2025-10-22 14:41:25 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-10-22 14:41:26 - INFO - Running in parallel on CPU mode with 8 threads
There are 1 sets of sky models to run
2025-10-22 14:41:26 - INFO - Reading set 0 sky models
2025-10-22 14:41:26 - INFO - From sky set 0 thread num 0 reading 16 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:41:26 - INFO - From sky set 0 thread num 1 reading 16 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:41:26 - INFO - From sky set 0 thread num 2 reading 16 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:41:26 - INFO - From sky set 0 thread num 3 reading 16 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:41:26 - INFO - From sky set 0 thread num 5 reading 16 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:41:26 - INFO - From sky set 0 thread num 4 reading 16 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:41:26 - INFO - From sky set 0 thread num 7 reading 9 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:41:26 - INFO - From sky set 0 thread num 6 reading 16 points, 0 gauss, 0 shape, 0 shape coeffs
2025-10-22 14:41:27 - INFO - Sending Sky set 0 chunk 0 to the CPU
2025-10-22 14:41:27 - INFO - Sending Sky set 0 chunk 1 to the CPU
2025-10-22 14:41:27 - INFO - Sending Sky set 0 chunk 2 to the CPU
2025-10-22 14:41:27 - INFO - Sending Sky set 0 chunk 3 to the CPU
2025-10-22 14:41:28 - INFO - Sending Sky set 0 chunk 4 to the CPU
2025-10-22 14:41:28 - INFO - Sending Sky set 0 chunk 5 to the CPU
2025-10-22 14:41:28 - INFO - Sending Sky set 0 chunk 6 to the CPU
2025-10-22 14:41:28 - INFO - Sending Sky set 0 chunk 7 to the CPU
2025-10-22 14:41:29 - INFO - Sky set 0 chunk 0 has returned from the CPU after 2.2 seconds
2025-10-22 14:41:29 - INFO - Sky set 0 chunk 1 has returned from the CPU after 2.2 seconds
2025-10-22 14:41:29 - INFO - Sky set 0 chunk 2 has returned from the CPU after 2.2 seconds
2025-10-22 14:41:29 - INFO - Sky set 0 chunk 3 has returned from the CPU after 2.2 seconds
2025-10-22 14:41:30 - INFO - Sky set 0 chunk 7 has returned from the CPU after 1.4 seconds
2025-10-22 14:41:30 - INFO - Sky set 0 chunk 4 has returned from the CPU after 2.2 seconds
2025-10-22 14:41:30 - INFO - Sky set 0 chunk 5 has returned from the CPU after 2.2 seconds
2025-10-22 14:41:30 - INFO - Sky set 0 chunk 6 has returned from the CPU after 2.2 seconds
2025-10-22 14:41:31 - INFO - Have completed 121 of 121 components calcs (100.0%)
2025-10-22 14:41:31 - INFO - Finished all rounds of processing
2025-10-22 14:41:31 - INFO - Deleting /home/jack-line/software/WODEN_dev/test_installation/everybeam/pointed_check_positions_everybeam_LOFAR_band01.ms...
2025-10-22 14:41:31 - INFO - Done
2025-10-22 14:41:31 - INFO - Full run took 0:00:06.663804
2025-10-22 14:41:31 - INFO - WODEN is done. Closing the log. S'later
0
cmd = "woden_uv2ms.py "
cmd += " --uvfits_prepend=check_positions_everybeam_LOFAR_band "
cmd += " --band_nums=1 "
call(cmd, shell=True)
cmd = "wsclean -name check_positions_everybeam_LOFAR -size 2048 2048 -niter 2000 "
cmd += " -auto-threshold 0.5 -auto-mask 3 "
cmd += " -pol I -multiscale -weight briggs 0 -scale 0.002 -j 12 -mgain 0.85 "
cmd += " -no-update-model-required "
cmd += " check_positions_everybeam_LOFAR_band*.ms "
call(cmd, shell=True)
/home/jack-line/software/WODEN_dev/woden_dev/bin/woden_uv2ms.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')
The telescope frame is set to '????', which generally indicates ignorance. Defaulting the frame to 'itrs', but this may lead to other warnings or errors.
The uvw_array does not match the expected values given the antenna positions. The largest discrepancy is 7.1919323818235625 meters. This is a fairly common situation but might indicate an error in the antenna positions, the uvws or the phasing.
The uvw_array does not match the expected values given the antenna positions. The largest discrepancy is 7.1919323818235625 meters. This is a fairly common situation but might indicate an error in the antenna positions, the uvws or the phasing.
Writing in the MS file that the units of the data are uncalib, although some CASA process will ignore this and assume the units are all in Jy (or may not know how to handle data in these units).
WSClean version 3.4 (2023-10-11)
This software package is released under the GPL version 3.
Author: André Offringa (offringa@gmail.com).
No corrected data in first measurement set: tasks will be applied on the data column.
=== IMAGING TABLE ===
# Pol Ch JG ²G FG FI In Freq(MHz)
| Independent group:
+-+-J- 0 I 0 0 0 0 0 0 50-51 (100)
Reordering check_positions_everybeam_LOFAR_band01.ms into 1 x 1 parts.
Reordering: 0%....10%....20%....30%....40%....50%....60%....70%....80%....90%....100%
Initializing model visibilities: 0%....10%....20%....30%....40%....50%....60%....70%....80%....90%....100%
== Constructing PSF ==
Precalculating weights for Briggs'(0) weighting...
Opening reordered part 0 spw 0 for check_positions_everybeam_LOFAR_band01.ms
Detected 62.7 GB of system memory, usage not limited.
Opening reordered part 0 spw 0 for check_positions_everybeam_LOFAR_band01.ms
Determining min and max w & theoretical beam size... DONE (w=[0.000267848:104.927] lambdas, maxuvw=14348 lambda)
Theoretic beam = 14.38''
Small inversion enabled, but inversion resolution already smaller than beam size: not using optimization.
Loading data in memory...
Gridding 6660 rows...
Gridded visibility count: 663663
Fitting beam... major=26.65'', minor=10.44'', PA=-114.18 deg, theoretical=14.38''.
Writing psf image... DONE
== Constructing image ==
Opening reordered part 0 spw 0 for check_positions_everybeam_LOFAR_band01.ms
Loading data in memory...
Gridding 6660 rows...
Gridded visibility count: 663663
Writing dirty image...
== Deconvolving (1) ==
Estimated standard deviation of background noise: 62.29 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.136688
- Scale 8, bias factor=1.7, psfpeak=0.21916, gain=0.456287, kernel peak=0.047099
- Scale 16, bias factor=2.8, psfpeak=0.121282, gain=0.824527, kernel peak=0.0142695
- Scale 32, bias factor=4.6, psfpeak=0.0710628, gain=1.40721, kernel peak=0.00396542
- Scale 64, bias factor=7.7, psfpeak=0.0466004, gain=2.1459, kernel peak=0.00104807
- Scale 128, bias factor=12.9, psfpeak=0.0291355, gain=3.43224, kernel peak=0.000269603
- Scale 256, bias factor=21.4, psfpeak=0.0116829, gain=8.55949, kernel peak=6.83828e-05
- Scale 511, bias factor=35.7, psfpeak=0.00346552, gain=28.8557, kernel peak=1.72211e-05
- Scale 1022, bias factor=59.5, psfpeak=0.00113027, gain=88.4741, kernel peak=4.32226e-06
RMS per scale: {0: 72.93 mJy, 8: 36.92 mJy, 16: 28.86 mJy, 32: 20.64 mJy, 64: 13.36 mJy, 128: 8.2 mJy, 256: 4.64 mJy, 511: 4.1 mJy, 1022: 3.31 mJy}
Starting multi-scale cleaning. Start peak=1.01 Jy, major iteration threshold=186.87 mJy (final)
Iteration 31, scale 0 px : 804.48 mJy at 1267,825
Iteration 150, scale 0 px : 646.76 mJy at 1395,427
Iteration 331, scale 0 px : 517 mJy at 666,1227
Iteration 603, scale 0 px : 416.49 mJy at 1273,225
Iteration 914, scale 0 px : 338.5 mJy at 565,2029
Iteration 1252, scale 32 px : 283.93 mJy at 1986,127
Iteration 1267, scale 0 px : 272.99 mJy at 1150,25
Iteration 1617, scale 64 px : 290.82 mJy at 74,925
Iteration 1631, scale 64 px : -231.22 mJy at 138,1469
Subminor loop is near minor loop threshold. Initiating countdown.
(12) Iteration 1719, scale 32 px : -249.09 mJy at 16,1271
Iteration 1742, scale 1022 px : 244.68 mJy at 1411,1535
Iteration 1749, scale 0 px : 223.83 mJy at 906,1424
(11) Iteration 2000, scale 511 px : 232.38 mJy at 256,256
Cleaning finished because maximum number of iterations was reached.
Auto-masking threshold reached; continuing next major iteration with deeper threshold and mask.
Maximum number of minor deconvolution iterations was reached: not continuing deconvolution.
Assigning from 1 to 1 channels...
Writing model image...
== Converting model image to visibilities ==
Opening reordered part 0 spw 0 for check_positions_everybeam_LOFAR_band01.ms
Predicting 6660 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for check_positions_everybeam_LOFAR_band01.ms
Loading data in memory...
Gridding 6660 rows...
Gridded visibility count: 663663
1 major iterations were performed.
Rendering sources to restored image (beam=10.44''-26.65'', PA=-114.18 deg)... DONE
Writing restored image... DONE
Multi-scale cleaning summary:
- Scale 0 px, nr of components cleaned: 1853 (69.51 Jy)
- Scale 8 px, nr of components cleaned: 0 (0 Jy)
- Scale 16 px, nr of components cleaned: 0 (0 Jy)
- Scale 32 px, nr of components cleaned: 38 (140.08 mJy)
- Scale 64 px, nr of components cleaned: 102 (2.47 Jy)
- Scale 128 px, nr of components cleaned: 0 (0 Jy)
- Scale 256 px, nr of components cleaned: 0 (0 Jy)
- Scale 511 px, nr of components cleaned: 0 (0 Jy)
- Scale 1022 px, nr of components cleaned: 7 (2.29 Jy)
Total: 2000 components (74.41 Jy)
Inversion: 00:00:01.732557, prediction: 00:00:00.401218, deconvolution: 00:00:07.957851
Cleaning up temporary files...
0
with fits.open('check_positions_everybeam_LOFAR-image.fits') as hdu:
image = np.squeeze(hdu[0].data)
wcs = WCS(hdu[0].header).celestial
fig, axs = plt.subplots(1, 1, figsize=(12, 8), subplot_kw={'projection': wcs})
im = axs.imshow(image, origin='lower', cmap='Blues_r', vmin=-0.1, vmax=0.5)
plt.colorbar(im, ax=axs, label='Jy/beam')
# half_width = 600
# axs.set_xlim(1024-half_width, 1024+half_width)
# axs.set_ylim(1024-half_width, 1024+half_width)
plt.grid(alpha=0.5)
plt.show()
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 60512.149363 from DATE-OBS'. [astropy.wcs.wcs]
Fairly ratty image, but we defo have a grid of point source aligned with the celestial sphere, so we’re looking fine. I think lba.MS is a minimal example of the LBA array so the \(uv\) coverage probably isn’t great. Actually we can have a look at that:
with fits.open('check_positions_everybeam_LOFAR_band01.uvfits') as hdu:
uu = (hdu[0].data['UU']*c)/1e+3
vv = (hdu[0].data['VV']*c)/1e+3
fig, axs = plt.subplots(1, 1, figsize=(8, 8))
axs.plot(uu, vv, 'k.', ms=1)
axs.plot(-uu, -vv, 'k.', ms=1)
axs.set_xlabel('U (km)')
axs.set_ylabel('V (km)')
plt.show()
Yeah pretty patchy coverage. Amazing resolution though, blimey!