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, and WSClean to do the imaging.
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'
##use subprocess to run the command
call(command, shell=True)
LST: 100.31812162623258, RA: 98.99999999999999
/home/jack-line/software/anaconda3/envs/woden/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.3.0')
You are using WODEN commit: No git describe nor __version__ avaible
LOADING IN /home/jack-line/software/WODEN/wodenpy/libwoden_double.so
Setting phase centre RA,DEC 99.00000deg -20.60000deg
Obs epoch initial LST was 100.3222997096 deg
Setting initial J2000 LST to 100.0765133194 deg
Setting initial mjd to 60496.1666782406
After precession initial latitude of the array is -26.6814845743 deg
We are precessing the array
Doing the initial reading/mapping of sky model into chunks
Mapping took 0.0 mins
Have read in 1 components
After cropping there are 1 components
After chunking there are 1 chunks
Reading chunks skymodel chunks 0:50
YO ['PRIMARY', 'MAIN', 'V_LIST_FLUXES']
Reading chunks 0:50 took 0.0 mins
WODEN is using DOUBLE precision
NO PRIMARY BEAM HAS BEEN SELECTED
Will run without a primary beam
Simulating band 01 with bottom freq 1.80000000e+08
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 0
Number of components in chunk are: P 1 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
GPU calls for band 1 finished
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'
##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/anaconda3/envs/woden/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.3.0')
You are using WODEN commit: No git describe nor __version__ avaible
LOADING IN /home/jack-line/software/WODEN/wodenpy/libwoden_double.so
Setting phase centre RA,DEC 99.00000deg -20.60000deg
Obs epoch initial LST was 100.3222997096 deg
Setting initial J2000 LST to 100.0765133194 deg
Setting initial mjd to 60496.1666782406
After precession initial latitude of the array is -26.6814845743 deg
We are precessing the array
Doing the initial reading/mapping of sky model into chunks
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
Mapping took 0.0 mins
Have read in 3 components
After cropping there are 3 components
After chunking there are 1 chunks
Reading chunks skymodel chunks 0:50
YO ['PRIMARY', '']
Reading chunks 0:50 took 0.0 mins
WODEN is using DOUBLE precision
NO PRIMARY BEAM HAS BEEN SELECTED
Will run without a primary beam
Simulating band 01 with bottom freq 1.80000000e+08
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 0
Number of components in chunk are: P 3 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
GPU calls for band 1 finished
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. 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/anaconda3/envs/woden/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.3.0')
You are using WODEN commit: No git describe nor __version__ avaible
LOADING IN /home/jack-line/software/WODEN/wodenpy/libwoden_float.so
Setting phase centre RA,DEC 60.00000deg -27.00000deg
Obs epoch initial LST was 63.0347905707 deg
Setting initial J2000 LST to 62.8702463523 deg
Setting initial mjd to 57396.5495717593
After precession initial latitude of the array is -26.7414173539 deg
We are precessing the array
Doing the initial reading/mapping of sky model into chunks
INFO: couldn't find second table containing shapelet information, so not attempting to load any shapelets.
Mapping took 0.0 mins
Have read in 788432 components
After cropping there are 394238 components
After chunking there are 59 chunks
Reading chunks skymodel chunks 0:50
YO ['PRIMARY', '']
Reading chunks 0:50 took 0.1 mins
Reading chunks skymodel chunks 50:100
WODEN is using FLOAT precision
Setting up Gaussian primary beam settings
pointing at HA, Dec = -0.21717deg, -26.74111deg
setting beam FWHM to 100.00000deg and ref freq to 150.000MHz
Simulating band 01 with bottom freq 1.67035000e+08
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 0
Number of components in chunk are: P 6892 G 0 S_coeffs 0
YO ['PRIMARY', '']
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Reading chunks 50:100 took 0.0 mins
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 1
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 2
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 3
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 4
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 5
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 6
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 7
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 8
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 9
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 10
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 11
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 12
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 13
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 14
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 15
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 16
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 17
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 18
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 19
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 20
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 21
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 22
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 23
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 24
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 25
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 26
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 27
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 28
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 29
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 30
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 31
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 32
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 33
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 34
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 35
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 36
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 37
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 38
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 39
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 40
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 41
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 42
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 43
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 44
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 45
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 46
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 47
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 48
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 49
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
GPU calls for band 1 finished
WODEN is using FLOAT precision
Setting up Gaussian primary beam settings
pointing at HA, Dec = -0.21717deg, -26.74111deg
setting beam FWHM to 100.00000deg and ref freq to 150.000MHz
Simulating band 01 with bottom freq 1.67035000e+08
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 0
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 1
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 2
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 3
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 4
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 5
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 6
Number of components in chunk are: P 6892 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 7
Number of components in chunk are: P 866 G 0 S_coeffs 0
Doing point components
Extrapolating fluxes and beams...
Doing Gaussian Beam
Extrapolating fluxes and beams done.
Doing visi kernel...
Visi kernel done
About to copy the chunked source to the GPU
Have copied across the chunk to the GPU
Processing chunk 8
Number of components in chunk are: P 0 G 528 S_coeffs 0
Doing gaussian components
Doing Gaussian Beam
GPU calls for band 1 finished
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/anaconda3/envs/woden/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.3.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.
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:07.834198, prediction: 00:00:02.039755, deconvolution: 00:00:19.873050
== 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 1905,2109
Iteration 10, scale 1547 px : -9.72 Jy at 1906,2109
Iteration 15, scale 1547 px : -7.49 Jy at 1827,2167
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,2099
(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.83 mJy
Iteration 73, scale 774 px : -2.89 Jy at 1710,1621
Iteration 77, scale 1547 px : -2.78 Jy at 2419,2099
Iteration 91, scale 1547 px : -2.21 Jy at 1559,2399
Iteration 109, scale 774 px : -1.95 Jy at 1710,1621
Iteration 113, scale 1547 px : -1.79 Jy at 2419,2099
Iteration 127, scale 1547 px : 1.43 Jy at 1491,819
Iteration 145, scale 774 px : -1.37 Jy at 1710,1621
Iteration 149, scale 774 px : -1.09 Jy at 1710,1621
Iteration 153, scale 1547 px : -1.17 Jy at 2419,2099
Iteration 165, scale 1547 px : -930.58 mJy at 2423,2090
Iteration 183, scale 1547 px : 881.46 mJy at 1906,2109
Iteration 215, scale 774 px : -736.58 mJy at 1710,1621
Iteration 219, scale 1547 px : -714.56 mJy at 1559,2399
Iteration 274, scale 1547 px : 571.09 mJy at 2085,1985
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 317, scale 1547 px : -511.91 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: 329.8 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.9 mJy, 387: 335.53 mJy, 774: 286.8 mJy, 1547: 194.24 mJy}
Starting multi-scale cleaning. Start peak=-611.6 mJy, major iteration threshold=164.9 mJy (final)
Iteration 332, scale 1547 px : -537.84 mJy at 2419,2099
Iteration 394, scale 1547 px : -426.79 mJy at 2414,2103
Iteration 488, scale 1547 px : 340.63 mJy at 2085,1985
Iteration 581, scale 1547 px : -270.45 mJy at 2414,2103
Iteration 675, scale 1547 px : -214.08 mJy at 1564,2400
Iteration 770, scale 1547 px : 171.16 mJy at 2085,1985
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 782, scale 1547 px : -173.58 mJy at 1564,2400
(6) Iteration 783, scale 1547 px : -164.76 mJy at 1564,2400
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.32 Jy)
- Scale 1547 px, nr of components cleaned: 760 (-1.95 KJy)
Total: 783 components (-2 KJy)
Inversion: 00:00:13.401442, prediction: 00:00:05.080051, deconvolution: 00:02:21.242978
== 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 1905,2109
Iteration 10, scale 1547 px : 2.25 Jy at 1905,2109
Iteration 15, scale 1547 px : 1.74 Jy at 1822,2171
Iteration 23, scale 1547 px : 1.35 Jy at 1641,2308
Iteration 33, scale 1547 px : 1.07 Jy at 1596,2348
Iteration 48, scale 774 px : 899.22 mJy at 1710,1621
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 51, scale 1547 px : 846.82 mJy at 2419,2100
(6) Iteration 56, scale 1547 px : 786.3 mJy at 2025,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: 134.29 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: 128.49 mJy, 387: 122.32 mJy, 774: 110.95 mJy, 1547: 84.84 mJy}
Starting multi-scale cleaning. Start peak=799.79 mJy, major iteration threshold=119.97 mJy
Iteration 73, scale 774 px : 671.14 mJy at 1710,1621
Iteration 77, scale 1547 px : 646.6 mJy at 2419,2100
Iteration 91, scale 1547 px : 513.85 mJy at 2419,2100
Iteration 109, scale 774 px : 451.26 mJy at 1710,1621
Iteration 113, scale 1547 px : 417.72 mJy at 2419,2100
Iteration 127, scale 1547 px : -333.91 mJy at 1493,817
Iteration 145, scale 774 px : 315.94 mJy at 1710,1621
Iteration 149, scale 774 px : 251.04 mJy at 1710,1621
Iteration 153, scale 1547 px : 272.2 mJy at 2419,2100
Iteration 165, scale 1547 px : 216.63 mJy at 2423,2091
Iteration 183, scale 1547 px : -204.21 mJy at 1906,2109
Iteration 215, scale 774 px : 169.96 mJy at 1710,1621
Iteration 219, scale 1547 px : 166.3 mJy at 1559,2400
Iteration 273, scale 1547 px : 132.86 mJy at 2419,2100
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 312, scale 1547 px : 119.96 mJy at 1887,1660
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.9 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: 87.1 mJy, 387: 77.89 mJy, 774: 67.12 mJy, 1547: 45.45 mJy}
Starting multi-scale cleaning. Start peak=140.88 mJy, major iteration threshold=38.45 mJy (final)
Iteration 329, scale 1547 px : 127.16 mJy at 2419,2100
Iteration 383, scale 1547 px : 101.71 mJy at 2414,2104
Iteration 479, scale 1547 px : 80.35 mJy at 1887,1660
Iteration 575, scale 1547 px : 63.77 mJy at 1887,1660
Iteration 669, scale 1547 px : 50.82 mJy at 1563,2400
Iteration 763, scale 1547 px : -40.08 mJy at 2088,1984
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 788, scale 1547 px : 38.3 mJy at 1563,2400
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 (12.36 Jy)
- Scale 1547 px, nr of components cleaned: 765 (451.69 Jy)
Total: 788 components (464.04 Jy)
Inversion: 00:00:18.859954, prediction: 00:00:08.086930, deconvolution: 00:04:21.338979
== 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.7 mJy, 387: 317.21 mJy, 774: 285.52 mJy, 1547: 219.76 mJy}
Starting multi-scale cleaning. Start peak=4.17 Jy, major iteration threshold=841.7 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 1949,1779
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.53 mJy at 2456,2048
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 108, scale 774 px : 843.9 mJy at 2247,1420
(6) Iteration 109, scale 1547 px : 842.98 mJy at 1568,1706
(5) Iteration 110, scale 387 px : -845.64 mJy at 2043,1774
(4) Iteration 111, scale 774 px : 822.13 mJy at 2200,2329
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.29 mJy, major iteration threshold=124.69 mJy
Iteration 122, scale 387 px : -850.94 mJy at 2043,1774
Iteration 126, scale 774 px : 823.11 mJy at 2243,1419
Iteration 145, scale 774 px : 672.37 mJy at 1442,1680
Iteration 178, scale 387 px : -614.92 mJy at 2043,1774
Iteration 182, scale 387 px : -486.24 mJy at 2043,1774
Iteration 186, scale 1547 px : 531.57 mJy at 2456,2048
Iteration 198, scale 774 px : 547.1 mJy at 2247,1420
Iteration 211, scale 774 px : 437.29 mJy at 2683,1951
Iteration 244, scale 1547 px : 357.69 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.75 mJy at 1442,1680
Iteration 311, scale 1547 px : 248.13 mJy at 1568,1706
Iteration 319, scale 387 px : -244.22 mJy at 2043,1774
Iteration 323, scale 774 px : 242.67 mJy at 2247,1420
Iteration 342, scale 774 px : 195.92 mJy at 1442,1680
Iteration 375, scale 387 px : -177.24 mJy at 2043,1774
Iteration 379, scale 1547 px : 170.35 mJy at 1568,1706
Iteration 386, scale 774 px : 159.19 mJy at 2247,1420
Iteration 409, scale 774 px : 124.76 mJy at 1706,2210
Subminor loop is near minor loop threshold. Initiating countdown.
(7) Iteration 410, scale 387 px : -141.9 mJy at 2043,1774
(6) Iteration 413, scale 1547 px : 136.01 mJy at 1568,1706
(5) Iteration 415, scale 774 px : 124.75 mJy at 2195,2310
(4) Iteration 416, scale 774 px : 125.02 mJy at 2683,1951
(3) Iteration 417, scale 1547 px : 123.5 mJy at 1568,1706
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.19 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.45 mJy, 387: 64.65 mJy, 774: 47.84 mJy, 1547: 29.26 mJy}
Starting multi-scale cleaning. Start peak=142.9 mJy, major iteration threshold=24.1 mJy (final)
Iteration 422, scale 774 px : 122.74 mJy at 1706,2210
Iteration 447, scale 387 px : -120.44 mJy at 2043,1774
Iteration 451, scale 1547 px : 104.05 mJy at 1568,1706
Iteration 458, scale 387 px : -99.51 mJy at 2043,1774
Iteration 462, scale 774 px : 99.08 mJy at 2247,1420
Iteration 484, scale 774 px : 79.24 mJy at 2683,1951
Iteration 519, scale 1547 px : -79.34 mJy at 1950,1780
Iteration 536, scale 1547 px : 66.36 mJy at 2456,2048
Iteration 564, scale 774 px : 85.36 mJy at 2247,1420
Iteration 574, scale 774 px : 65.92 mJy at 1706,2210
Iteration 594, scale 774 px : 52.71 mJy at 1868,1383
Iteration 619, scale 1547 px : -79.29 mJy at 1950,1780
Iteration 624, scale 1547 px : -61.07 mJy at 1950,1780
Iteration 643, scale 1547 px : 52.56 mJy at 2456,2048
Iteration 669, scale 1547 px : 41.52 mJy at 1568,1706
Iteration 704, scale 774 px : 72.82 mJy at 2103,1409
Iteration 713, scale 774 px : 57.92 mJy at 2044,2165
Iteration 745, scale 774 px : 46.28 mJy at 2044,2165
Iteration 832, scale 774 px : 45.75 mJy at 2683,1951
Iteration 842, scale 774 px : 39.94 mJy at 2449,1746
Iteration 986, scale 774 px : -31.64 mJy at 2511,1817
Iteration 1000, scale 1547 px : -56.93 mJy at 1950,1780
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: 742 (266.34 Jy)
- Scale 1547 px, nr of components cleaned: 218 (276.64 Jy)
Total: 1000 components (535.64 Jy)
Inversion: 00:00:24.572327, prediction: 00:00:11.293068, deconvolution: 00:08:03.705232
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.549572 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 directinos, 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.