SIMIO

Create a Template

Step 1: Get the ms file of the observation you want to mimic.

For this tutorial, we will create a template based on Elias 24, one of the bright DSHARP planet-forming disks (Huang et al. 2018). The original ms file can be obtained in the DSHARP database.

Download the continuum of Elias 24 (0.8GB). This file contains all the infortmation of the observation, including antenna, visibilities, times, etc.

Step 2: Create a folder for your template.

Go to you SIMIO directory, and open your templates folder. Create an empty folder with the name of your template. In this case, the template will be called Elias24.

Step 3: Create sub folders and scripts.

Create 3 empty folders, named "images", "msfiles" and "uvtables". Additionally, create 2 empty scripts named "template_casa_handling.py" and "template_info.py". In this example the template is Elias24, therefore the templates will be called: "Elias24_casa_handling.py" and "Elias24_info.py".

Move the ms file downloaded in step 1 into the folder "msfiles", and extract the ms file from the .tgz.

# Source Geometry, Elias 2-24 in Huang et al. (2018), DSHARP II

dRa = 0. # arcsec

dDec = 0. # arcsec

inc = 29. # degrees

pa = 45.7 # degrees

geom_params = [dRa, dDec, inc, pa]

# Distance

parallax = 7.180888047224439 # GAIA EDR3 mas

dist = 1000. / parallax # pc

# Sky info

temp_center_ra = '16h26m24.0782s' # J2000

temp_center_dec = '-24d16m13.884s' # J2000

temp_semimaj = 1.15 # arcsec

# Clean info

temp_robust = 0.0 # From Andrews et al. (2021)

temp_cellsize = '0.003arcsec' # DSHARP Standard

temp_imsize = 2000 # Big enough for noise estimation

Step 4: Template information.

The file "template_info.py", in this case "Elias24_info.py", contains important information to generate the images for this template. You also can copy/paste the info file of another template and replace the numbers with the ones of your template.

Source Geometry: Obtained from Huang et al. (DSHARP II, 2018). You can add a spatial offset if needed, and the inclination and position angle must be input in degrees. The ms file of Elias 27 has a phase center different from the center of the source, but we will change that in the next step. Therefore, we set dRa and dDec to 0.

Distance: We estimate the distance as the inverse of the parallax. We use GAIA EDR3 (2016, 2021) to obtain the parallax (SIMBAD name of the source is Elias 2-24).

Sky info: Here you should write the coordinates of the source center, in J2000 epoch. The coordinate we chose for Elias 24 is peak brightness pixel. The parameter temp_semimaj is the radial extent of the source in arcsec. This parameter is used to create the mask to generate its images.

Clean info: This parameters are needed to generate the template images. The robust parameters controls the weight that is given to the baselines. In this case, we know the original study was made on images generated with a robust parameter of 0.0 (Andrews et al. DSHARP I, 2021). The cellsize should be around 7x to 10x smaller than the beam size, and the image size should be big enough such that we can estimate the noise level with an annulus around the source.

# Append the au functions from other directory.

sys.path.append('../../codes/analysis_scripts/')

import analysisUtils as au

# Import uv-handling functions

execfile('../../codes/simio_ms2ascii.py')

# Import imaging functions

execfile('../../codes/simio_clean.py')

# Info

prefix = 'Elias24'

ms_cont = 'msfiles/'+prefix+'_avg.ms'

listobs('msfiles/Elias24_continuum.ms')

Step 5.1: Open CASA to prepare ms file.

Generally speaking, high angular resolution observations of disks can result in very heavy ms files. In SIMIO, we are mainly interested in checking the observability of your model under the uv-coverage of a template observation, therefore we can reduce the ms files volume by averaging in frequency and time, without considerably modifying the uv-coverage. The goal of this step is to reduce the data volume of the original observation, thus speeding up all the upcoming processes.

Open CASA in your terminal, from the folder /home/path_to_simio/simio/templates/Elias24/ . From this location, we will call the CASA analysis utils, and also the SIMIO scripts to extract uvtables and generate images.

Set the name of the template as prefix, and create the name for the template measurement set. Execute the task listobs on the downloaded ms file.

Step 5.2: Check the ms file

Here I show a small part of the information given by listobs. Here you can see that the visibilities of Elias 24 are stored in the ms file in 16 spectral windows (numbered from 0 to 15). We will average each spectral window to only 1 channel, and we will also average in time to 24s. This is an aggressive averaging, but it is ok for our comparison goal.

# Remove previous auxiliary files

os.system('rm -rf msfiles/Elias24_avg.ms')

os.system('rm -rf msfiles/Elias24_avg_aux1.ms')

os.system('rm -rf msfiles/Elias24_avg_aux2.ms')

os.system('rm -rf msfiles/Elias24_avg_aux3.ms')

# Split the data to average and combine

split(vis='msfiles/Elias24_continuum.ms', \

outputvis='msfiles/Elias24_avg_aux1.ms', \

timebin='30s', \

spw='0,1,3,4,5,6', \

width=8, \

datacolumn='data', \

intent='OBSERVE_TARGET#ON_SOURCE')

split(vis='msfiles/Elias24_continuum.ms', \

outputvis='msfiles/Elias24_avg_aux2.ms', \

timebin='30s', \

spw='11,15', \

datacolumn='data', \

intent='OBSERVE_TARGET#ON_SOURCE')

split(vis='msfiles/Elias24_continuum.ms', \

outputvis='msfiles/Elias24_avg_aux3.ms', \

timebin='30s', \

spw='2,7,8,9,10,12,13,14', \

width=16, \

datacolumn='data', \

intent='OBSERVE_TARGET#ON_SOURCE')

Step 5.3: Split to average

We start by removing previous ms files, just in case. Then, we split the original ms file we downloaded into 3 different ms files. The first ms file contains the spectral windows (spws) with 4 channels per spw, the second ms file contains the spws with 8 channels, and the third with 16 channels. In each split, we combine those channels into 1, and we also average in time.

This step can also be done with mstransform.

# Join together

concat(vis=['msfiles/Elias24_avg_aux1.ms', \

'msfiles/Elias24_avg_aux2.ms', \

'msfiles/Elias24_avg_aux3.ms'], \

concatvis='msfiles/Elias24_avg.ms', \

dirtol = '0.1arcsec', copypointing = False)

# Remove auxiliary files

os.system('rm -rf msfiles/Elias24_avg_aux1.ms')

os.system('rm -rf msfiles/Elias24_avg_aux2.ms')

os.system('rm -rf msfiles/Elias24_avg_aux3.ms')

# Check averaged Elias24

listobs('msfiles/Elias24_avg.ms')

Step 5.4: Concatenate

Now we take those 3 ms files, and we concatenate them into only 1 file. To save memory, we remove the intermediate auxiliary files.

When checking the new file with listobs, you will see that each spw has only 1 channel.

# Move phase center

fixvis(vis=ms_cont, \

outputvis=ms_cont, \

phasecenter='J2000 16h26m24.0782s -24d16m13.884s')

Step 5.5: Shift phase-center

The phase center of this observation is not in the center of the disk (you can confirm this downloading the continuum fits file form the DSHARP data release webpage).

We will move the phase center of the visibilities to match the location of the peak emission of the disk. This should be a good approximation for the disk center.

It is not mandatory to change the phase-center, and you can skip this step in your own templates.

# Mask geometry

inc, pa = 29., 45.7 # Huang et al. (2018), DSHARP II

pa_mask = str(pa)

# Mask semimajor axis size

mask_semimaj = 1.15

semimajor = str(mask_semimaj)+'arcsec'

semiminor = str(mask_semimaj*np.cos(np.deg2rad(inc)))+'arcsec'

# Mask center

center_ra = '16h26m24.0782s'

center_dec = '-24d16m13.884s'

center = center_ra + ' ' + center_dec

# Create masks

mask_LB = 'ellipse[[%s, %s], [%s, %s], %sdeg]' % \

(center_ra, center_dec, \

semimajor, semiminor, pa_mask)

res_mask_LB = 'annulus[[%s, %s], [%s, %s]]' % \

(center_ra, \

center_dec, \

'2.5arcsec', '3.5arcsec')

Step 5.6: Create a mask for your template.

Using the same information you put in the info script, create a mask for your template and to estimate the noise level. This masks are CASA regions.

# Clean

tclean_wrapper(vis='msfiles/'+prefix+'_avg.ms', \

imagename='images/'+prefix+'_im', \

imsize=2000, \

cellsize='0.003arcsec', \

mask=mask_LB, \

scales=[0, 3, 8, 13], \

robust=0.0, \

gain=0.05, \

smallscalebias=0.32, \

cyclefactor=1.75, \

threshold=str(4*2.05e-02)+'mJy', \

savemodel='modelcolumn', \

niter=10000, \

interactive=True)

Step 6.1: Image your template

Use the task tclean_wrapper to image the averaged ms file. Here you have the freedom to set the CLEAN parameters by yourself, but you should use the info script parameters (same cellsize, imsize, robust).

The parameter smallscalebias will force CLEAN to prefer the bigger scales before cleaning with point sources. The gain=0.05 and cyclefactor=1.75 are for more conservative iterations.

We will clean until we reach 4 sigma, such that afterwards we can apply the JvM correction.

Step 6.2: CLEAN

The CLEAN interactive viewer should look like this, with the mask already set. Press the green arrow iteratively, or press the blue arrow and wait until its finished.

Be patient. This step will take a long time.

Step 6.3: Estimate SNR and apply JvM correction

You have succesfully generated the images. Execute the estimate_SNR function and the JvM correction yo your newly created images. This will give you the angular resolution of your images, and the noise information.

Check your JvM corrected image with the CASA viewer task. Here we have changed the brightness scale.

exportfits(imagename='images/'+prefix+'_im.image', \

fitsimage='images/'+prefix+'_im_noJvM.fits', \