Protoplanetary Disk Simulation (CASA 6.6.1): Difference between revisions
Created page with " Category: Simulations == About == This is an advanced simulation tutorial. New users are recommended to begin with reading the [https://casadocs.readthedocs.io/en/v6.6.1/notebooks/simulation.html CASA Docs pages on Simulations]. To create a script of the Python code on this page see Extracting scripts from these tutorials. == Data == For this CASA Guide we will use a protoplanetary disk model from S. Wolf. Get the data [http://www.cv.nrao.edu/~awootten/mma..." |
|||
| Line 52: | Line 52: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
ra_hms = qa.formxxx(str(ra_radians)+"rad",format='hms',prec=5) | ra_hms = qa.formxxx(str(ra_radians)+"rad",format='hms',prec=5) | ||
dec_dms = qa.formxxx(str(dec_radians)+"rad",format='dms',prec=5) | dec_dms = qa.formxxx(str(dec_radians)+"rad",format='dms',prec=5) | ||
| Line 58: | Line 59: | ||
Let's call our project "psim2". This defines the root prefix for any output files from simobserve. | Let's call our project "psim2". This defines the root prefix for any output files from simobserve. | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
project = "psim2" | project = "psim2" | ||
</source> | </source> | ||
| Line 63: | Line 65: | ||
We'll set '''skymodel''' to the FITS file downloaded, above, and leave all '''skymodel''' subparameters at their default values. simobserve will create CASA image <tt>psim2.skymodel</tt>. | We'll set '''skymodel''' to the FITS file downloaded, above, and leave all '''skymodel''' subparameters at their default values. simobserve will create CASA image <tt>psim2.skymodel</tt>. | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
skymodel = "ppdisk672_GHz_50pc.fits" | skymodel = "ppdisk672_GHz_50pc.fits" | ||
</source> | </source> | ||
| Line 69: | Line 72: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
setpointings = True | setpointings = True | ||
direction = "J2000 18h00m00.031s -22d59m59.6s" | direction = "J2000 18h00m00.031s -22d59m59.6s" | ||
| Line 77: | Line 81: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
obsmode = "int" | obsmode = "int" | ||
totaltime = "1200s" | totaltime = "1200s" | ||
| Line 84: | Line 89: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
antennalist = "alma.out20.cfg" | antennalist = "alma.out20.cfg" | ||
</source> | </source> | ||
| Line 89: | Line 95: | ||
We do not want to simulate observations with thermal noise in this example, so we set: | We do not want to simulate observations with thermal noise in this example, so we set: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
thermalnoise = '' | thermalnoise = '' | ||
</source> | </source> | ||
| Line 95: | Line 102: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
simobserve() | simobserve() | ||
</source> | </source> | ||
| Line 101: | Line 109: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
default ("simanalyze") | default ("simanalyze") | ||
project = "psim2" | project = "psim2" | ||
| Line 109: | Line 118: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
modelimage = "ppdisk672_GHz_50pc.fits" | modelimage = "ppdisk672_GHz_50pc.fits" | ||
vis = project + ".alma.out20.ms" | vis = project + ".alma.out20.ms" | ||
| Line 117: | Line 127: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
niter = 10000 | niter = 10000 | ||
threshold = "1e-7Jy" | threshold = "1e-7Jy" | ||
| Line 125: | Line 136: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
analyze = True | analyze = True | ||
showuv = False | showuv = False | ||
| Line 134: | Line 146: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
graphics = "both" | graphics = "both" | ||
verbose = True | verbose = True | ||
| Line 142: | Line 155: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | |||
simanalyze() | simanalyze() | ||
</source> | </source> | ||
Revision as of 19:39, 8 August 2025
About
This is an advanced simulation tutorial. New users are recommended to begin with reading the CASA Docs pages on Simulations.
To create a script of the Python code on this page see Extracting scripts from these tutorials.
Data
For this CASA Guide we will use a protoplanetary disk model from S. Wolf. Get the data here. See the Simulation Inputs page for this and other sample images. If you use this FITS data for anything more than learning CASA, please cite Wolf & D'Angelo 2005.
Script with Explanation
Set simobserve as the current task and reset all parameters.
# In CASA
default("simobserve")
Review the image coordinate system using task imhead.
# This reports image header parameters in the Log Messages window
imhead("ppdisk672_GHz_50pc.fits")
We now use the ia (image analysis) and qa (units and quantities) tools from the CASA Toolkit to find the image center. In comparison to tasks, tools are a more advanced way of manipulating data in CASA. You can learn more about tools using the tool reference manual.
When data are being manipulated with tools the data file must be explicitly opened and closed.
# In CASA
ia.open("ppdisk672_GHz_50pc.fits")
Next, get the right ascension and declination of the image center. We get the number of pixels along each axis using ia.shape. Then, we get the RA and Dec values for the center pixel using ia.toworld.
# In CASA
axesLength = ia.shape()
# Divide the first two elements of axesLength by 2.
center_pixel = [ x / 2.0 for x in axesLength[:2] ]
# Feed center_pixel to ia.toworld and and save the RA and Dec to ra_radians and dec_radians
(ra_radians, dec_radians) = ia.toworld( center_pixel )['numeric'][:2]
ia.close()
Use the qa tool to convert the image center from radians to sexagesimal coordinates.
# In CASA
ra_hms = qa.formxxx(str(ra_radians)+"rad",format='hms',prec=5)
dec_dms = qa.formxxx(str(dec_radians)+"rad",format='dms',prec=5)
Let's call our project "psim2". This defines the root prefix for any output files from simobserve.
# In CASA
project = "psim2"
We'll set skymodel to the FITS file downloaded, above, and leave all skymodel subparameters at their default values. simobserve will create CASA image psim2.skymodel.
# In CASA
skymodel = "ppdisk672_GHz_50pc.fits"
We will specify the sky position for the center of the observation and set the map size to the size of the model image. Since the model image is 2/3 arcseconds across, we should only need one pointing. In this case, pointingspacing and maptype can be left at their default values.
# In CASA
setpointings = True
direction = "J2000 18h00m00.031s -22d59m59.6s"
mapsize = "0.76arcsec"
We do want to simulate an interferometric observation, so we set obsmode accordingly. We'll set totaltime to a 20-minute snapshot observation.
# In CASA
obsmode = "int"
totaltime = "1200s"
We want to use the appropriate antenna configuration for the desired angular resolution. Configuration 20, alma.out20.cfg, is the largest "compact" configuration. A list of configuration files available in CASA 5.4 is available here.
# In CASA
antennalist = "alma.out20.cfg"
We do not want to simulate observations with thermal noise in this example, so we set:
# In CASA
thermalnoise = ''
Now run simobserve.
# In CASA
simobserve()
Now that we've simulated the visibility measurements, we want to generate an image from the simulated data. simanalyze makes this easy. We begin by setting project to the same prefix used in simobserve and setting image to True.
# In CASA
default ("simanalyze")
project = "psim2"
image = True
We set modelimage to use the input FITS image when cleaning the simulated visibilities. We set the image size to 192 pixels square.
# In CASA
modelimage = "ppdisk672_GHz_50pc.fits"
vis = project + ".alma.out20.ms"
imsize = [192, 192]
Specify the number of iterations for cleaning, with proper threshold and weighting.
# In CASA
niter = 10000
threshold = "1e-7Jy"
weighting = "natural"
We'd like to calculate a difference and fidelity image, but we don't need to see the uv coverage.
# In CASA
analyze = True
showuv = False
showresidual = True
showconvolved = True
Plot both to the screen and PNG files with lots of messages.
# In CASA
graphics = "both"
verbose = True
overwrite = True
Run simanalyze.
# In CASA
simanalyze()
Simulation Output
Input:
|
simobserve:
|
simanalyze image:
|
simanalyze output:
|
