SIMIO-continuum - SIMulatIons to Observations

SIMIO-continuum is a collection of codes designed to help you compare your radiative transfer images with existing ALMA observations. With SIMIO-continuum you can generate your synthetic observations with different geometries, distances, and even including white noise.

Please check www.nicolaskurtovic.com/simio and the Contents section for a detailed description of the package. Although SIMIO-continuum is the official name, this page uses SIMIO-continuum or SIMIO indistinctively.

SIMIO-continuum

What is SIMIO-continuum?

SIMIO-continuum is a collection of functions and wrappers for CASA 5.6.2, designed to take a synthetic image representing a sky intensity distribution in millimeter and sub-millimeter wavelengths (such as a radiative transfer image from RADMC3D) and return the synthetic ALMA observation based on an existing dataset. By using an existing observation as a template, SIMIO keeps all the technical properties (number of antennas, time of observation, sky coordinates) and replaces the observed data with the visibilities of the synthetic image.

As the technical setup between the template and synthetic observation is the same, SIMIO allows for direct comparison. Therefore, with SIMIO we can answer the question: How would a synthetic model have looked if an existing ALMA observation had observed it?

b0a7874a882d4e01a500e62896c9a5e3

SIMIO was designed to compare ALMA observations of planet-forming disks with synthetic models. Nevertheless, the package works for any combination of synthetic model and ALMA continuum observations. SIMIO supports single synthetic images representing the intensity distribution of continuum emission. Dust polarization models or variable intensity as a function of frequency are not yet supported.

SIMIO is an alternative to generate synthetic observations that is as easy to use as convolving an image with a Gaussian beam but as robust as using a simulator tool like simobserve. This way, astronomers with little or no observational background can use it to generate their synthetic observations (see why going visibility based is more robust). Even though SIMIO can be used for non-ALMA experts, it also returns all the possible data products from a continuum observation, so experienced observers can also take them and analyze them with their own tools.

Compare directly to any ALMA observation:

SIMIO has several templates that you can use for your synthetic observation, and you can also add templates of existing ALMA observations by yourself. Each template will contain an archival observation of a planet-forming disk and the information about the source, to enable direct comparison if needed.

By default, SIMIO will take your synthetic image and calculate its visibilities as if your disk model had the same geometry (inclination and position angle) and the same distance of the template. The geometry and the source’s distance can be modified (see tutorials 3 and 4). You can also disable the geometry modification, and input your radiative transfer image with the desired geometry (see tutorial 2).

SIMIO generates the measurement set of the synthetic observation, but it can also generate further products such as: - Measurement Set (ms file). - Images in fits format (Beam convolved, PSF, model, residuals). - Visibility tables in .txt format. - Additional data figures, comparing your source to the template.

An example of how the young Solar System would have looked if DSHARP had observed it as HD163296 is shown in Figure 1, and also explained in tutorial 1.

331e4d0a8f4d465ebf4abbe8d33498a9

From sky to images

What does ALMA observes?

Let us say there is an interesting object in the sky, and we wish to study its spatial brightness distribution. For simplicity, let us assume we are observing a portion of the sky small enough such that the sky can be described as a flat surface (this is the case for most ALMA observations of planet-forming disks). Our object of study will have an intensity distribution over the 2D surface defined as $I(l,m)$ , where $(l,m)$ are the sky spatial coordinates. From an observational point of view, our goal is to recover the function $I$ as accurately as possible.

As ALMA is an interferometer, it does not directly observe the sky brightness distribution (as most optical and near-infrared cameras). Instead, ALMA observes the visibility function $V$ of the sky intensity, which is the Fourier Transform of the intensity distribution:

${V}(u, v) = \int \int {I}(l,m) \, e^{-2\pi i (ul + vm)}\,\text{d}l\,\text{d}m$

If we knew $V$ for every possible $(u,v)$ , then recovering $I$ would be as simple as calculating the inverse Fourier Transform of $V$ . With interferometers, that is not the case. For a given set of antennas, the baselines of an observation are all the possible combinations of antenna pairs. Each baseline $i$ will measure the value of $V_i$ in a given coordinate $(u_i,v_i)$ , therefore, an interferometer with $N$ baselines will sample the values $\{V(u_i, v_i)\}_{i=0}^N$ only for a discrete set of positions $\left\{(u_i,v_i)\right\}_{i=0}^N$ . This set of coordinates in the uv-plane will be called uv-coverage (see panel (c) of Figure 1 for an example of an incompletely sampled visibility space).

Reconstructing the sky brightness distribution

As we only have an incomplete measurement of the visibility function $V$ , we cannot directly recover the sky intensity distribution $I$ with an inverse Fourier Transform. Instead, we will need to assume values of $V$ in the spatial frequencies $(u,v)$ where we do not have measurements.

Although there are several different algorithms to reconstruct $I$ from an incompletely sampled $V$ , all of them rely on creating a model $I_{\text{mod}}$ . The visibility function of this model matches $\{V(u_i, v_i)\}_{i=0}^N$ by construction or by the minimization of a given variable (such as $\chi^2$ ). The whole process from observing to model reconstruction is shown in the following Figure 1:

c651ca9f725447a28f08ed9ceedb7dbb

Figure 1: An object in the sky has a particular visibility representation, given by its Fourier Transform. ALMA only samples the visibility representation in a discrete set of spatial frequencies, and we need to reconstruct a model to recover the sky brightness distribution.

Why do we convolve the models with Gaussians?

The longest baselines of an interferometer set the maximum angular resolution of a particular observation. An interferometer does not sample details in spatial frequencies higher than the longest baselines can access; thus, a given observation will have no information about them.

On the other hand, the reconstructed models have information stored in pixels (the minimum spatial unit of an image). These pixels, however, are much smaller than the physical angular resolution limit an observation can achieve, meaning that the information of intensity variation from one pixel to another is contained in a region of the visibility function that was not sampled by the observation. In other words, the observation does not strongly constrain what individual pixels do in a model.

One way to deal with this super-resolution behavior of the pixels in a model is to convolve it with a Gaussian representative of the observation PSF (point spread function). In this way, the information is deleted in spatial scales smaller than the angular resolution, and the features in the Gaussian convolved image represent the observation physical angular resolution.

The above-described process is what the CLEAN algorithm does. It creates a model from the observation, subtracts the visibility function of the model from the data to obtain the residuals, and then adds the residuals to the Gaussian convolved model. A diagram of this image reconstruction algorithm is shown in Figure 2:

0b87c76fd0f843ad8a0bdda16b0abf12

Figure 2: A model is convolved with a Gaussian to erase the spatial information of scales smaller than the angular resolution. Afterward, it is added to the inverse Fourier Transform of the data minus the model, creating the final interferometric image. For a simulated model, there is no residual to be added; therefore, it can immediately be convolved with a Gaussian and generate the synthetic observation.

The angular resolution is not uniquely determined

The angular resolution of an observation is not uniquely determined, as it is calculated from a weighted average of the baselines. Depending on how you weigh the baselines of your observation (more weight to high S/N baselines or more weight to extended baselines), the angular resolution of your reconstructed image can change dramatically.

One way to change the weight of the baselines is to change a variable called robust parameter, which is a parameter of the CLEAN algorithm. An example of the angular resolution change is the observation of the circumbinary disk CS Cha (Kurtovic et al. 2022). ALMA observed this disk for 6.5hrs in total. Here in Figure 3, you can see different images for the same disk, with the same observation, but with different robust parameters.

f1fdaf3c18654a429e528aa8a2ad9207

Figure 3: Images of CS Cha continuum emission from the same ALMA observation, imaged with different robust parameters shown in the upper left corner of each panel. The white ellipse at the bottom left corner shows the angular resolution of each image, while the numbers next to it show the half-width of the representative Gaussian used to convolve the model in units of milliarcsec. The scale bar represents 20au at the distance of the source.

Let us say you have a radiative transfer model of this circumbinary disk, and you want to compare it with the observation. The angular resolution changes by a factor of $\approx 3$ between the images shown in the gallery; all of them are valid images of the observation. If you were convolving your model with a Gaussian, which angular resolution would you use to compare?

Model reconstruction is a fundamental step of image reconstruction.

When comparing a simulation or a synthetic model to an interferometric observation, the Gaussian convolution is a good approximation of one of the steps of the imaging reconstruction algorithm, as shown in Figure 2. However, a fundamental step is the construction of the model itself, which will be dependent on the uv-coverage and the quality of the data.

The most robust way to compare simulations to observations is by calculating the Fourier Transform of the simulated image at the same frequencies measured by the uv-coverage of a given observation. Afterward, we can apply an image reconstruction algorithm and unveil what structures and properties would have been observed by a particular observational setup and image reconstruction.

Simulating with simobserve vs. SIMIO

The CASA software, which supports the ALMA observations, already has a built-in function to generate simulated observations starting from an intensity distribution image called simobserve. Even though simobserve works very well to test a particular observational setup, it can be very challenging to use it to compare with existing observations, as every single detail of an observation needs to be matched in order to obtain the same uv-coverage: The exact number and position of antennas, time on source, position of the source in the sky, time of observation, frequency coverage, and others.

Setting all the observational details manually in simobserve is useful when the science goal has to be tested in an observation that has never been done before. However, if the goal is to compare with an existing observation, manually setting all the details can be tedious, especially considering that the observation datasets already store all those technical properties. SIMIO-continuum is an alternative to simobserve to compare models with existing observations.

SIMIO takes an existing observation and then replaces the visibility data of the template with the visibilities of the input image. The Fourier Transform used in SIMIO is the same as CASA, so it is as robust as simobserve. The technical details are matched perfectly between synthetic and archival observation, and the newly generated synthetic dataset can be imaged with any algorithm.

An example of the advantage of SIMIO for comparing with existing observations comes from considering the observation log of TW Hya from Huang et al. (2018) shown in Table 1. With simobserve, we would have to generate each observation to compare with TW Hya. With SIMIO, we can take the existing observation and replace the data.

b827aa736d744e09a95ecd19fda3bb41

Table 1: Observation log of TW Hya from Huang et al. (2018).

How does SIMIO generates the observations?

Templates are the base of SIMIO

Instead of creating a new observation from scratch, SIMIO takes existing observations and replaces the data with the visibilities of a given model. By doing so, the synthetic observation of the input model will have the same visibility coverage as the template. In practice, the model is positioned in the sky at the exact coordinates of the template observation.

Generating a synthetic observation

The steps to replace the template visibilities with the input model visibilities are as follows:

A SIMIO object is created: The class simio_object sets all the necessary details for the synthetic observation. The object will get all the technical details of the selected template, and other observable parameters will be set too, such as geometry modifications, the distance of the observation, and observed flux. The tutorial section contains examples of how to modify those values when defining a simio_object.
Generate synthetic observation: A synthetic observation is created based on the template chosen for the simio_object. The workflow is shown in Figure 1, and it goes as follows:

2a) An empty image is created with tclean using the template’s information and with the same pixel size and image size as the input model. The input model values are copied into this empty image, which has the CASA format of an image.model.

2b) If a geometry change is needed, the uv-points of the template will be projected with the requested inclination and position angle. The Fourier Transform will be calculated over this new set of uv-points, thus avoiding modifying the image to apply a geometry modification. After that, the uv-points are deprojected with the same inclination and position angle, returning to their original position, but effectively having measured the Fourier Transform of the image.model in the projected space.

2c) The Fourier transform of the image.model is written into the synthetic observation, which is identical to the template observation except for the visibility data. Each spectral window is written separately.

0729c2e7a7b147178b637ee5b3073ec3

Figure 1: Steps followed by SIMIO to go from the input image to the synthetic observation. Instead of modifying the image geometry, the uv-points are projected/deprojected before and after the Fourier Transform calculation.

Imaging a synthetic observation

The image reconstruction of the synthetic observation can be done within SIMIO using the function easy_mod_tclean. This function is designed so people unfamiliar with CASA or ALMA data can generate their own images as easily as possible. Just tell easy_mod_tclean what simio_object will be imaged, set a stopping threshold if needed, and get the images done.

CASA API

Here is the description of all the functions used to generate a synthetic observation with SIMIO-continuum. These functions must be executed in a `CASA 5.6.2<https://casa.nrao.edu/casa_obtaining.shtml>`_ terminal interface.

Note

Please, check the tutorials to see how to utilize these functions. For simplicity, most of the SIMIO-continuum functions are made to run inside the simio_object, and therefore you will not interact directly with them.

Main functions

class simio_object(object_name, out_file_name, template, use_geom=True, distance=None, rescale_flux=None, pxsize_au=None, add_inc=0, add_pa=0, add_dRa=0, add_dDec=0)

Location: codes/simio_obj.py

The simio_object is the main object of the SIMIO package. It contains the functions and properties needed to generate the synthetic visibilities and images from a simulation.

Parameters

object_name – (str) Name of the project.
out_file_name – (str) Name of the RADMC3D .out file, or .npy file name.
template – (str) Template to be used as observation base.
use_geom – (bool) Set to True if you want to use the geometry of the template. If you set it to False then the parameters add_inc, add_pa, add_dRa, add_dDec are activated. Default: True.
distance – (float) Distance at which your model has to be positioned, in parsecs. If set to None, then the distance of the template will be used. Default: None.
rescale_flux – (float) Your model image is rescaled by a scalar, so that the total flux is rescale_flux. The units are Jy. If set to None, no flux rescaling is applied. Default: None.
pxsize_au – (float) Pixel size in au. If your input model is a .npy file, then this parameter is mandatory. It is not used if your file format is .out. Default: None.
add_inc – (float) Incline the source by this value, in degrees. Default: 0.
add_pa – (float) Rotate the source by this value, in degrees. Default: 0.
add_dRa – (float) Shift the source by this value in RA, in arcsec. Default: 0.
add_dDec – (float) Shift the source by this value in Dec, in arcsec. Default: 0.

add_noise(mod_ms, level='10.2mJy')

Location: codes/simio_clean.py

Wrapper for sm.setnoise from CASA. This function receives the name of the model measurement set (mod_ms from SIMIO tutorials), and returns a measurement set with the same name, but with added simple thermal noise.

Warning

The noise level in the measurement set will not be the same as you input in level. After succesful execution, generate an image to measure the noise level in the residuals image, and then run SIMIO again to iteratively find the correct level for the noise desired.

Parameters

mod_ms – (str) Name of the measurement set to be modified.
level – (str) Level of noise to be given to sm.setnoise, and passed directly to simplenoise.

Returns:

(int) Returns 1 if everything worked correctly. The noiseless measurement set will be copied into a file with the same name but ending in _no_noise.ms, while the mod_ms file will be modified to include the requested noise.

Imaging functions

easy_mod_tclean(simobj, interactive=False, remove_end=True, manual_threshold=str(0.024) + 'mJy')

Location: codes/simio_clean.py

Function wrapper of tclean, estimate SNR, JvM correction and delete wrapper. It uses the values from the template and simobj to fill the tclean_wrapper parameters. For a more customized clean, see custom_clean function, or tclean_wrapper.

Parameters

simobj – (simio_object) A simio object that already went through the get_mod_ms_ft function.
interactive – (boolean) Interactive clean. Recommended to set True. Default: False.
remove_end – (Boolean) If True, will remove the folder files after finishing the imaging. Default: True.
manual_threshold – Set the threshold for tclean. By default it cleans to 2sigma of DSHARP-like rms. Default: '2.4e-02mJy'.

Returns:

Fits files containing the reconstructed images, including the residuals, psf, JvM corrected image, and non-JvM corrected images.

custom_tclean(simobj, imsize, cellsize, robust, mask, threshold, scales=[0, 3, 8], gain=0.05, smallscalebias=0.45, cyclefactor=1.75, niter=10000, imagename=None, interactive=False, remove_end=True)

Location: codes/simio_clean.py

Function wrapper of tclean, estimate SNR, JvM correction and delete wrapper. It allows for a more customized clean compared to easy_mod_tclean. For more details on some of these parameters, check the tclean task in tclean documentation

Parameters

simobj – (simio_object) A simio object that already went through the get_mod_ms_ft function.
imsize – (int) Image size in pixels.
cellsize – (float) Pixel size, must be input in arcsec.
mask – (str) Mask for cleaning the emission, must be a CASA region format.
threshold – (float) Threshold for how deep the CLEAN should go, in mJy. For JvM corrected images, set the threshold to be 4 times the rms of the image. For model comparison with other models, you should clean up to 2 or 1 sigma.
scales – (list of int) Scales to use in multiscale, in pixels. Default: [0, 3, 8]
gain – (float) Fraction of the source flux to subtract out of the residual image for the CLEAN algorithm. Default: 0.05
smallscalebias – (float) Controls the bias towards smaller scales. Default: 0.45
cyclefactor – (float) Computes the minor-cycle stopping threshold. Default: 1.75
niter – (int) Total number of iterations. Default: 10000
imagename – (str) Sufix name for the images, it will be saved in the same folder as in default. Default: None
interactive – (boolean) Interactive clean. Recommended to set True. Default: False
remove_end – (boolean) If True, will remove the folder files after finishing the imaging. Default: None.

Returns:

Fits files containing the reconstructed images, including the residuals, psf, JvM corrected image, and non-JvM corrected images.

Additional Imaging functions

delete_wrapper(imagename)

Location: codes/simio_clean.py

Wrapper to delete the images generated by tclean.

Parameters: imagename – (str) Base name for the images to be deleted.

write_fits(im_base_name)

Location: codes/simio_clean.py

Given the im_base_name from tclean, it takes the products and write fits files of them.

Parameters: im_base_name – (str) Base name for the images to be written in fits format.

estimate_SNR(imagename, disk_mask, noise_mask)

Location: codes/simio_clean.py

Original from DSHARP.

Estimate peak SNR of source, given a mask that encompasses the emission and another annulus mask to calculate the noise properties.

Parameters

imagename – (str) Image name ending in .image.
disk_mask – (str) must be a CASA region format.
noise_mask – (str) Annulus to measure image rms, in the CASA region format, e.g. annulus[['0arcsec', '0arcsec'],['1arcsec', '2arcsec']].

create_dotmodel(simobj, imagename=None)

Location: codes/simio_clean.py

Function to create a .model image that mimics the .out or .npy input, with the coordinate information of the template.

Parameters

simobj – (simio_object) SIMIO object that will be used to generate the synthetic observation.
imagename – (str) Name of the image model to be generated.

Returns:

im_mod: (str) The name of the .model image generated.

Additional Visibility functions

change_geom(ms_file, inc=0.0, pa=0.0, dRa=0.0, dDec=0.0, datacolumn1='DATA', datacolumn2='DATA', inverse=False)

Location: codes/simio_ms2ascii.py

Changes the geometry of an observation, by inclining and rotating the uv-points themselfs. This function modifies the input ms_file.

Parameters

ms_file – (str) Name of the measurement set you want to incline, rotate or shift in physical space.
inc – (float) Inclination, in degrees. Default: 0.
pa – (float) Position angle, measured from north to east, in degrees. Default: 0.
dRa – (float) Shift in RA to be applied to the visibilities, in arcsec. Default: 0.
dDec – (float) Shift in Dec to be applied to the visibilities. in arcsec. Default: 0.
datacolumn1 – DATA or MODEL_DATA, column from where the data must be read. Default: DATA.
datacolumn1 – DATA or MODEL_DATA, column from where the data must be written. Default:DATA.
(bool) (inverse) – Set False to deproject, or False to project. Default: False.

Returns:

Returns True if everything worked correctly. The ms_file will have been modified with the new visibility geometry.

Masking functions

simio_object.get_mask(mask_semimajor=None, inc=None, pa=None)

Location: codes/simio_obj.py

Elliptical mask for CLEAN. The emission inside this mask will be cleaned. If no input is specified, the parameters of the template will be used. The output is a CASA region. See CASA Regions format for more information

Parameters

mask_semimajor – (float) Semimajor axis of the ellipse in arcsec. Default: None.
inc – (float) inclination of the ellipse in degrees. Default: None.
pa – (float) position angle of the ellipse, measured from the north to the east, or counter-clock wise, in degrees. Default: None.

Returns:

mask_obj: (str) elliptical mask. This is a CASA region.

simio_object.get_residual_mask(mask_rin=None, mask_rout=None)

Location: codes/simio_obj.py

Annulus mask to calculate the residuals properties. This mask is a circular annulus centered on the phase-center. The inner and outer radius should be set such that the mask does not include any real emission.

Parameters

mask_rin – (float) Inner radius of the annulus in arcsec. Default: None.
mask_rout – (float) Outer radius of the annulus in arcsec. Default: None.

Returns:

mask_res: (str) Annulus mask. This is a CASA region.

Installation

Requirements: Common Astronomy Software Applications (CASA), version 5.X or 6.X.

SIMIO has to be executed from the CASA terminal, therefore you need to install a CASA. Starting from verson 1.2, SIMIO works in CASA version 5.X or 6.X, although the most recommended versions are CASA 5.4.X or CASA 5.6.X.

SIMIO comes in a self-contained folder, and after installing CASA, no further installation of any package is needed. To use the SIMIO functions, execute the CASA software from the SIMIO folder, and then execute the SIMIO codes in the CASA terminal.

Check the download section of SIMIO in this page, or download it from the git-hub page. Once you have SIMIO in your computer, check tutorial 1 for step by step instructions of how to use it.

Download the package

Download the SIMIO-continuum package (from now on, referred to as SIMIO), and open the SIMIO folder. You will see the following:

codes: Where the functions and wrappers are located.
plots: After creating your synthetic observation, you can use the functions within this folder to generate figures and radial profiles. Check these codes if you do not have experience with fits files.
projects: Store your projects in this folder.
templates: Store your templates in this folder.
casa_examples: Example codes to run SIMIO are stored in this folder.
simio_casa.py: Example code to use SIMIO. Any SIMIO code should be run from this location (~/path_to_simio/simio/).

5d11cb37cb9949f49624390b66776218

Install a template

Each template is an archival ALMA observation adapted to run with SIMIO. By selecting a template, you will match the same uv-coverage, exposure time, position of the object in the sky, frequency bandwidth, and all the technical parameters of such observation.

You can download all the publicly available templates from this page. Move your template into the templates/ folder, such as the next example showing the HD163296 template.

cb76954cc1b249949460a22d68d15d83

Test if SIMIO is working correctly

This test uses data from tutorial 1 , which you can download from here

If you customize your SIMIO functions and you want to test if your modifications have affected the code, there are two major properties you should check: Flux conservation from the input model and size scaling as a function of distance. Let us use the Solar System from Bergez-Casalou et al. 2022 as the testing case for SIMIO functionalities.

5dd1aa769e8247c48171efdb65f2bb4c

Figure 1: The 1.3mm continuum emission from the young Solar System as simulated by Bergez-Casalou et al. 2022.

Flux conservation

SIMIO takes the input image and translates it to a “.model” image of CASA, before calculating the Fourier Transform with ft. During this translation the image is copied and flux-scaled following the inverse square law.

Let us test if this translation worked correctly. The original young Solar System model has a flux of 427.591784Jy when positioned at 1pc, which is the default distance from RADMC3D outputs and is the assumed distance for the input models of SIMIO (when the flux is not overriden, see tutorial 5). Let us generate the synthetic observation at the distance of HD163296, which is 100.966pc from GAIA DR3. At such distance, we expect a flux of:

$\frac{427.591784 \,\text{Jy}}{100.966^2} = 0.041944 \,\text{Jy} = 41.944 \,\text{mJy}$

# Import needed python packages
import sys
import os
import numpy as np
import matplotlib.pyplot as plt

# Get the current directory path
current_dir = os.getcwd()+'/'

####################################

# Import the analysis utils functions
sys.path.append(current_dir+'codes/analysis_scripts/')
import analysisUtils as au

# Import the simio object
execfile(current_dir+'codes/simio_obj.py')
# Import functions for uv-handling
execfile(current_dir+'codes/simio_ms2ascii.py')
# Import functions for imaging
execfile(current_dir+'codes/simio_clean.py')

###########################
# Solar System as HD163296
###########################

# Create a simio object.
simobj = simio_object(object_name  = 'SolarS_HD163296',
                      im_file_name = 'image_1300micron.out',
                      template     = 'HD163296',
                      use_tempgeom = True)

# Create the measurement file of your simio object
mod_ms = get_mod_ms_ft(simobj)

To generate mod_ms, SIMIO generated an image called “SolarS_HD163296_orig_model_im.fits”. The visibilities are calculated from this image. Let us use python to read the image and calculate the flux.

# Import astropy.io.fits to read the image in python
from astropy.io import fits

# Path to the image
model_path = 'projects/SolarS_HD163296/images/SolarS_HD163296_orig_model_im.fits'
# Open image
model_im  = fits.open(model_path)[0].data

# Calculate flux
print (np.sum(model_im))

# 0.04194438763274089

This value matches the expected value from the inverse square law. As the Fourier Transform is calculated from this image, we confirm that distance correction conserved the flux. Always check the project_orig_model_im.fits for this test.

Distance scaling

The same image from the previous part will be used to check distance scaling. The original young Solar System image has 800 pixels and 0.4au of pixel size, therefore the whole image is 320au in length. Let us check the translated image:

# Import astropy.io.fits to read the image in python
from astropy.io import fits

# Distance to HD163296
parallax = 9.90426503087695 # GAIA EDR3 mas
dist     = 1000. / parallax # pc
# Size of the input image
imsize_input = 320.

# Path to the image
model_path = 'projects/SolarS_HD163296/images/SolarS_HD163296_orig_model_im.fits'
# Read header
header  = fits.open(dir_cont)[0].header
# Calculate image size in arcsec
ra_ext  = 3600. * header_cont['CDELT1'] * header_cont['NAXIS1']
dec_ext = 3600. * header_cont['CDELT2'] * header_cont['NAXIS2']

# Calculate image size in au
imsize_au = dec_ext * dist
print (dec_ext * dist)
# 320.0045547174389

# Compare to real image size
ratio_im = dec_ext * dist / imsize_input # 320 is the size in au of the input image
print (ratio_im)
# 1.0000142334919966

The difference between the input image size and the model image is of the order of $10^{-5}$ in ratio, and it probably comes from the numerical systematics of transforming the pixel size in radians contained in header_cont['CDELT1'] to arcsec units and then to astronomical units. This way, you can test the distance scaling of your input model.

Considerations

Input image details

For each template, aim to have a pixel size at least $\times6\sim10$ times smaller than the angular resolution and image sizes larger than $\times6\sim10''$ . For optimal results, generate a different image for observations at different distances. In short, smaller pixel sizes and larger image sizes are always better. For a discussion about the input image size and the assumptions for a correct Fourier Transform calculation, we refer you to Tazzari et al. (2018).

CASA Warnings

CASA will write several warnings in the terminal while executing SIMIO. You can ignore them if they are included in the following list:

Leap Second: This SEVERE warning does not affect the results, unless you are working with VLBI or extremely high time precision data. Please check this page.

# SEVERE    MeasTable::dUTC(Double) (file ../../measures/Measures/MeasTable.cc, line 4290)  Leap second table TAI_UTC seems out-of-date. Until the table is updated (see the CASA documentation or your system admin), times and coordinates derived from UTC could be wrong by 1s or more.

Non-optimal architecture for synthetic measurement sets: As the templates are a combination of several observations, different spectral windows of the measurement sets have different frequency coverage and number of scans. Therefore, the Fourier Transform of the input model is calculated for each one separately (using the function split). The final measurement set is a concatenation of all the single spectral windows. The WARN will appear every time a new spectral window is concatenated.

The issue of a non-optimal architecture for the synthetic observation has no impact on the visibilities or the imaging products. A future version of SIMIO-continuum will explore a more efficient procedure to concatenate the synthetic observation.

# WARN  MSConcat::concatenate (file ../../ms/MSOper/MSConcat.cc, line 825)  Zero or negative scan numbers in MS. May lead to duplicate scan numbers in concatenated MS.

Published Papers that used SIMIO

Pinilla et al. (2022): Distributions of gas and small and large grains in the LkHa 330 disk trace a young planetary system.
Garrido-Deutelmoser et al. (2023): A Gap-sharing Planet Pair Shaping the Crescent in HD 163296: A Disk Sculpted by a Resonant Chain.
Gárate et al. (2023): Millimeter emission in photoevaporating disks is determined by early substructures.

Contribute or Report bugs

If you want to contribute to the code, please contact me at nicokurtovic at gmail.com. Similarly, if you encounter a bug or malfunction, contact me at my email or raise an issue in the github page.

Tutorial 1: From model to observations

Solar System as HD163296

Download this tutorial contents from here, including the project, models and script.

One key question to connect the origin of our Solar System to the general understanding of the planet-formation field is to estimate how would the Solar System have looked like if ALMA had observed it. In Bergez-Casalou et al. 2022, hydro-models were complemented with dust evolution models to answer this question.

By running radiative transfer codes with the dust distribution of the hypothetical Solar System, it is possible to obtain a prediction of the brightness of the Solar System planet-forming disk at 1.3mm wavelengths. The following figure shows such prediction, with each line showing the approximate orbit position of Jupiter, Saturn, Uranus, and Neptune, from the inside out, respectively.

80fd21213bea4e44ab6476ad62bcf3cd

Figure 1: The 1.3mm continuum emission from the young Solar System as simulated by Bergez-Casalou et al. 2022.

Radiative transfer images are the fundamental input for SIMIO-continuum. The package will take this image and generate a synthetic observation based on its brightness distribution. The details of how this is done is covered in “How does SIMIO generates the observations?” In the following steps, we will cover in detail the procedure to go from a model image to an ALMA observation.

Step 1: The package

Download the SIMIO-continuum package (from now on, referred to as SIMIO), and open the SIMIO folder. You will see the following:

codes: Where the functions and wrappers are located.
plots: After creating your synthetic observation, you can use the functions within this folder to generate figures and radial profiles. Check these codes if you do not have experience with fits files.
projects: Store your projects in this folder.
templates: Store your templates in this folder.
casa_examples: Example codes to run SIMIO are stored in this folder.
simio_casa.py: Example code to use SIMIO. Any SIMIO code should be run from this location (~/path_to_simio/simio/).

ca723b23c4bb45bf85e61534c424e5a7

Step 2: Include your template

We want to know how the Solar System would have looked like if it had been observed as HD163296 from DSHARP. Go into the templates/ folder and add the HD163296 template. You can download all the publicly available templates from this page.

Each template is an archival ALMA observation adapted to run with SIMIO. By selecting a template, you will match the same uv-coverage, exposure time, position of the object in the sky, frequency bandwidth, and all the technical parameters of such observation.

0c6fcc8c0e554c1886d667d32be981a3

Step 3: Create your project

Create a folder with the name of your project in the folder projects/. You can use any name you want for this folder. In this example, the name of our project will be “SolarS_HD163296” since we want to generate an observation of how the Solar System would look like if DSHARP had observed it at the distance and geometry of HD163296.

f2b3ff24ec954b2096ede6a8ffa9f4ee

Step 4: Prepare your project

Go inside your project folder. In this example, we are inside the folder SolarS_HD163296/. Create the folders “images”, “msfiles”, “uvtables”, and leave them empty. Add your radiative transfer image either in “.out” format (standard output format from RADMC3D), or “.npy” format. In this example, our radiative transfer image is “image_1300micron.out”, shown at the beginning of this tutorial. Alternatively, a file can be saved to “.npy” format by saving a NumPy array with the function np.save(array). Therefore, any image stored as a NumPy matrix can also be input in SIMIO.

Note: The input image must be as big as the field of view, and the pixel size must be at least five times smaller than the highest angular resolution. Bigger images and smaller pixel sizes will produce a more stable Fourier Transform of the models by reducing the short baselines artifacts and smoothing the brightness difference from pixel to pixel.

SIMIO will use the folders you just created to store:

images: The fits files generated with SIMIO will be stored in this folder. Check it after running the code in Step 5.
msfiles: Your synthetic observation will have its own measurement file, where the visibilities are stored. You will be able to find that file in this folder.
uvtables: (Being implemented) The visibilities will also be given in “.txt” format, which you can further use to analyze with alternative tools, such as frank or galario.

e91840fb9403404d8aed3109aa500630

Step 5: Run SIMIO

Go back to the initial SIMIO folder. Open simio_casa.py and open CASA 5.6.2 in a terminal. Note: SIMIO should work on any CASA 5.X version, but 5.4.X or 5.6.X should be preferred, as most of the testing has been done in those. You can find the links to download the CASA software here.

The first part of the simio_casa.py code will import all the necessary python packages in CASA and set the path to the SIMIO folder. These libraries are already included in the python of CASA, and you do not need to install them separately.

# Import needed python packages
import sys
import os
import numpy as np
import matplotlib.pyplot as plt

# Get the current directory path
current_dir = os.getcwd()+'/'

The second part of the code is needed to load the SIMIO functions. Each subcode needs functions from the previous, and so it is necessary to execute them in the correct order. The analysis utils of CASA are included in the SIMIO package for self-containing purposes.

You should not need to change anything in the import and execfile blocks.

# Import the analysis utils functions
sys.path.append(current_dir+'codes/analysis_scripts/')
import analysisUtils as au


# Import the simio object
execfile(current_dir+'codes/simio_obj.py')
# Import functions for uv-handling
execfile(current_dir+'codes/simio_ms2ascii.py')
# Import functions for imaging
execfile(current_dir+'codes/simio_clean.py')

Create your simio_object, the primary object of the SIMIO package. This object will contain all the needed information to generate your synthetic observation.

object_name: Write the name of your project, which you created in steps 3 and 4.
im_file_name: Name of your radiative transfer image, including the format (.out or .npy). Note: If your image file is .npy, you will need to specify the pixel size. Check the other tutorials for an example.
template: Name of the template you want to mimic. In this example, HD163296.
use_tempgeom: Set to True if you want SIMIO to incline and rotate your image with the inclination and position angle of the template, under the assumption that your input is face-on. If your image already has the geometry you desire, and you do not want SIMIO to apply any geometric change, then set it to False. Note You can also use just a rotation or incline. Please check the other tutorials.

###########################
# Solar System as HD163296
###########################

# Create a simio object.
simobj = simio_object(object_name  = 'SolarS_HD163296',
                      im_file_name = 'image_1300micron.out',
                      template     = 'HD163296',
                      use_tempgeom = True)

The function get_mod_ms_ft will take your image and generate the measurement set (which is the object that contains all the information of an ALMA observation) as if it had been observed with the same technical setup of your template. Running this function will give you the path to the generated measurement set in the variable mod_ms. The time this function takes to run will depend on the computer and template. As a reference, most templates take ~4min on an average laptop.

This line is where the Fourier Transform of your model is computed. For each spectral window of the real observation, SIMIO will replace the DATA COLUMN with the visibilities of your model, therefore matching exactly the uv-coverage. In other words, you will have the same angular resolution. The function get_mod_ms_ft will use the Fourier Transform from CASA.

# Create the measurement file of your simio object
# Can take several minutes
mod_ms = get_mod_ms_ft(simobj)

Now that the measurement set has been generated, the following step is to create the images of your synthetic observation. We will do this with the CLEAN algorithm.

Create the masks to CLEAN the synthetic observation. These masks will be loaded into the simobjand also returned as a string. They will have the standard of CASA Regions.

mask_obj: An elliptical mask with the geometry of the template. You can set the semi-major axis of the ellipse in units of arcsec.
mask_res: Generates an annulus mask needed to calculate the background properties of your image.

# Create Masks
mask_obj = simobj.get_mask(mask_semimajor=0.65) # in arcsec
mask_res = simobj.get_residual_mask()

Now we need to generate the images of your synthetic observation. The function easy_mod_tclean is a wrapper of the function tclean for simobj. It runs the CLEAN algorithm over your generated observation. Note: If you downloaded this project in the link at the beginning of the tutorial, you will notice that the images for the Solar System are already stored in SolarS_HD163296/images/. Running this tutorial will overwrite them, and th new images will look exactly the same.

Set interactive to True to check if the mask includes all the emission. Press the green arrow (in the Next Action section) to start a cycle of the cleaning process, or press the blue arrow and wait until it is done. Note: Use the blue arrow only if you have tested convergence of the current threshold.

d7432500918647529cd97bec37713615

The image looks structured and noisy, as this is the “dirty image” resulting from convolving the input image with the dirty beam. Please check the imaging tutorial of CASA for more information. Depending on your computer and the template, cleaning the image can take a few minutes to a fraction of an hour. Be patient.

Step 6: Check the results

After finishing the cleaning process, go back to the project folder to check the SIMIO products. The execution of get_mod_ms_ft generated a measurement set in the msfiles folder (the one you created in step 4). After executing easy_mod_tclean at the end of Step 5, you will get the products in the images folder.

The images will be named by your project name, plus a suffix. Each image is:

project_im.fits: Beam convolved image, with the JvM correction (Czekala et al. 2021). This image is how your source would look if ALMA had observed it with the same observational setup as the template.
project_im_model.fits: The model image, generated by the CLEAN algorithm as a description of the visibilities of your source in the sky plane.
project_im_noJvM.fits: Beam convolved image, without the JvM correction.
project_im_psf.fits: The PSF of the observation.
project_im_residual.fits: The residuals of the CLEAN algorithm. This residual image will be very structured, as the generated ms file does not contain noise. However, it should be negligible in total flux.

3494a11dd1024ae2936051789cf1c6ef

The standard imaging product of ALMA is the beam convolved image, which is called “project_im.fits” or “SolarS_HD163296_im.fits” in this tutorial. You can visualize that image with the task casaviewer from CASA or the software DS9. You can also open those files in python with the package astropy (Check the plots folder for examples).

If the Solar System was 1Myr old, located at the distance of HD163296, and had been observed by the DSHARP survey, then it would look like the following figure, where it is compared to the template disk.

5bd88d7a686d4950bf1d96fd0341d235

That’s all! :D

Now you have a synthetic observation of the Solar System, if it was located at the same distance of HD163296, at the same position in the sky, and if DSHARP had observed it. You can see that the Solar System would have looked like a compact disk with a cavity when observed through the eyes of ALMA.

You can take all the products and study them however you want: In the visibility plane or the image plane. You can also use the measurement set to generate images with different angular resolutions. Check the other tutorials to see what else you can do with SIMIO.

Tutorial 2: Modify a model geometry

Download this tutorial contents from here, including the project, models and script.

The radiative transfer model of a disk can be dependent on the inclination of the disk relative to the observer’s line of sight. In the previous tutorial, the initial model was face-on, and SIMIO later modified the geometry (inclination and position angle), which is only correct if we assume that the changes in inclination will conserve the flux.

It is also possible to change individual geometrical parameters or change none. In this tutorial, we want to compare a disk from Garate et al. (subm.) with Elias 24 without changing the inclination but modifying the position angle.

An inclined disk model

The model of the disk was generated with the same inclination of Elias 24 and with a semi-major axis along the x-axis. We would like to have the disk with the same position angle as Elias 24.

87af4cded88a4fc19bb2c4a551718254

Create SIMIO object

First of all, load the codes for SIMIO

# Import needed python packages
import sys
import os
import numpy as np
import matplotlib.pyplot as plt

# Get the current directory path
current_dir = os.getcwd()+'/'

####################################

# Import the analysis utils functions
sys.path.append(current_dir+'codes/analysis_scripts/')
import analysisUtils as au

# Import the simio object
execfile(current_dir+'codes/simio_obj.py')
# Import functions for uv-handling
execfile(current_dir+'codes/simio_ms2ascii.py')
# Import functions for imaging
execfile(current_dir+'codes/simio_clean.py')

Now, we will create the SIMIO object. Here, the name of the project is “photo_elias24_13”. We want to use the observation of Elias 24 as a template, but we do not want to use the template geometry. Therefore, we set use_tempgeom to False, and we tell SIMIO to add the model a position angle to match the one of Elias 24.

After that, we generate the observation with get_mod_ms_ft.

# Create a simio object.
simobj = simio_object(object_name  = 'photo_elias24_13',
                      im_file_name = 'continuum_13_inc29.out',
                      template     = 'Elias24',
                      add_pa       = 45.7+90,
                      use_tempgeom = False)

# Create the measurement file of your simio object, and get the path.
# Can take several minutes
mod_ms = get_mod_ms_ft(simobj)

Generate the images

With the observation already created, now the only thing left is to generate the images. We will create an elliptical mask slightly larger than the size of our model (0.52arcsec of radius) and with the geometry of Elias24.

After setting the masks, easy_mod_tclean will generate the images.

# Create a mask for your system, and one to measure the residuals
mask_obj = simobj.get_mask(mask_semimajor=0.52, inc=29., pa=45.7)
mask_res = simobj.get_residual_mask()

# Generate image for your simio object.
# Can take several minutes, maybe an hour. Depends on your computer
easy_mod_tclean(simobj, interactive=True)

964103fd6d764518b1963624e6bcde0d

Tutorial 3: Add noise to your observation

Download this tutorial contents from here, including the project, models and script.

Observational noise is a fundamental property to consider when generating predictions and comparing the detectability of different models. The CASA software has algorithms to include simple thermal noise, which are compatible with the SIMIO generated measurement sets.

Include thermal noise in your observation

We will start from the same model as in tutorial 2.

0ad4346105274f54a68c704c212046d3

Create SIMIO-continuum object

Load the codes for SIMIO-continuum and rotate the model to match the position angle of Elias24

# Import needed python packages
import sys
import os
import numpy as np
import matplotlib.pyplot as plt

# Get the current directory path
current_dir = os.getcwd()+'/'

####################################

# Import the analysis utils functions
sys.path.append(current_dir+'codes/analysis_scripts/')
import analysisUtils as au

# Import the simio object
execfile(current_dir+'codes/simio_obj.py')
# Import functions for uv-handling
execfile(current_dir+'codes/simio_ms2ascii.py')
# Import functions for imaging
execfile(current_dir+'codes/simio_clean.py')

####################################

# Create a simio object.
simobj = simio_object(object_name  = 'photo_elias24_13',
                      im_file_name = 'continuum_13_inc29.out',
                      template     = 'Elias24',
                      add_pa       = 45.7+90,
                      use_tempgeom = False)

# Create the measurement file of your simio object, and get the path.
# Can take several minutes
mod_ms = get_mod_ms_ft(simobj)

The function get_mod_ms_ft generated the observation. Now, we will use add_noise to add thermal noise to your observation. This SIMIO function is just a wrapper of sm.setnoise from CASA, and we refer you to that function for further details.

# Add thermal noise
add_noise(mod_ms, level='10.2mJy')

The function add_noise will generate a new measurement set with noise. Be careful with your level input. The number you give in level is directly passed to sm.setnoise, and is not exactly the same you will get as background rms in your reconstructed image. You will probably need to try a couple of times before finding the correct input number that will return the desired noise level.

Generate the images

The images are generated in the same way as the noiseless images. Just run easy_mod_tclean as we did in the previous tutorials.

# Create a mask for your system, and one to measure the residuals
mask_obj = simobj.get_mask(mask_semimajor=0.52, inc=29., pa=45.7)
mask_res = simobj.get_residual_mask()

# Generate image for your simio object.
# Can take several minutes, maybe an hour. Depends on your computer
easy_mod_tclean(simobj, interactive=True)

715004e3379b42f0a0e9b8dcaf2db10b

Here the brightness scale was modified to emphasize the thermal noise of the generated images.

Tutorial 4: Change the distance of your source

Download this tutorial contents from here, including the project, models and script.

In tutorial 1, we tested how the synthetic Solar System would have looked if it had been observed at the same distance of HD163296 (101pc from GAIA DR3). At that distance, the disk is observed with a cavity and a single ring. What features are recoverable if the Solar System is located farther away? How far away can we place the Solar System before we stop detecting the cavity?

Or, more generally: How far away can your disk be from Earth before a feature of your model is no longer detectable?

Changing the distance of your observation

Radiative transfer models are usually calculated as if the object was at a distance of 1pc. Inside the SIMIO functions, the flux from your model is scaled following the inverse square law, and this is done with the distance of the chosen template.

It could be the case that you want to test your substructure recovery at a specific distance, not necessarily the distance of the template. For example, let us choose HD163296 as a representative high angular resolution observation. How would the Solar System look with this observational setup if it were at different distances from Earth.

Let us begin by calling SIMIO.

# Import needed python packages
import sys
import os
import numpy as np
import matplotlib.pyplot as plt

# Get the current directory path
current_dir = os.getcwd()+'/'

####################################

# Import the analysis utils functions
sys.path.append(current_dir+'codes/analysis_scripts/')
import analysisUtils as au

# Import the simio object
execfile(current_dir+'codes/simio_obj.py')
# Import functions for uv-handling
execfile(current_dir+'codes/simio_ms2ascii.py')
# Import functions for imaging
execfile(current_dir+'codes/simio_clean.py')

Let us create the simobj with a specific distance

# Create a simio object.
simobj = simio_object(object_name  = 'SolarS_100pc',
                      im_file_name = 'image_1300micron.out',
                      template     = 'HD163296',
                      use_tempgeom = False,
                      distance=100.)

Here, we specified that we do not want to use the template geometry cause we would like to see the Solar System in a face-on configuration, and the distance must be input in float format with pc unit. In this example, we chose 100pc as the distance.

After that, just run the SIMIO noise function and imaging function.

# Create the measurement file of your simio object, and get the path.
mod_ms = get_mod_ms_ft(simobj)

# Add noise
add_noise(mod_ms, level='10.2mJy')

# Create a mask for your system, and one to measure the residuals
mask_obj = simobj.get_mask(mask_semimajor=0.52, inc=29., pa=45.7)
mask_res = simobj.get_residual_mask()

# Generate image for your simio object.
easy_mod_tclean(simobj, interactive=True, manual_threshold='4.8e-02mJy')

The noise level included will return an image noise of about 0.019mJy/beam, similar to DSHARP observations.

Image at any distance you want

The 100pc can be changed to any distance you want. The next figure shows the Solar System reconstructed with SIMIO at distances ranging from 100pc to 800pc. You will need an individual project for each one of those synthetic observations.

1792819157ae4168a4580c497a0b670d

SIMIO also returns the original model used to calculate the visibilities. Here you can check the reconstructed images against the original model:

d28ff8d24838468cac89d4ca1c225936

You can see that the noise level makes it very challenging to recover the inner cavity ring of the Solar System. You can try changing the noise level and analyzing the observation in the visibilities to see what kind of observation you would need to recover such features.

Tutorial 5: Change the flux of your source

Download this tutorial contents from here, including the project, models and script.

Sometimes a radiative transfer model will give you the correct contrast between emitting regions, but fine-tuning the total flux to match an observation can be challenging. SIMIO allows you to set the observed flux of your source independently from the distance at which the source is located.

Changing the observed flux of your observation

The observed flux can be specified when setting up the simobj. Let us begin by importing the SIMIO functions.

# Import needed python packages
import sys
import os
import numpy as np
import matplotlib.pyplot as plt

# Get the current directory path
current_dir = os.getcwd()+'/'

####################################

# Import the analysis utils functions
sys.path.append(current_dir+'codes/analysis_scripts/')
import analysisUtils as au

# Import the simio object
execfile(current_dir+'codes/simio_obj.py')
# Import functions for uv-handling
execfile(current_dir+'codes/simio_ms2ascii.py')
# Import functions for imaging
execfile(current_dir+'codes/simio_clean.py')

Let us create the simobj with a specific flux.

# Create a simio object.
simobj = simio_object(object_name  = 'SolarS_100pc',
                      im_file_name = 'image_1300micron.out',
                      template     = 'HD163296',
                      rescale_flux = 0.080)

The parameter rescale_flux receives the desired observed flux at the distance of the source in units of Jy. SIMIO will scale the model flux to match this value before the Fourier Transform is calculated.

# Create the measurement file of your simio object, and get the path.
mod_ms = get_mod_ms_ft(simobj)

# Add noise
add_noise(mod_ms, level='10.2mJy')

# Create a mask for your system, and one to measure the residuals
mask_obj = simobj.get_mask(mask_semimajor=0.65)
mask_res = simobj.get_residual_mask()

# Generate image for your simio object.
easy_mod_tclean(simobj, interactive=True, manual_threshold='5.7e-02mJy')

Here we are cleaning the image until we reach a threshold of 3 sigma. The rescaled image will have a flux of 80mJy, higher than the nominal flux of 42mJy when the image has no rescaling applied. You can check the flux with the function estimate_SNR.

2060635f74de46e691777dc64e0510a9

Support

If you are having issues, please open a issue on the GitHub page, or send me a message to my email address.

License

The project is licensed under the MIT license.