Protoplanetary Disk Simulation (CASA 3.4): Difference between revisions

From CASA Guides
Jump to navigationJump to search
No edit summary
 
(8 intermediate revisions by the same user not shown)
Line 2: Line 2:
[[Category: Simulations]]
[[Category: Simulations]]


* '''This is an advanced simulation tutorial.  New users are recommended to begin with the [[Simulation Guide for New Users (CASA 3.4)]].'''
* '''This guide is applicable to CASA version 3.4.  For older versions of CASA see [[PPdisk_simdata_(CASA_3.3)]].'''
* '''To create a script of the Python code on this page see [[Extracting scripts from these tutorials]].'''
* '''To create a script of the Python code on this page see [[Extracting scripts from these tutorials]].'''


== Protoplanetary disk ==
== Data ==
*[ftp://ftp.cv.nrao.edu/NRAO-staff/rindebet/input50pc_672GHz.fits This fits file] is a model of a protoplanetary disk from S. Wolf (If you use it for anything more than learning CASA, please cite [http://adsabs.harvard.edu/abs/2005ApJ...619.1114W Wolf & D'Angelo 2005]).


*simobserve and simanalyze version for CASA 3.4
For this CASA Guide we will use [ftp://ftp.cv.nrao.edu/NRAO-staff/rindebet/input50pc_672GHz.fits a protoplanetary disk model] from S. Wolf.  ''If you use this FITS data for anything more than learning CASA, please cite [http://adsabs.harvard.edu/abs/2005ApJ...619.1114W Wolf & D'Angelo 2005].''
__NOTOC__
 
<br>
== Script with Explanation ==
====Explanation of the script====
 
Set '''simobserve''' as the current task and reset all parameters.


=====Set sim_observe as current task and reset all parameters=====
<source lang="python">
<source lang="python">
# Setting everything in simobserve to original defaults
# In CASA
default("simobserve")
default("simobserve")
</source>  
</source>
=====Image coordinate system can be verified=====
 
Review the image coordinate system using task {{imhead}}.
 
<source lang="python">
<source lang="python">
# This reports image header parameters in the Log Messages window
# This reports image header parameters in the Log Messages window
imhead("input50pc_672GHz.fits")
imhead("input50pc_672GHz.fits")
</source>   
</source>   
=====Image center can be identified=====
 
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 [http://casa.nrao.edu/docs/CasaRef/CasaRef.html tool reference manual].
 
When data are being manipulated with tools the data file must be explicitly opened and closed.
 
<source lang="python">
<source lang="python">
# These are tools found in the CASA toolkit
# In CASA
# They are very useful, but the interface is not as straightforward as the tasks
# You can find the tool reference manual here: http://casa.nrao.edu/docs/CasaRef/CasaRef.html
# The following command is used to open an image (this is part of the image analysis toolkit)
# Whenever data are being viewed/manipulated by tools, what is being operated on needs to be explicitly opened
# and closed (i.e. an image, a table, etc.)
ia.open("input50pc_672GHz.fits")
ia.open("input50pc_672GHz.fits")
#  Out[9]: True
</source>
#
 
# Reports the length of each axis in the opened image
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'''.
ia.shape()
 
# Out[11]: [257L, 257L, 1L, 1L]
<source lang="python">
#
# In CASA
# This command converts from pixel (our source file) to world coordinates (something usable by simdata)
axesLength = ia.shape()
ia.toworld([128.5,128.5])
# Divide the first two elements of axesLength by 2.
#  Out[12]:
center_pixel = [ x / 2.0 for x in axesLength[:2] ]
#{'numeric': array([  4.71239120e+00,  -4.01423802e-01,  1.00000000e+00,
# Feed center_pixel to ia.toworld and and save the RA and Dec to ra_radians and dec_radians
#        6.72000001e+11])}
(ra_radians, dec_radians) = ia.toworld( center_pixel )['numeric'][:2]
#
# Formats the coordinate just converted into hms
qa.formxxx("4.71239120rad",format='hms',prec=5)
# Out[13]: '18:00:00.03052'
#
# Formats one of the other coordinates into dms
qa.formxxx("-0.401423802rad",format='dms',prec=5)
#  Out[14]: '-022.59.59.602743'
#
# Final housekeeping by closing the image tool
# The image tool will now be detached from the image
ia.close()
ia.close()
</source>
</source>


=====Brightness scale can be viewed with 'imstat' task=====
Use the '''qa''' tool to convert the image center from radians to sexagesimal coordinates.
 
<source lang="python">
<source lang="python">
# Default parameters are adequate for this
ra_hms  = qa.formxxx(str(ra_radians)+"rad",format='hms',prec=5)
imstat("input50pc_672GHz.fits")
dec_dms = qa.formxxx(str(dec_radians)+"rad",format='dms',prec=5)
# ...
</source>
#  'max': array([  6.52469971e-05]),
 
# ...
Let's call our project "psim2".  This defines the root prefix for any output files from simobserve.
# that's 0.0652 mJy/pixel.
</source>
=====Let's call our project psim2=====
<source lang="python">
<source lang="python">
default("simobserve")
project = "psim2"
# This defines the root prefix for any output files from simobserve
project           = "psim2"
</source>
</source>


=====We'll leave the sky model the way it is: simobserve will create psim2.skymodel CASA image since this model is a fits file, and most but not all of CASA routines can operate directly on fits=====
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">
skymodel           = "input50pc_672GHz.fits"
skymodel = "input50pc_672GHz.fits"
</source>  
</source>
=====We need to decide where to point the telescopeThe image is 2/3 arcsec in size, so we only need one pointing. We could put that in a text file ourself, or let simobserve create the ascii pointing file for us.=====
We will specify the sky position for the center of the observation and set the map size to the size of the model imageSince 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.
 
<source lang="python">
<source lang="python">
setpointings      =  True
setpointings      =  True
Line 82: Line 71:
mapsize            =  "0.76arcsec"
mapsize            =  "0.76arcsec"
</source>  
</source>  
=====The default pointingspacing is fine: we'll only fit one pointing in the small mapsize the default calculation maptype hexagonal is ok too since only one will fit anyway.=====


=====We do want to calculate visibilities in a measurement set: let's do a 20 min snapshot observation using the "out20" ALMA antenna configuration:=====
We do want to simulate an interferometric observation, so we set '''obsmode''' accordingly.  We'll set '''totaltime''' to a 20-minute snapshot observation.
 
<source lang="python">
<source lang="python">
obsmode            =  "int"
obsmode            =  "int"
Line 90: Line 79:
</source>
</source>


=====Use appropriate antenna configurations based on desired angular resolution (configuration 20 - alma.out20.cfg in this case - is the largest "compact" configuration)=====
We want to use the appropriate antenna configuration for the desired angular resolution.  Configuration 20, <tt>alma.out20.cfg</tt>, is the largest "compact" configuration.
 
<source lang="python">
<source lang="python">
antennalist        =  "alma.out20.cfg"
antennalist        =  "alma.out20.cfg"
</source>
Now run simobserve.
<source lang="python">
simobserve()
simobserve()
</source>
</source>


===== Deconvolve the visibilities back into an image=====
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.
 
<source lang="python">
<source lang="python">
default ("simanalyze")
default ("simanalyze")
project           = "psim2"
project = "psim2"
image             = True
image = True
# Prior image to use in clean
modelimage        =  "input50pc_672GHz.fits"
vis                =  project+".alma.out20.ms"
imsize            =  [192, 192] 
</source>
</source>


=====Specify number of iteration of cleaning task with proper threshold and weighting=====
We set '''modelimage''' to use the input FITS image when cleaning the simulated visibilities.  We set the image size to 192 pixels square.
 
<source lang="python">
<source lang="python">
niter              = 10000
modelimage = "input50pc_672GHz.fits"
threshold          = "1e-7Jy"
vis = project + ".alma.out20.ms"
weighting          "natural"   
imsize = [192, 192]  
</source>
</source>


=====We'd like to calculate a difference and fidelity image, and see some diagnostics:=====
Specify the number of iterations for cleaning, with proper threshold and weighting.
 
<source lang="python">
<source lang="python">
analyze            = True 
niter = 10000
threshold = "1e-7Jy"
weighting = "natural"   
</source>
</source>


=====And see the array but not the UV coverage:=====
We'd like to calculate a difference and fidelity image, but we don't need to see the uv coverage.
 
<source lang="python">
<source lang="python">
showuv             = False
analyze = True 
showresidual       = True   
showuv = False
showconvolved     = True
showresidual = True   
showconvolved = True
</source>
</source>


=====Plot both to the screen and the png files with lots of messages:=====
Plot both to the screen and PNG files with lots of messages.
 
<source lang="python">
<source lang="python">
graphics           = "both"
graphics = "both"
verbose           = True
verbose = True
overwrite = True
overwrite = True
</source>
</source>


===Run simdata===
Run simanalyze.
 
<source lang="python">
<source lang="python">
# This commands CASA to execute simanalyze
simanalyze()
simanalyze()
</source>
</source>
<br>
*Output results:


==Simulation Output==


{| style="border:1px solid #3366FF; " cellspacing=2
{| cellspacing=2
|Input:<br> [[File:Psim2.skymodel.png|300px]]
|Input:<br> [[File:Psim2.skymodel.png|300px]]
|Predict:<br> [[File:Psim2.alma.out20.observe.png|300px]]
|simobserve:<br> [[File:Psim2.alma.out20.observe.png|300px]]
|-
|-
|Image:<br> [[File:Psim2.alma.out20.image.png|300px]]
|simanalyze image:<br> [[File:Psim2.alma.out20.image.png|300px]]
|Analyze:<br> [[File:Psim2.analysis.png|300px]]
|simanalyze output:<br> [[File:Psim2.analysis.png|300px]]
|}
|}


{{Simulations Intro}}
{{Simulations Intro}}
{{Checked 3.4.0}}
{{Checked 3.4.0}}

Latest revision as of 15:37, 12 June 2012

Simulating Observations in CASA

Data

For this CASA Guide we will use a protoplanetary disk model from S. Wolf. 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("input50pc_672GHz.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("input50pc_672GHz.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.

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.

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.

skymodel = "input50pc_672GHz.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.

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.

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.

antennalist        =  "alma.out20.cfg"

Now run simobserve.

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.

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.

modelimage = "input50pc_672GHz.fits"
vis = project + ".alma.out20.ms"
imsize = [192, 192]

Specify the number of iterations for cleaning, with proper threshold and weighting.

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.

analyze = True  
showuv = False
showresidual = True  
showconvolved = True

Plot both to the screen and PNG files with lots of messages.

graphics = "both"
verbose = True
overwrite = True

Run simanalyze.

simanalyze()

Simulation Output

Input:
Psim2.skymodel.png
simobserve:
Psim2.alma.out20.observe.png
simanalyze image:
Psim2.alma.out20.image.png
simanalyze output:
Psim2.analysis.png

Simulating Observations in CASA

Last checked on CASA Version 3.4.0.