Polarisation examples¶
Below are a couple of examples of how to use WODEN to simulation polarisation data. Both examples include details on how to generate the FITS-style sky model, how t run the simulation, and check the results. This documentation is simply the notebook found in WODEN/examples/polarisation/polarisation_examples.ipynb. If you want to run it, you’ll need a working WODEN installation, WSClean to do the imaging, and the pygdsm package (I install via pip install git+https://github.com/telegraphic/pygdsm) to generate the sky model.
Big shout out to Emil for all the help in getting Jack to (kind of) understanding polarisation, and for the RM synthesis code.
Apologies, I haven’t worked out how to get outputs from cells to be scrollable in the online documentation, so please just scroll past the reams of output text to get to plots.
Southern Hotspot of PKS J0636−2036¶
This is an example of a real source, details of which can be found in O’Sullivan et al. (2018). Assuming that the southern lobe is a single point source, at 216MHz, the component has the parameters (from Table 3 of this paper):
Stokes \(I = 13.9\) Jy (at reference 216MHz)
Polarisation fraction \(\Pi = 0.092\)
Spectral index \(\alpha = -0.815\)
Faraday rotation \(\phi_\textrm{RM} = 50.2 \) rad \(m^{-2}\)
Ignoring the jump in frequency and instrument in the measurement (it’s an example so we can do what we want) we can also say Intrinsic rotation angle \(\chi_0 = -71.1^{\circ}\)
Furthermore, to make this more useful as an example, we’ll say the source actually has a Stokes \(V\) component that bounces around, so we can look at list-type polarisation data.
First, we’ll make a FITS table with these parameters. Then we’ll stick it through WODEN, and process the outputs.
from astropy.table import Table, Column
from astropy.coordinates import SkyCoord
import numpy as np
from subprocess import call
from astropy.coordinates import EarthLocation
from astropy import units as u
from astropy.time import Time
from astropy.io import fits
import matplotlib.pyplot as plt
from pygdsm import GlobalSkyModel16
import healpy as hp
from astropy.constants import c
C = c.value
##Stick the source at location of it's name J0636−2036, good enough
coord = SkyCoord(ra='6h36m', dec="-20d36m")
ra0 = coord.ra.deg
dec0 = coord.dec.deg
##These are the parameters for the source
alpha = -0.815
##hyperdrive/WODEN/LoBES catalogues are all referenced to 200MHz, so we need
##to scale the stokesI flux
stokesI = 13.9 * (200/216)**(alpha)
pol_frac = 0.092
phi_RM = 50.2
chi_0 = np.radians(71.1)
##The LoBES/hyperdrive/WODEN catalogues need a unique source ID (UNQ_SOURCE_ID),
##a then edge component of that source needs a name (NAME)
c_ids = Column(data=np.array(['PKS_J0636_source']), name='UNQ_SOURCE_ID', dtype='|S20')
c_names = Column(data=np.array(['PKS_J0636_source_C00']), name='NAME', dtype='|S20')
##Component position
c_ras = Column(data=np.array([ra0]), name='RA')
c_decs = Column(data=np.array([dec0]), name='DEC')
##This says we have a point source
c_comp_types = Column(data=np.array(['P'], 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'], dtype='|S3'), name="MOD_TYPE", dtype='|S3')
##Stokes I power law parameters
c_stokes_I_ref = Column(data=np.array([stokesI]), name='NORM_COMP_PL')
c_stokes_SI = Column(data=np.array([alpha]), name='ALPHA_PL')
##This says we are using a polarisation fraction linear polarisation model
c_lin_mod_type = Column(data=np.array(['pf'], dtype='|S3'), name='LIN_MOD_TYPE', dtype='|S3')
##linear polarisation parameters
c_lin_pol_frac = Column(data=np.array([pol_frac]), name='LIN_POL_FRAC')
c_lin_pol_angle = Column(data=np.array([chi_0]), name='INTR_POL_ANGLE')
c_rm = Column(data=np.array([phi_RM]), name='RM')
##This says we are using a list-type circular polarisation model. The list info
##goes into another table, see below.
c_v_mod_type = Column(data=np.array(['nan'], dtype='|S3'), name='V_MOD_TYPE', dtype='|S3')
##Stick all the columns in a list
columns = [c_ids, c_names, c_ras, c_decs, c_comp_types, c_mod_types, c_stokes_I_ref, c_stokes_SI, c_lin_mod_type,c_lin_pol_frac, c_lin_pol_angle, c_rm, c_v_mod_type]
##Make the main table
main_table = Table(columns)
##Now, make a table for Stokes V list-type SED. This goes in it's own table,
##where the columns have the frequency in the title. We dupliate the NAME
##column so that we can link the two tables
v_cols = [c_names]
v_freqs = np.arange(150, 250, 0.2) ##MHz
v_fluxes = np.random.uniform(-0.5, 0.5, len(v_freqs))
for freq, flux in zip(v_freqs, v_fluxes):
col_name = f'V_INT_FLX{freq:.1f}'
v_cols.append(Column(data=np.array([flux]), name=col_name))
v_table = Table(v_cols)
##Create an HDU list from the Tables, so we can easily write it to a FITS file
hdu_list = [fits.PrimaryHDU(),
fits.table_to_hdu(main_table),
fits.table_to_hdu(v_table)]
##Name the tables so WODEN can understand which is which
hdu_list[1].name = 'MAIN'
hdu_list[2].name = 'V_LIST_FLUXES'
##convert and write out to a FITS table file
hdu_list = fits.HDUList(hdu_list)
cat_name = 'PKS_J0636.fits'
hdu_list.writeto(cat_name, overwrite=True)
Now that we have our single source, let’s setup a simple simulation command. For this example, we’re just testing the polarisation reponse of WODEN, so make a very very basic simulatin. We’ll run without a primary beam, put the source at phase centre, and run with just a few baselines via a custom array layout.
Importantly, we’ll set the --IAU_order flag, which means that XX = north-south, YY = east-west. This is the IAU convention for Stokes parameters.
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 = 10
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)
##stick our array in the MWA location.
mwa_location = EarthLocation(lat=-26.703319405555554*u.deg,
lon=116.67081523611111*u.deg,
height=377.827)
##pick a time/date that sticks our source overhead
date = "2024-07-05T04:00:00"
observing_time = Time(date, scale='utc', location=mwa_location)
##Grab the LST
LST = observing_time.sidereal_time('mean')
LST_deg = LST.value*15
print(f"LST: {LST_deg}, RA: {ra0}")
##parameters for the simulation; setting the date sets an LST of about 0 deg
##to match the phase centre
uvfits_name = "PKS_J0636"
freq_reso = 80e+3
low_freq = 180e+6
high_freq = 210e+6
num_freq_chans = int((high_freq - low_freq) / freq_reso)
##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} --array_layout={array_name} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam=none '
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 += '--num_threads=1'
##use subprocess to run the command
call(command, shell=True)
LST: 100.3181217509486, RA: 98.99999999999999
/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-14 17:11:17 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-14 17:11:17 - INFO - Input arguments after parsing:
Phase centre: 99.00000, -20.60000 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: 80.000 kHz
Coarse band bandwidth: 1.280 MHz
Num channels per coarse band: 375
Start date: 2024-07-05T04:00:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 10 antenna positions from array layout file: eg_array.txt
Will write outputs to: /home/jack-line/software/WODEN_dev/examples/polarisation
2025-03-14 17:11:18 - INFO - Obs epoch initial LST was 100.3222998343 deg
2025-03-14 17:11:18 - INFO - Setting initial J2000 LST to 100.0765134441 deg
2025-03-14 17:11:18 - INFO - Setting initial mjd to 60496.1666782406
2025-03-14 17:11:18 - INFO - After precession initial latitude of the array is -26.6814845740 deg
2025-03-14 17:11:18 - INFO - Precessing array layout to J2000
2025-03-14 17:11:18 - INFO - No primary beam selected, no beam attenuation will be applied.
2025-03-14 17:11:18 - INFO - Doing the initial mapping of sky model
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
2025-03-14 17:11:18 - INFO - Sky model mapping took 0.0 mins
2025-03-14 17:11:18 - INFO - Have read in 3 components
2025-03-14 17:11:18 - INFO - After cropping there are 3 components
2025-03-14 17:11:18 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-14 17:11:18 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-03-14 17:11:18 - INFO - Reading set 0 sky models
2025-03-14 17:11:18 - INFO - From sky set 0 thread num 0 reading 3 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 17:11:18 - INFO - Sending Sky set 0 to the GPU
2025-03-14 17:11:19 - INFO - Sky set 0 has returned from the GPU after 0.1 seconds
2025-03-14 17:11:19 - INFO - Have completed 3 of 3 components calcs (100.0%)
2025-03-14 17:11:19 - INFO - Full run took 0:00:01.302348
2025-03-14 17:11:19 - INFO - WODEN is done. Closing the log. S'later
0
Righto, now we have an output, let’s read it in, convert it to Stokes parameters, and plot the SED.
def read_uvfits(uvfits_name):
"""Real shorthand function to read in visibilities from a WODEN UVFITS file"""
with fits.open(uvfits_name) as hdus:
data = np.squeeze(hdus[0].data.data)
##Yup, this is the order of things in a UVFITS file
XX = data[:, :, 0, 0] + 1j*data[:, :, 0, 1]
YY = data[:, :, 1, 0] + 1j*data[:, :, 1, 1]
XY = data[:, :, 2, 0] + 1j*data[:, :, 2, 1]
YX = data[:, :, 3, 0] + 1j*data[:, :, 3, 1]
return XX, XY, YX, YY
XX, XY, YX, YY = read_uvfits(f'{uvfits_name}_band01.uvfits')
##pick a random baseline to plot, they should all be the same
baseline = np.random.randint(0, XX.shape[0])
print(baseline)
##converts instrumental pols into Stokes params. We didn't use a primary beam
##at we're right at zenith so this is good enough without corrections
recover_I = 0.5*(XX[baseline] + YY[baseline])
recover_Q = 0.5*(XX[baseline] - YY[baseline])
recover_U = 0.5*(XY[baseline] + YX[baseline])
recover_V = -0.5j*(XY[baseline] - YX[baseline])
##these are the freqs we put into the simulation
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()
15
Lovely, we have some nice sinusoidal Q/U lines, a power-law Stokes I, and a noisy-like Stokes V. Let’s see if we can recover the RM from this data. All credit to Emil for the RM synthesis code.
def getFDF(dataQ, dataU, freqs, startPhi, stopPhi, dPhi, dType='float32'):
"""
# Perform RM-synthesis on Stokes Q and U data
#
# dataQ, dataU and freqs - contains the Q/U data at each frequency (in Hz) measured.
# startPhi, dPhi - the starting RM (rad/m^2) and the step size (rad/m^2)
Author: Emil Lenc
"""
# Calculate the RM sampling
phiArr = np.arange(startPhi, stopPhi, dPhi)
# Calculate the frequency and lambda sampling
lamSqArr = np.power(C / np.array(freqs), 2.0)
# Calculate the dimensions of the output RM cube
nPhi = len(phiArr)
# Initialise the complex Faraday Dispersion Function (FDF)
FDF = np.ndarray((nPhi), dtype='complex')
# Assume uniform weighting
wtArr = np.ones(len(lamSqArr), dtype=dType)
K = 1.0 / np.nansum(wtArr)
# Get the weighted mean of the LambdaSq distribution (B&dB Eqn. 32)
lam0Sq = K * np.nansum(lamSqArr)
# Mininize the number of inner-loop operations by calculating the
# argument of the EXP term in B&dB Eqns. (25) and (36) for the FDF
a = (-2.0 * 1.0j * phiArr)
b = (lamSqArr - lam0Sq)
arg = np.exp( np.outer(a, b) )
# Create a weighted complex polarised surface-brightness cube
# i.e., observed polarised surface brightness, B&dB Eqns. (8) and (14)
Pobs = (np.array(dataQ) + 1.0j * np.array(dataU))
# Calculate the Faraday Dispersion Function
# B&dB Eqns. (25) and (36)
FDF = K * np.nansum(Pobs * arg, 1)
return FDF, phiArr
def findpeaks(freqs, fdf, phi, rmsf, rmsfphi, nsigma):
"""Find peaks in the FDF
Author: Emil Lenc"""
# Create the Gaussian filter for reconstruction
lam2 = (C / freqs) ** 2.0
lam02 = np.mean(lam2)
minl2 = np.min(lam2)
maxl2 = np.max(lam2)
width = (2.0 * np.sqrt(3.0)) / (maxl2 - minl2)
Gauss = np.exp((-rmsfphi ** 2.0) / (2.0 * ((width / 2.355) ** 2.0)))
components = np.zeros((len(phi)), np.float32)
peaks = []
phis = []
std = 0.0
rmsflen = int((len(rmsf) - 1) / 2)
fdflen = len(phi) + rmsflen
while True:
std = np.std(np.abs(fdf))
peak1 = np.max(np.abs(fdf))
pos1 = np.argmax(np.abs(fdf))
val1 = phi[pos1]
if peak1 < nsigma * std :
break
fdf -= rmsf[rmsflen - pos1:fdflen - pos1] * fdf[pos1]
peaks.append(peak1)
phis.append(val1)
components[pos1] += peak1
fdf += np.convolve(components, Gauss, mode='valid')
return phis, peaks, std
startPhi = -100.0
dPhi = 1.0
stopPhi = -startPhi+dPhi
lambda2 = np.power(C / np.array(freqs), 2.0)
df = []
for f in range(1, len(freqs)):
df.append(freqs[f] - freqs[f-1])
chanBW = np.min(np.array(df))
fmin = np.min(freqs)
fmax = np.max(freqs)
bw = fmax - fmin
dlambda2 = np.power(C / fmin, 2) - np.power(C / (fmin + chanBW), 2)
Dlambda2 = np.power(C / fmin, 2) - np.power(C / (fmin + bw), 2)
phimax = np.sqrt(3) / dlambda2
dphi = 2.0 * np.sqrt(3) / Dlambda2
phiR = dphi / 5.0
Nphi = 2 * phimax / phiR
##These are things that mean something to polarisation peoples
print("Input resolutions going into RM-synthesis:")
print("\tFrequency range: %7.3f MHz - %7.3f MHz" %(fmin / 1.0e6, (fmin + bw) / 1.0e6))
print("\tBandwidth: %7.3f MHz" %(bw / 1.0e6))
print("\tChannel width: %.1f KHz" %(chanBW / 1.0e3))
print("\tdlambda2: %7.3f" %(dlambda2))
print("\tDlambda2: %7.3f" %(Dlambda2))
print("\tphimax: %7.3f" %(phimax))
print("\tdphi: %7.3f" %(dphi))
print("\tphiR: %7.3f" %(phiR))
print("\tNphi: %7.3f" %(Nphi))
fwhm = dphi
# Determine the FDF using the q, u and ferquency values read from the file.
dirty, phi = getFDF(recover_Q, recover_U, freqs, startPhi, stopPhi, dPhi)
FDFqu, phi = getFDF(recover_Q, recover_U, freqs, startPhi, stopPhi, dPhi)
rstartPhi = startPhi * 2
rstopPhi = stopPhi * 2 - dPhi
RMSF, rmsfphi = getFDF(np.ones((len(recover_Q))), np.zeros((len(recover_Q))), freqs, rstartPhi, rstopPhi, dPhi)
phis, peaks, sigma = findpeaks(np.array(freqs), FDFqu, phi, RMSF, rmsfphi, 6.0)
snr = peaks / sigma
phierr = fwhm / (2 * snr)
# print(phis, phierr, peaks, sigma)
imean = np.mean(np.real(recover_I))
print("Recovered paramaters:")
print("\tS = %.3f" %(imean))
if len(peaks) > 0:
print("\tPol frac from recovered peak = %.3f%%" %(100.0 * peaks[0] / imean))
##What is the polarised flux as a function of frequency
# polarised_flux = np.sqrt(recover_Q**2 + recover_U**2)
polarised_flux = np.abs(recover_Q + 1j*recover_U)
recovered_pol_frac = (polarised_flux.real / recover_I.real)
print(f'\tMean recovered pol frac: {np.mean(recovered_pol_frac):.3f}, Expected: {pol_frac:.3f}')
fig, axs = plt.subplots(1, 1, figsize=(6,4))
axs.plot(phi, np.real(FDFqu), 'k-', label='FDF')
axs.set_xlabel("${\phi}$ (rad m$^{-2}$)")
axs.set_ylabel("Polarised Flux (Jy RMSF$^{-1}$)")
axs.set_xlim(-100, 100)
axs.axvline(x=phi_RM, color='C1', linestyle='--', label='Expected RM')
axs.legend()
plt.show()
Input resolutions going into RM-synthesis:
Frequency range: 180.000 MHz - 209.920 MHz
Bandwidth: 29.920 MHz
Channel width: 80.0 KHz
dlambda2: 0.002
Dlambda2: 0.734
phimax: 702.920
dphi: 4.717
phiR: 0.943
Nphi: 1490.189
Recovered paramaters:
S = 15.133
Pol frac from recovered peak = 9.167%
Mean recovered pol frac: 0.092, Expected: 0.092
Using Emil’s magic RM code, we can see a nice strong peak at the expected RM of 50.2 rad/m^2, and we recover the expected polarisation fraction of 9.2%.
Multiple components same direction¶
Now let’s try sticking a few different components in the same direction, to see if we can disentangle the RMs.
Make a new sky model:
num_comps = 3
##Make some simple component params. Keep everything the same, but vary RM
alpha = np.full(num_comps, -0.815)
stokesI = np.full(num_comps, 1.0)
pol_frac = np.full(num_comps, 0.3)
phi_RM = np.array([10, 50, 90])
chi_0 = np.full(num_comps, 0.0)
##The LoBES/hyperdrive/WODEN catalogues need a unique source ID (UNQ_SOURCE_ID),
##a then edge component of that source needs a name (NAME)
c_ids = Column(data=np.full(num_comps, 'multi_comp', dtype='|S20'), name='UNQ_SOURCE_ID', dtype='|S20')
c_names = Column(data=[f'multi_comp_C{comp:02d}' for comp in range(num_comps)], name='NAME', dtype='|S20')
##Component position
c_ras = Column(data=np.full(num_comps, ra0), name='RA')
c_decs = Column(data=np.full(num_comps, dec0), name='DEC')
##This says we have a point source
c_comp_types = Column(data=np.full(num_comps, 'P', dtype='|S1'), name="COMP_TYPE", dtype='|S1')
##This says we have a Stokes I power-law SED
c_mod_types = Column(data=np.full(num_comps, 'pl', dtype='|S3'), name="MOD_TYPE", dtype='|S3')
##Stokes I power law parameters
c_stokes_I_ref = Column(data=stokesI, name='NORM_COMP_PL')
c_stokes_SI = Column(data=alpha, name='ALPHA_PL')
##This says we are using a polarisation fraction linear polarisation model
c_lin_mod_type = Column(data=np.full(num_comps, 'pf', dtype='|S3'), name='LIN_MOD_TYPE', dtype='|S3')
##linear polarisation parameters
c_lin_pol_frac = Column(data=pol_frac, name='LIN_POL_FRAC')
c_lin_pol_angle = Column(data=chi_0, name='INTR_POL_ANGLE')
c_rm = Column(data=phi_RM, name='RM')
##Stick all the columns in a list
columns = [c_ids, c_names, c_ras, c_decs, c_comp_types, c_mod_types, c_stokes_I_ref, c_stokes_SI, c_lin_mod_type,c_lin_pol_frac, c_lin_pol_angle, c_rm]
##Make a table
table = Table(columns)
##write the table
cat_name = 'multi_comp.fits'
table.write(cat_name, overwrite=True)
Run it through WODEN, convert outputs into I,Q,U,V, have a look:
##The command to run WODEN
command = f'run_woden.py --ra0={ra0} --dec0={dec0} --array_layout={array_name} '
command += f'--date={date} --output_uvfits_prepend={uvfits_name} '
command += f'--cat_filename={cat_name} --primary_beam=none '
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'
##use subprocess to run the command
call(command, shell=True)
XX, XY, YX, YY = read_uvfits(f'{uvfits_name}_band01.uvfits')
##pick a random baseline to plot, they should all be the same
baseline = np.random.randint(0, XX.shape[0])
print(baseline)
##converts instrumental pols into Stokes params. We didn't use a primary beam
##at we're right at zenith so this is good enough without corrections
recover_I = 0.5*(XX[baseline] + YY[baseline])
recover_Q = 0.5*(XX[baseline] - YY[baseline])
recover_U = 0.5*(XY[baseline] + YX[baseline])
recover_V = -0.5j*(XY[baseline] - YX[baseline])
##these are the freqs we put into the simulation
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()
/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-14 17:11:46 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-14 17:11:46 - INFO - Input arguments after parsing:
Phase centre: 99.00000, -20.60000 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: 80.000 kHz
Coarse band bandwidth: 1.280 MHz
Num channels per coarse band: 375
Start date: 2024-07-05T04:00:00
Time resolution: 2.0 (s)
Num time steps: 1
Have read 10 antenna positions from array layout file: eg_array.txt
Will write outputs to: /home/jack-line/software/WODEN_dev/examples/polarisation
2025-03-14 17:11:47 - INFO - Obs epoch initial LST was 100.3222998343 deg
2025-03-14 17:11:47 - INFO - Setting initial J2000 LST to 100.0765134441 deg
2025-03-14 17:11:47 - INFO - Setting initial mjd to 60496.1666782406
2025-03-14 17:11:47 - INFO - After precession initial latitude of the array is -26.6814845740 deg
2025-03-14 17:11:47 - INFO - Precessing array layout to J2000
2025-03-14 17:11:47 - INFO - No primary beam selected, no beam attenuation will be applied.
2025-03-14 17:11:47 - 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-14 17:11:47 - INFO - Sky model mapping took 0.0 mins
2025-03-14 17:11:47 - INFO - Have read in 3 components
2025-03-14 17:11:47 - INFO - After cropping there are 3 components
2025-03-14 17:11:47 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_double.so
2025-03-14 17:11:47 - INFO - Running in serial on GPU.
Will read sky model using 1 threads
There are 1 sets of sky models to run
2025-03-14 17:11:47 - INFO - Reading set 0 sky models
2025-03-14 17:11:47 - INFO - From sky set 0 thread num 0 reading 3 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 17:11:48 - INFO - Sending Sky set 0 to the GPU
2025-03-14 17:11:48 - INFO - Sky set 0 has returned from the GPU after 0.1 seconds
2025-03-14 17:11:48 - INFO - Have completed 3 of 3 components calcs (100.0%)
2025-03-14 17:11:48 - INFO - Full run took 0:00:01.249432
2025-03-14 17:11:48 - INFO - WODEN is done. Closing the log. S'later
24
Wiggly! Alrighty, what do we have in RM space?
# Determine the FDF using the q, u and ferquency values read from the file.
FDFqu, phi = getFDF(recover_Q, recover_U, freqs, startPhi, stopPhi, dPhi)
fig, axs = plt.subplots(1, 1, figsize=(6,4))
axs.plot(phi, np.real(FDFqu), 'k-', label='FDF')
axs.set_xlabel("${\phi}$ (rad m$^{-2}$)")
axs.set_ylabel("Polarised Flux (Jy RMSF$^{-1}$)")
for phi in phi_RM:
axs.axvline(x=phi, color='C1', linestyle='--', label='Expected RM')
axs.legend()
plt.show()
I don’t know nothing about polarisation, so I don’t know if a negative peak in the FDF is bad. But we have three peaks where we expected them to be, which is nice.
Linear polarisation diffuse sky¶
In this example, we’ll make a crazy made-up sky. Then we’ll push it through an instrument which has symmetric X and Y primary beams, which means we can make I/Q/U/V images without beam correction (as there won’t be direction dependent differences).
First up, make the sky model. This will be a little involved as we’re going to make the Stokes I, linear polarised, and Stokes V sky separately. We’ll make Stokes I simple point sources. The linear sky will be the Haslam map (via the pygdsm package), and the Stokes V sky will be a bunch of Gaussian blobs. Everything will have a simple power-law SED.
##We're going to use the EDA2 array
##layout, so we can get away with low-resolution sky model. So set the
##healpixel nside to 256
nside = 256
num_diffuse = hp.nside2npix(nside)
##chuck in 1000 each for stokes I and V
num_stokesI = 1000
num_stokesV = 1000
total_num_comps = num_diffuse + num_stokesI + num_stokesV
##first, the diffuse sky linear polarised sky. Generate it at 200MHz,
##convert it to the nside we want, and then convert it to Jy
gsm_2016 = GlobalSkyModel16(freq_unit='MHz', data_unit='MJysr')
output = gsm_2016.generate(200)
map_nside = hp.npix2nside(len(output))
output *= hp.nside2pixarea(map_nside)*1e+6
##Set some power law properties, and set the RM to 2 ra/m^2 everywhere
diffuse_ref_flux = hp.ud_grade(output, nside)
diffuse_SI = np.full(num_diffuse, -2.7)
diffuse_rm = np.full(num_diffuse, 2)
##coords are natively in galactic, so we need to convert to ICRS
l, b = hp.pix2ang(nside,np.arange(hp.nside2npix(nside)), lonlat=True)
gal_coords = SkyCoord(l*u.deg, b*u.deg, frame='galactic')
diff_ras = gal_coords.icrs.ra.value
diff_decs = gal_coords.icrs.dec.value
##next, make a simple point source stokes I sky
##make 'em power laws
comp_ras = np.random.uniform(0, 360, num_stokesI)
comp_decs = np.random.uniform(-90.0, 90.0, num_stokesI)
comp_ref_fluxes = np.random.uniform(1, 10, num_stokesI)
comp_ref_SIs = np.random.uniform(-1, 0.2, num_stokesI)
##finally, make some weird large gaussian source stokes V sky
##also make 'em power laws
v_ras = np.random.uniform(0, 360, num_stokesV)
v_decs = np.random.uniform(-90.0, 90.0, num_stokesV)
v_ref_fluxes = np.random.uniform(1, 50, num_stokesV)
v_ref_SIs = np.random.uniform(-1, 0.2, num_stokesV)
##these are the gaussian params; all these sizes are in degrees
v_major_axes = np.random.uniform(10, 12, num_stokesV)
v_minor_axes = np.random.uniform(0.25, 0.5, num_stokesV)
v_pas = np.random.uniform(0,360, num_stokesV)
##now we take all these individual parameters and stick them in a table
all_stokes_I_ref = np.zeros(total_num_comps)
all_stokes_I_ref[:num_stokesI] = comp_ref_fluxes
all_stokes_SI = np.zeros(total_num_comps)
all_stokes_SI[:num_stokesI] = comp_ref_SIs
all_lin_mod_type = np.full(total_num_comps, '', dtype='|S3')
all_lin_mod_type[num_stokesI:num_stokesI+num_diffuse] = 'pl'
all_linpol_ref = np.zeros(total_num_comps)
all_linpol_ref[num_stokesI:num_stokesI+num_diffuse] = diffuse_ref_flux
all_linpol_SI = np.zeros(total_num_comps)
all_linpol_SI[num_stokesI:num_stokesI+num_diffuse] = diffuse_SI
all_rm = np.zeros(total_num_comps)
all_rm[num_stokesI:num_stokesI+num_diffuse] = diffuse_rm
all_v_mod_type = np.full(total_num_comps, '', dtype='|S3')
all_v_mod_type[-num_stokesV:] = 'pl'
all_v_ref = np.zeros(total_num_comps)
all_v_ref[-num_stokesV:] = v_ref_fluxes
all_v_SI = np.zeros(total_num_comps)
all_v_SI[-num_stokesV:] = v_ref_SIs
all_majors = np.full(total_num_comps, np.nan)
all_majors[-num_stokesV:] = v_major_axes
all_minors = np.full(total_num_comps, np.nan)
all_minors[-num_stokesV:] = v_minor_axes
all_pas = np.full(total_num_comps, np.nan)
all_pas[-num_stokesV:] = v_pas
##make everything except the stokes V sources gaussian
all_comp_types = np.full(total_num_comps, 'P', dtype='|S1')
all_comp_types[-num_stokesV:] = 'G'
##create Columns to go into the table
c_ids = Column(data=np.array(['{:08d}'.format(i) for i in range(total_num_comps)]), name='UNQ_SOURCE_ID', dtype='|S20')
c_names = Column(data=np.array(['{:08d}_C00'.format(i) for i in range(total_num_comps)]), name='NAME', dtype='|S20')
c_ras = Column(data=np.concatenate([comp_ras, diff_ras, v_ras]), name='RA')
c_decs = Column(data=np.concatenate([comp_decs, diff_decs, v_decs]), name='DEC')
c_comp_types = Column(data=all_comp_types, name="COMP_TYPE", dtype='|S1')
c_mod_types = Column(data=np.full(total_num_comps, 'pl', dtype='|S3'), name="MOD_TYPE", dtype='|S3')
c_stokes_I_ref = Column(data=all_stokes_I_ref, name='NORM_COMP_PL')
c_stokes_SI = Column(data=all_stokes_SI, name='ALPHA_PL')
c_lin_mod_type = Column(data=all_lin_mod_type, name='LIN_MOD_TYPE', dtype='|S3')
c_linpol_ref = Column(data=all_linpol_ref, name='LIN_NORM_COMP_PL')
c_linpol_SI = Column(data=all_linpol_SI, name='LIN_ALPHA_PL')
c_rm = Column(data=all_rm, name='RM')
c_v_mod_type = Column(data=all_v_mod_type, name='V_MOD_TYPE', dtype='|S3')
c_v_ref = Column(data=all_v_ref, name='V_NORM_COMP_PL')
c_v_SI = Column(data=all_v_SI, name='V_ALPHA_PL')
c_major = Column(data=all_majors, name='MAJOR_DC')
c_minor = Column(data=all_majors, name='MINOR_DC')
c_pas = Column(data=all_pas, name='PA_DC')
##Make a da table
columns = [c_ids, c_names, c_ras, c_decs, c_comp_types, c_mod_types,
c_stokes_I_ref, c_stokes_SI, c_lin_mod_type, c_linpol_ref,
c_linpol_SI, c_rm, c_v_mod_type, c_v_ref, c_v_SI, c_major,
c_minor, c_pas]
table = Table(columns)
table.write('polarised_sky.fits', overwrite=True)
OK, now we have a sky model, let’s run the simulation (which takes 14 mins on my desktop). We’ll use the EDA2 array layout, and set the --IAU_order flag again. We’ll also set both X and Y beams to be a Gaussian beam with FWHM of 100 deg. This will let us see most of the sky, but chop the horizon off to make imaging nicer. This beam also has zero leakage, which means we keep the I/Q/U/V all nice and separate. We’ll phase things up to EoR1 so the galactic plane is down. Do it with float precision as we don’t really care too much about accuracy, just want some fun images
metafits="../metafits/1136380296_metafits_ppds.fits"
uvfits_name = "scary_sky"
array_layout = "../../test_installation/array_layouts/EDA2_layout_255.txt"
command = "run_woden.py "
command += " --ra0=60.0 --dec0=-27.0"
command += " --num_freq_channels=16 --num_time_steps=14"
command += " --freq_res=80e+3 --time_res=8.0"
command += " --cat_filename=polarised_sky.fits"
command += f" --metafits_filename={metafits}"
command += " --band_nums=1"
command += f" --output_uvfits_prepend={uvfits_name}"
command += " --primary_beam=Gaussian --gauss_beam_FWHM=100"
command += f" --array_layout={array_layout}"
command += " --sky_crop_components"
command += " --chunking_size=5e+10"
command += " --precision=float"
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-14 16:44:19 - INFO -
) ( )
( ( ( /( )\ ) ( /(
)\))( ')\())(()/( ( )\())
((_)()\ )((_)\ /(_)) )\ ((_)\
_(())\_)() ((_)(_))_ ((_) _((_)
\ \((_)/ // _ \ | \ | __|| \| |
\ \/\/ /| (_) || |) || _| | .` |
\_/\_/ \___/ |___/ |___||_|\_|
You are using wodenpy version/git hash: a43c383
2025-03-14 16:44:19 - INFO - Input arguments after parsing:
Phase centre: 60.00000, -27.00000 deg
Array central latitude: -26.703 deg
Array central longitude: 116.671 deg
Array central height: 377.827 m
Lowest channel frequency: 167.035 MHz
Channel frequency resolution: 80.000 kHz
Coarse band bandwidth: 1.280 MHz
Num channels per coarse band: 16
Start date: 2016-01-09T13:11:19
Time resolution: 8.0 (s)
Num time steps: 14
Have read 255 antenna positions from array layout file: ../../test_installation/array_layouts/EDA2_layout_255.txt
Will write outputs to: /home/jack-line/software/WODEN_dev/examples/polarisation
2025-03-14 16:44:20 - INFO - Obs epoch initial LST was 63.0347905707 deg
2025-03-14 16:44:20 - INFO - Setting initial J2000 LST to 62.8702463523 deg
2025-03-14 16:44:20 - INFO - Setting initial mjd to 57396.5495717593
2025-03-14 16:44:20 - INFO - After precession initial latitude of the array is -26.7414173539 deg
2025-03-14 16:44:20 - INFO - Precessing array layout to J2000
2025-03-14 16:44:20 - INFO - Using Gaussian primary beam with the following parameters:
Locked to pointing: HA -0.2 deg, Dec -26.74111165708009 deg
FWHM: 100.0 at reference frequency: 150.0 MHz
2025-03-14 16:44: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-03-14 16:44:23 - INFO - Sky model mapping took 0.0 mins
2025-03-14 16:44:23 - INFO - Have read in 788432 components
2025-03-14 16:44:23 - INFO - After cropping there are 394238 components
2025-03-14 16:44:23 - INFO - Will load libwoden from /home/jack-line/software/WODEN_dev/wodenpy/libwoden_float.so
2025-03-14 16:44:24 - INFO - Running in parallel on GPU.
Will read sky model using 8 threads
There are 3 sets of sky models to run
2025-03-14 16:44:24 - INFO - Reading set 0 sky models
2025-03-14 16:44:25 - INFO - From sky set 0 thread num 0 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:26 - INFO - From sky set 0 thread num 1 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:26 - INFO - From sky set 0 thread num 2 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:27 - INFO - From sky set 0 thread num 3 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:27 - INFO - From sky set 0 thread num 4 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:28 - INFO - From sky set 0 thread num 5 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:29 - INFO - From sky set 0 thread num 6 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:29 - INFO - From sky set 0 thread num 7 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:30 - INFO - Reading set 1 sky models
2025-03-14 16:44:31 - INFO - From sky set 1 thread num 0 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:31 - INFO - Sending Sky set 0 to the GPU
2025-03-14 16:44:31 - INFO - From sky set 1 thread num 1 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:32 - INFO - From sky set 1 thread num 2 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:33 - INFO - From sky set 1 thread num 3 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:33 - INFO - From sky set 1 thread num 4 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:34 - INFO - From sky set 1 thread num 5 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:35 - INFO - From sky set 1 thread num 6 reading 24608 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:44:35 - INFO - From sky set 1 thread num 7 reading 24590 points, 0 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:51:30 - INFO - Sky set 0 has returned from the GPU after 418.4 seconds
2025-03-14 16:51:30 - INFO - Have completed 196864 of 394238 components calcs (49.9%)
2025-03-14 16:51:30 - INFO - Reading set 2 sky models
2025-03-14 16:51:31 - INFO - From sky set 2 thread num 0 reading 0 points, 66 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:51:31 - INFO - Sending Sky set 1 to the GPU
2025-03-14 16:51:32 - INFO - From sky set 2 thread num 1 reading 0 points, 66 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:51:32 - INFO - From sky set 2 thread num 2 reading 0 points, 66 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:51:33 - INFO - From sky set 2 thread num 3 reading 0 points, 66 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:51:33 - INFO - From sky set 2 thread num 4 reading 0 points, 66 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:51:34 - INFO - From sky set 2 thread num 5 reading 0 points, 66 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:51:35 - INFO - From sky set 2 thread num 6 reading 0 points, 66 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:51:35 - INFO - From sky set 2 thread num 7 reading 0 points, 66 gauss, 0 shape, 0 shape coeffs
2025-03-14 16:58:19 - INFO - Sky set 1 has returned from the GPU after 407.2 seconds
2025-03-14 16:58:19 - INFO - Have completed 393710 of 394238 components calcs (99.9%)
2025-03-14 16:58:19 - INFO - Sending Sky set 2 to the GPU
2025-03-14 16:58:23 - INFO - Sky set 2 has returned from the GPU after 3.9 seconds
2025-03-14 16:58:23 - INFO - Have completed 394238 of 394238 components calcs (100.0%)
2025-03-14 16:58:23 - INFO - Finished all rounds of processing
2025-03-14 16:58:26 - INFO - Full run took 0:14:06.392452
2025-03-14 16:58:26 - INFO - WODEN is done. Closing the log. S'later
0
Now we can image it using WSClean. To be clear, normally you’d need to beam correct the images, and direction dependent differences between XX,XY,YX,YY mess with the conversion from instrumental to Stokes parameters. But as the primary beam is the same for both polarisations, you just get Stokes I/Q/U/V times a Gaussian beam after running through WSClean.
command = "woden_uv2ms.py "
command += f" --uvfits_prepend={uvfits_name}_band --band_nums=1"
call(command, shell=True)
command = f"wsclean -name {uvfits_name} -size 4096 4096 -niter 1000 "
command += " -auto-threshold 0.5 -auto-mask 3 "
command += " -pol I,Q,U,V -multiscale -weight briggs 0 -scale 0.03 -j 12 -mgain 0.85 "
command += " -no-update-model-required " #-join-polarizations
command += f" {uvfits_name}*.ms"
call(command, 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.5.0')
rm: cannot remove 'scary_sky_band01.ms': No such file or directory
The telescope frame is set to '????', which generally indicates ignorance. Defaulting the frame to 'itrs', but this may lead to other warnings or errors.
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 167-168 (16)
| Independent group:
+-+-J- 1 Q 0 1 1 1 0 0 167-168 (16)
| Independent group:
+-+-J- 2 U 0 2 2 2 0 0 167-168 (16)
| Independent group:
+-+-J- 3 V 0 3 3 3 0 0 167-168 (16)
Reordering scary_sky_band01.ms into 1 x 4 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 scary_sky_band01.ms
Detected 62.7 GB of system memory, usage not limited.
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Determining min and max w & theoretical beam size... DONE (w=[2.02049e-06:0.99262] lambdas, maxuvw=19.749 lambda)
Theoretic beam = 2.9 deg
Minimal inversion size: 102 x 102, using optimal: 108 x 108
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
Fitting beam... major=119.87', minor=119.53', PA=48.03 deg, theoretical=2.9 deg.
Writing psf image... DONE
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
Writing dirty image...
== Deconvolving (1) ==
Estimated standard deviation of background noise: 159.09 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 407.35 mJy, 387: 116.39 mJy, 774: 83.59 mJy, 1547: 57.43 mJy}
Starting multi-scale cleaning. Start peak=10.93 Jy, major iteration threshold=1.64 Jy
Iteration 4, scale 0 px : 8.28 Jy at 1952,2044
Iteration 13, scale 0 px : 6.49 Jy at 2003,1791
Iteration 28, scale 0 px : 5.08 Jy at 1793,1908
Iteration 52, scale 0 px : 4.05 Jy at 2417,2394
Iteration 94, scale 0 px : 3.24 Jy at 2374,2175
Iteration 147, scale 0 px : 2.59 Jy at 2508,2450
Iteration 220, scale 0 px : 2.06 Jy at 2278,1964
Iteration 304, scale 0 px : 1.64 Jy at 1517,1441
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 306, scale 0 px : 1.64 Jy at 2504,2446
Minor loop finished, continuing cleaning after inversion/prediction round.
Assigning from 1 to 1 channels...
== Converting model image to visibilities ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
== Deconvolving (2) ==
Estimated standard deviation of background noise: 74.81 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 171.59 mJy, 387: 45.45 mJy, 774: 26.22 mJy, 1547: 15.88 mJy}
Starting multi-scale cleaning. Start peak=1.69 Jy, major iteration threshold=253.76 mJy
Iteration 383, scale 0 px : 1.35 Jy at 2090,3070
Iteration 497, scale 0 px : 1.08 Jy at 2617,2665
Iteration 613, scale 0 px : 859.64 mJy at 1136,2509
Iteration 749, scale 0 px : 687.48 mJy at 2098,1258
Iteration 921, scale 0 px : 559.55 mJy at 2475,1168
Iteration 1000, scale 0 px : 499.7 mJy at 1541,1979
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 scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
2 major iterations were performed.
Rendering sources to restored image (beam=119.53'-119.87', PA=48.03 deg)... DONE
Writing restored image... DONE
Multi-scale cleaning summary:
- Scale 0 px, nr of components cleaned: 1000 (154.87 Jy)
- Scale 387 px, nr of components cleaned: 0 (0 Jy)
- Scale 774 px, nr of components cleaned: 0 (0 Jy)
- Scale 1547 px, nr of components cleaned: 0 (0 Jy)
Total: 1000 components (154.87 Jy)
Inversion: 00:00:08.256919, prediction: 00:00:02.163137, deconvolution: 00:00:20.704123
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Determining min and max w & theoretical beam size... DONE (w=[2.02049e-06:0.99262] lambdas, maxuvw=19.749 lambda)
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
Writing dirty image...
== Deconvolving (1) ==
Estimated standard deviation of background noise: 1.15 Jy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 1.2 Jy, 387: 1.17 Jy, 774: 1.09 Jy, 1547: 877.75 mJy}
Starting multi-scale cleaning. Start peak=-16.38 Jy, major iteration threshold=3.44 Jy (final)
Iteration 5, scale 1547 px : -12.62 Jy at 1906,2109
Iteration 10, scale 1547 px : -9.72 Jy at 1906,2109
Iteration 15, scale 1547 px : -7.49 Jy at 1826,2168
Iteration 23, scale 1547 px : -5.84 Jy at 1642,2307
Iteration 33, scale 1547 px : -4.59 Jy at 1597,2348
Iteration 48, scale 774 px : -3.87 Jy at 1710,1621
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 51, scale 1547 px : -3.64 Jy at 2419,2100
(6) Iteration 56, scale 1547 px : -3.38 Jy at 2024,1622
Cleaning finished because the final threshold was reached.
Auto-masking threshold reached; continuing next major iteration with deeper threshold and mask.
Assigning from 1 to 1 channels...
== Converting model image to visibilities ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
== Deconvolving (2) ==
Estimated standard deviation of background noise: 578.96 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 553.64 mJy, 387: 522.39 mJy, 774: 479.75 mJy, 1547: 363.44 mJy}
Starting multi-scale cleaning. Start peak=-3.45 Jy, major iteration threshold=516.81 mJy
Iteration 73, scale 774 px : -2.89 Jy at 1710,1621
Iteration 77, scale 1547 px : -2.78 Jy at 2419,2100
Iteration 91, scale 1547 px : -2.21 Jy at 2419,2100
Iteration 109, scale 774 px : -1.95 Jy at 1710,1621
Iteration 113, scale 1547 px : -1.8 Jy at 2419,2100
Iteration 127, scale 1547 px : 1.43 Jy at 1491,819
Iteration 145, scale 774 px : -1.36 Jy at 1710,1621
Iteration 149, scale 774 px : -1.08 Jy at 1710,1621
Iteration 153, scale 1547 px : -1.17 Jy at 2419,2100
Iteration 165, scale 1547 px : -930.79 mJy at 2419,2100
Iteration 183, scale 1547 px : 890.55 mJy at 1906,2109
Iteration 201, scale 1547 px : -840.1 mJy at 1888,1660
Iteration 231, scale 774 px : -694.65 mJy at 1710,1621
Iteration 235, scale 1547 px : -670.28 mJy at 1559,2399
Iteration 299, scale 1547 px : -531.11 mJy at 2419,2100
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 314, scale 1547 px : -512.48 mJy at 1564,2400
Minor loop finished, continuing cleaning after inversion/prediction round.
Assigning from 1 to 1 channels...
== Converting model image to visibilities ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
== Deconvolving (3) ==
Estimated standard deviation of background noise: 330.12 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 372.97 mJy, 387: 335.63 mJy, 774: 286.88 mJy, 1547: 194.29 mJy}
Starting multi-scale cleaning. Start peak=-612.04 mJy, major iteration threshold=165.06 mJy (final)
Iteration 329, scale 1547 px : -544.22 mJy at 2419,2100
Iteration 390, scale 1547 px : 434.82 mJy at 1491,819
Iteration 479, scale 1547 px : -347.51 mJy at 1888,1660
Iteration 570, scale 1547 px : 275.86 mJy at 2085,1985
Iteration 667, scale 1547 px : -220.34 mJy at 1888,1660
Iteration 757, scale 1547 px : -176.06 mJy at 1564,2400
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 787, scale 1547 px : 163.64 mJy at 2085,1985
Cleaning finished because the final threshold was reached.
Assigning from 1 to 1 channels...
Writing model image...
== Converting model image to visibilities ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
3 major iterations were performed.
Rendering sources to restored image (beam=119.53'-119.87', PA=48.03 deg)... DONE
Writing restored image... DONE
Multi-scale cleaning summary:
- Scale 0 px, nr of components cleaned: 0 (0 Jy)
- Scale 387 px, nr of components cleaned: 0 (0 Jy)
- Scale 774 px, nr of components cleaned: 23 (-53.1 Jy)
- Scale 1547 px, nr of components cleaned: 764 (-1.95 KJy)
Total: 787 components (-2 KJy)
Inversion: 00:00:13.748209, prediction: 00:00:05.396248, deconvolution: 00:02:15.731349
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Determining min and max w & theoretical beam size... DONE (w=[2.02049e-06:0.99262] lambdas, maxuvw=19.749 lambda)
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
Writing dirty image...
== Deconvolving (1) ==
Estimated standard deviation of background noise: 265.47 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 277.57 mJy, 387: 269.08 mJy, 774: 254.38 mJy, 1547: 202.53 mJy}
Starting multi-scale cleaning. Start peak=3.8 Jy, major iteration threshold=796.42 mJy (final)
Iteration 5, scale 1547 px : 2.92 Jy at 1906,2109
Iteration 10, scale 1547 px : 2.25 Jy at 1906,2109
Iteration 15, scale 1547 px : 1.74 Jy at 1988,2051
Iteration 23, scale 1547 px : 1.36 Jy at 2215,1915
Iteration 34, scale 1547 px : 1.04 Jy at 1592,2347
Iteration 50, scale 774 px : 877.01 mJy at 1706,1622
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 52, scale 1547 px : 828.27 mJy at 1564,2403
(6) Iteration 55, scale 1547 px : 789.35 mJy at 1563,2408
Cleaning finished because the final threshold was reached.
Auto-masking threshold reached; continuing next major iteration with deeper threshold and mask.
Assigning from 1 to 1 channels...
== Converting model image to visibilities ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
== Deconvolving (2) ==
Estimated standard deviation of background noise: 134.81 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 129.06 mJy, 387: 122.88 mJy, 774: 111.48 mJy, 1547: 85.22 mJy}
Starting multi-scale cleaning. Start peak=820.55 mJy, major iteration threshold=123.08 mJy
Iteration 71, scale 774 px : 710.69 mJy at 1706,1622
Iteration 75, scale 1547 px : 654.88 mJy at 2427,2103
Iteration 92, scale 1547 px : 505.15 mJy at 1869,1656
Iteration 111, scale 774 px : 458.2 mJy at 1706,1622
Iteration 115, scale 774 px : 364.08 mJy at 1706,1622
Iteration 119, scale 1547 px : 413.06 mJy at 2427,2103
Iteration 131, scale 1547 px : 328.85 mJy at 2427,2103
Iteration 150, scale 1547 px : 257.54 mJy at 1564,2403
Iteration 168, scale 774 px : 231.16 mJy at 1706,1622
Iteration 172, scale 774 px : 183.68 mJy at 1706,1622
Iteration 176, scale 1547 px : 209.58 mJy at 1564,2403
Iteration 206, scale 1547 px : 203.63 mJy at 1869,1656
Iteration 212, scale 1547 px : -186.38 mJy at 1988,2051
Iteration 252, scale 1547 px : 148.68 mJy at 1564,2403
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 313, scale 1547 px : 122.57 mJy at 1869,1656
Minor loop finished, continuing cleaning after inversion/prediction round.
Assigning from 1 to 1 channels...
== Converting model image to visibilities ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
== Deconvolving (3) ==
Estimated standard deviation of background noise: 76 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 86.21 mJy, 387: 77.01 mJy, 774: 66.52 mJy, 1547: 45.09 mJy}
Starting multi-scale cleaning. Start peak=143.59 mJy, major iteration threshold=38 mJy (final)
Iteration 330, scale 1547 px : 127.86 mJy at 2423,2109
Iteration 376, scale 1547 px : 102.04 mJy at 1869,1656
Iteration 447, scale 1547 px : -80.48 mJy at 1988,2051
Iteration 517, scale 1547 px : 64.15 mJy at 2423,2109
Iteration 591, scale 1547 px : 51 mJy at 1869,1656
Iteration 687, scale 1547 px : 40.71 mJy at 2423,2109
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 714, scale 1547 px : 37.99 mJy at 1564,2403
Cleaning finished because the final threshold was reached.
Assigning from 1 to 1 channels...
Writing model image...
== Converting model image to visibilities ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
3 major iterations were performed.
Rendering sources to restored image (beam=119.53'-119.87', PA=48.03 deg)... DONE
Writing restored image... DONE
Multi-scale cleaning summary:
- Scale 0 px, nr of components cleaned: 0 (0 Jy)
- Scale 387 px, nr of components cleaned: 0 (0 Jy)
- Scale 774 px, nr of components cleaned: 22 (11.67 Jy)
- Scale 1547 px, nr of components cleaned: 692 (453.11 Jy)
Total: 714 components (464.79 Jy)
Inversion: 00:00:19.245613, prediction: 00:00:08.507912, deconvolution: 00:04:07.420747
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Determining min and max w & theoretical beam size... DONE (w=[2.02049e-06:0.99262] lambdas, maxuvw=19.749 lambda)
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
Writing dirty image...
== Deconvolving (1) ==
Estimated standard deviation of background noise: 280.57 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 333.71 mJy, 387: 317.21 mJy, 774: 285.52 mJy, 1547: 219.77 mJy}
Starting multi-scale cleaning. Start peak=4.17 Jy, major iteration threshold=841.71 mJy (final)
Iteration 5, scale 774 px : 3.29 Jy at 2044,2166
Iteration 9, scale 1547 px : 3.01 Jy at 1949,1779
Iteration 14, scale 1547 px : 2.32 Jy at 1950,1780
Iteration 20, scale 774 px : 2.15 Jy at 2094,2222
Iteration 26, scale 774 px : 1.72 Jy at 2094,2222
Iteration 48, scale 1547 px : 1.39 Jy at 1719,1652
Iteration 57, scale 774 px : 1.32 Jy at 2225,1423
Iteration 77, scale 774 px : 1.05 Jy at 2490,1787
Iteration 107, scale 1547 px : 848.52 mJy at 2456,2048
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 108, scale 774 px : 843.91 mJy at 2247,1420
(6) Iteration 109, scale 1547 px : 842.97 mJy at 1568,1706
(5) Iteration 110, scale 387 px : -845.66 mJy at 2043,1774
(4) Iteration 111, scale 774 px : 822.16 mJy at 2200,2330
Cleaning finished because the final threshold was reached.
Auto-masking threshold reached; continuing next major iteration with deeper threshold and mask.
Assigning from 1 to 1 channels...
== Converting model image to visibilities ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
== Deconvolving (2) ==
Estimated standard deviation of background noise: 110.52 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 134.8 mJy, 387: 124.41 mJy, 774: 104.45 mJy, 1547: 73.96 mJy}
Starting multi-scale cleaning. Start peak=831.2 mJy, major iteration threshold=124.68 mJy
Iteration 122, scale 387 px : -850.99 mJy at 2043,1774
Iteration 126, scale 774 px : 823.13 mJy at 2243,1419
Iteration 145, scale 774 px : 672.34 mJy at 1442,1680
Iteration 178, scale 387 px : -614.95 mJy at 2043,1774
Iteration 182, scale 387 px : -486.27 mJy at 2043,1774
Iteration 186, scale 1547 px : 531.63 mJy at 2456,2048
Iteration 198, scale 774 px : 547.1 mJy at 2247,1420
Iteration 211, scale 774 px : 437.27 mJy at 2683,1951
Iteration 244, scale 1547 px : 357.51 mJy at 1568,1706
Iteration 253, scale 387 px : -407.91 mJy at 2043,1774
Iteration 257, scale 387 px : -322.55 mJy at 2043,1774
Iteration 261, scale 774 px : 353.73 mJy at 2247,1420
Iteration 281, scale 774 px : 293.83 mJy at 1442,1680
Iteration 311, scale 1547 px : 247.95 mJy at 1568,1706
Iteration 319, scale 387 px : -244.22 mJy at 2043,1774
Iteration 323, scale 774 px : 242.66 mJy at 2247,1420
Iteration 342, scale 774 px : 196.03 mJy at 1442,1680
Iteration 375, scale 387 px : -177.25 mJy at 2043,1774
Iteration 379, scale 1547 px : 170.2 mJy at 1568,1706
Iteration 386, scale 774 px : 159.27 mJy at 2247,1420
Iteration 409, scale 774 px : 124.97 mJy at 1705,2210
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 410, scale 387 px : -141.87 mJy at 2043,1774
(6) Iteration 413, scale 1547 px : 135.53 mJy at 1568,1706
(5) Iteration 415, scale 774 px : 124.91 mJy at 2195,2310
(4) Iteration 416, scale 774 px : 124.96 mJy at 2683,1951
(3) Iteration 417, scale 774 px : 123.64 mJy at 2247,1420
Minor loop finished, continuing cleaning after inversion/prediction round.
Assigning from 1 to 1 channels...
== Converting model image to visibilities ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
== Deconvolving (3) ==
Estimated standard deviation of background noise: 48.2 mJy
Scale info:
- Scale 0, bias factor=1, psfpeak=1, gain=0.1, kernel peak=0.000118546
- Scale 387, bias factor=1.7, psfpeak=0.109965, gain=0.909382, kernel peak=2.99227e-05
- Scale 774, bias factor=2.8, psfpeak=0.0273488, gain=3.65646, kernel peak=7.55434e-06
- Scale 1547, bias factor=4.6, psfpeak=0.00622096, gain=16.0747, kernel peak=1.89546e-06
RMS per scale: {0: 76.43 mJy, 387: 64.64 mJy, 774: 47.82 mJy, 1547: 29.25 mJy}
Starting multi-scale cleaning. Start peak=142.49 mJy, major iteration threshold=24.1 mJy (final)
Iteration 422, scale 774 px : 122.98 mJy at 1705,2210
Iteration 447, scale 387 px : -120.41 mJy at 2043,1774
Iteration 451, scale 1547 px : 103.99 mJy at 1568,1706
Iteration 458, scale 387 px : -99.49 mJy at 2043,1774
Iteration 462, scale 774 px : 99.12 mJy at 2247,1420
Iteration 484, scale 774 px : 79.26 mJy at 2044,2165
Iteration 519, scale 1547 px : -80 mJy at 1950,1780
Iteration 534, scale 1547 px : 65.44 mJy at 2456,2048
Iteration 567, scale 774 px : 86.34 mJy at 2247,1420
Iteration 575, scale 774 px : 69.06 mJy at 2247,1420
Iteration 594, scale 774 px : 55.03 mJy at 2449,1746
Iteration 616, scale 1547 px : -77.41 mJy at 1950,1780
Iteration 621, scale 1547 px : -59.63 mJy at 1950,1780
Iteration 641, scale 1547 px : 53.4 mJy at 2456,2048
Iteration 661, scale 1547 px : 42.56 mJy at 1568,1706
Iteration 695, scale 774 px : 72.87 mJy at 2103,1409
Iteration 704, scale 774 px : 57.55 mJy at 2044,2165
Iteration 734, scale 774 px : -45.84 mJy at 2511,1817
Iteration 822, scale 774 px : 48.56 mJy at 2683,1951
Iteration 826, scale 774 px : -41.81 mJy at 2511,1817
Iteration 849, scale 774 px : 44.18 mJy at 2449,1746
Iteration 853, scale 774 px : -40.94 mJy at 2511,1817
Iteration 910, scale 774 px : 45.83 mJy at 2683,1951
Iteration 914, scale 774 px : -37.54 mJy at 2511,1817
Iteration 936, scale 1547 px : -57.37 mJy at 1950,1780
Iteration 941, scale 1547 px : 51.87 mJy at 1568,1706
Iteration 955, scale 1547 px : -40.56 mJy at 1950,1780
Iteration 989, scale 774 px : -52.23 mJy at 1545,1616
Iteration 1000, scale 774 px : 42 mJy at 2103,1409
Cleaning finished because maximum number of iterations was reached.
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 scary_sky_band01.ms
Predicting 453390 rows...
Writing...
== Constructing image ==
Opening reordered part 0 spw 0 for scary_sky_band01.ms
Loading data in memory...
Gridding 453390 rows...
Gridded visibility count: 7.25424e+06
3 major iterations were performed.
Rendering sources to restored image (beam=119.53'-119.87', PA=48.03 deg)... DONE
Writing restored image... DONE
Multi-scale cleaning summary:
- Scale 0 px, nr of components cleaned: 0 (0 Jy)
- Scale 387 px, nr of components cleaned: 40 (-7.34 Jy)
- Scale 774 px, nr of components cleaned: 692 (266.19 Jy)
- Scale 1547 px, nr of components cleaned: 268 (277.32 Jy)
Total: 1000 components (536.17 Jy)
Inversion: 00:00:24.659906, prediction: 00:00:11.692319, deconvolution: 00:07:54.554818
Cleaning up temporary files...
0
Now we can plot the images
from astropy.wcs import WCS
with fits.open(f'{uvfits_name}-I-image.fits') as hdus:
w = WCS(hdus[0].header)
fig, axs = plt.subplots(2, 2, figsize=(10, 8), subplot_kw={'projection': w.celestial})
names = ['I', 'Q', 'U', 'V']
files = [f'{uvfits_name}-I-image.fits', f'{uvfits_name}-Q-image.fits',
f'{uvfits_name}-U-image.fits', f'{uvfits_name}-V-image.fits']
for name, file, ax in zip(names, files, axs.flatten()):
with fits.open(file) as hdus:
data = np.squeeze(hdus[0].data)
im = ax.imshow(data, origin='lower', cmap='inferno')
ax.set_title(name)
ax.set_xticks([])
ax.set_yticks([])
plt.colorbar(im, ax=ax, label='Jy/beam')
ax.grid(alpha=0.5)
plt.tight_layout()
plt.show()
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 57396.549571 from DATE-OBS'. [astropy.wcs.wcs]
Great success! We have boring point sources in Stokes I. As we have a constant rotation measure for all directions, and we’ve only done a 1.28MHz bandwidth, at these frequencies we get a negative Stokes Q and positive Stokes U; both look like the (low-res) Haslam map. The Stokes V is just a bunch of Gaussian blobs.