PPdisk simdata (CASA 3.1): Difference between revisions
From CASA Guides
Jump to navigationJump to search
Line 8: | Line 8: | ||
<br> | <br> | ||
====Explanation of the script==== | |||
Set simdata2 as current task and reset all parameters | =====Set simdata2 as current task and reset all parameters===== | ||
default("simdata2") | default("simdata2") | ||
Specify sky model image | =====Specify sky model image===== | ||
modelimage = "input50pc_672GHz.fits" | modelimage = "input50pc_672GHz.fits" | ||
Image coordinate system can be verified | =====Image coordinate system can be verified===== | ||
imhead("input50pc_672GHz.fits") | imhead("input50pc_672GHz.fits") | ||
Image center can be identified | =====Image center can be identified===== | ||
# ia.open("input50pc_672GHz.fits") | # ia.open("input50pc_672GHz.fits") | ||
# ia.shape() | # ia.shape() | ||
Line 27: | Line 27: | ||
# '-022.59.59.602743' | # '-022.59.59.602743' | ||
# ia.done() | # ia.done() | ||
Brightness scale can be viewed with 'imstat' task | =====Brightness scale can be viewed with 'imstat' task===== | ||
# imstat("input50pc_672GHz.fits") | # imstat("input50pc_672GHz.fits") | ||
# ... | # ... | ||
Line 33: | Line 33: | ||
# ... | # ... | ||
# that's 0.0652 mJy/pixel. | # that's 0.0652 mJy/pixel. | ||
Let's leave the brightness of the image as it is | =====Let's leave the brightness of the image as it is===== | ||
inbright = "unchanged" | inbright = "unchanged" | ||
Let's call our project psim2 | =====Let's call our project psim2===== | ||
project = "psim2" | project = "psim2" | ||
We'll leave the sky model the way it is -- simdata2 will create psim2.skymodel CASA image since this | =====We'll leave the sky model the way it is -- simdata2 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 | model is a fits file, and most but not all of CASA routines can operate directly on fits===== | ||
modifymodel = False | modifymodel = False | ||
skymodel = "input50pc_672GHz.fits" | skymodel = "input50pc_672GHz.fits" | ||
We need to decide where to point the telescope. The image is 2/3 arcsec in size, so we only need | =====We need to decide where to point the telescope. The image is 2/3 arcsec in size, so we only need | ||
one pointing. We could put that in a text file ourself, or let simdata2 create the ascii | one pointing. We could put that in a text file ourself, or let simdata2 create the ascii | ||
pointing file for us. | pointing file for us.===== | ||
setpointings = True | setpointings = True | ||
direction = "J2000 18h00m00.031s -22d59m59.6s" | direction = "J2000 18h00m00.031s -22d59m59.6s" | ||
mapsize = "0.76arcsec" | mapsize = "0.76arcsec" | ||
The default pointingspacing is fine - we'll only fit one pointing in the small mapsize | =====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. | the default calculation maptype hexagonal is ok too since only one will fit anyway.===== | ||
We do want to calculate visibilities in a measurement set: | =====We do want to calculate visibilities in a measurement set: | ||
let's do a 20 min snapshot observation using out20 configuration: | let's do a 20 min snapshot observation using out20 configuration:===== | ||
predict = True | predict = True | ||
totaltime = "1200s" | totaltime = "1200s" | ||
Use appropriate antenna configurations based on desired angular resolution | =====Use appropriate antenna configurations based on desired angular resolution===== | ||
repodir=os.getenv("CASAPATH").split(' ')[0] | repodir=os.getenv("CASAPATH").split(' ')[0] | ||
antennalist = repodir+"/data/alma/simmos/alma.out20.cfg" | antennalist = repodir+"/data/alma/simmos/alma.out20.cfg" | ||
Deconvolve the visibilities back into an image | =====Deconvolve the visibilities back into an image===== | ||
image = True | image = True | ||
vis = "$project.ms" | vis = "$project.ms" | ||
imsize = [192, 192] | imsize = [192, 192] | ||
Specify number of iteration of cleaning task with proper threshold and weighting | =====Specify number of iteration of cleaning task with proper threshold and weighting===== | ||
niter = 10000 | niter = 10000 | ||
threshold = "1e-7Jy" | threshold = "1e-7Jy" | ||
weighting = "natural" | weighting = "natural" | ||
We'd like to calculate a difference and fidelity image, and see some diagnostics: | =====We'd like to calculate a difference and fidelity image, and see some diagnostics:===== | ||
analyze = True | analyze = True | ||
And see the array but not the UV coverage: | =====And see the array but not the UV coverage:===== | ||
showarray = True | showarray = True | ||
showuv = False | showuv = False | ||
Plot both to the screen and the png files, and giving us lots of messages: | =====Plot both to the screen and the png files, and giving us lots of messages:===== | ||
graphics = "both" | graphics = "both" | ||
verbose = True | verbose = True | ||
overwrite = True | overwrite = True | ||
===To run the script, type ==== | |||
CASA<> execfile("Ppdisk.simdata2.txt") | CASA<> execfile("Ppdisk.simdata2.txt") |
Revision as of 22:34, 8 June 2010
↵ Simulating Observations in CASA
Protoplanetary disk
- 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 Wolf & D'Angelo 2005).
- Simdata2 version of script: File:Ppdisk.simdata2.txt
Explanation of the script
Set simdata2 as current task and reset all parameters
default("simdata2")
Specify sky model image
modelimage = "input50pc_672GHz.fits"
Image coordinate system can be verified
imhead("input50pc_672GHz.fits")
Image center can be identified
# ia.open("input50pc_672GHz.fits") # ia.shape() # [257L, 257L, 1L, 1L] # ia.toworld([128.5,128.5]) # {'numeric': array([ 4.71239120e+00, -4.01423802e-01, 1.00000000e+00, 6.72000001e+11])} # qa.formxxx("4.71239120rad",format='hms',prec=5) # '18:00:00.03052' # qa.formxxx("-0.401423802rad",format='dms',prec=5) # '-022.59.59.602743' # ia.done()
Brightness scale can be viewed with 'imstat' task
# imstat("input50pc_672GHz.fits") # ... # 'max': array([ 6.52469971e-05]), # ... # that's 0.0652 mJy/pixel.
Let's leave the brightness of the image as it is
inbright = "unchanged"
Let's call our project psim2
project = "psim2"
=====We'll leave the sky model the way it is -- simdata2 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===== modifymodel = False skymodel = "input50pc_672GHz.fits"
=====We need to decide where to point the telescope. The image is 2/3 arcsec in size, so we only need
one pointing. We could put that in a text file ourself, or let simdata2 create the ascii pointing file for us.===== setpointings = True direction = "J2000 18h00m00.031s -22d59m59.6s" mapsize = "0.76arcsec"
=====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 out20 configuration:=====
predict = True totaltime = "1200s"
Use appropriate antenna configurations based on desired angular resolution
repodir=os.getenv("CASAPATH").split(' ')[0] antennalist = repodir+"/data/alma/simmos/alma.out20.cfg"
Deconvolve the visibilities back into an image
image = True vis = "$project.ms" imsize = [192, 192]
Specify number of iteration of cleaning task with proper threshold and weighting
niter = 10000 threshold = "1e-7Jy" weighting = "natural"
We'd like to calculate a difference and fidelity image, and see some diagnostics:
analyze = True
And see the array but not the UV coverage:
showarray = True showuv = False
Plot both to the screen and the png files, and giving us lots of messages:
graphics = "both" verbose = True overwrite = True
To run the script, type =
CASA<> execfile("Ppdisk.simdata2.txt")
CASA<> go simdata2
- Output results:
Input: |
Predict: |
Image: |
Analyze: |