PPdisk simdata (CASA 3.1): Difference between revisions

From CASA Guides
Jump to navigationJump to search
No edit summary
 
(11 intermediate revisions by 3 users not shown)
Line 1: Line 1:
{{Simulations Intro}}
{{Simulations Intro 3.1}}
[[Category: Simulations]]
[[Category: Simulations]]


Line 5: Line 5:
*[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]).
*[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]).


*Simdata2 version of script: [[File:Ppdisk.simdata2.txt]]
*Simdata version for CASA 3.1
 
__NOTOC__
---------------
<br>
*Explanation of the script:
<br>
<br>
====Explanation of the script====


Set simdata2 as current task and reset all parameters
=====Set simdata as current task and reset all parameters=====
   default("simdata2")   
   default("simdata")   
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 29: 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 35: 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: simdata 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 simdata create the ascii pointing file for us.=====
  one pointing. We could put that in a text file ourself, or let simdata2 create the ascii
  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: let's do a 20 min snapshot observation using out20 configuration:=====
We do want to calculate visibilities in a measurement set:
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 with lots of messages:=====
   graphics          =  "both"
   graphics          =  "both"
   verbose            =  True
   verbose            =  True
   overwrite = True
   overwrite = True


*To run the script, type
===To run the script ===


   CASA<> execfile("Ppdisk.simdata2.txt")
   CASA<> execfile("Ppdisk.simdata.txt")


   CASA<> go simdata2
   CASA<> go simdata


 
<br>
--------
--------
*Output results:
*Output results:



Latest revision as of 18:28, 28 April 2011

Simulating Observations in CASA 3.1

Protoplanetary disk

  • Simdata version for CASA 3.1


Explanation of the script

Set simdata as current task and reset all parameters
 default("simdata")   
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: simdata 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 simdata 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 with lots of messages:
 graphics           =  "both"
 verbose            =  True
 overwrite = True

To run the script

 CASA<> execfile("Ppdisk.simdata.txt")
 CASA<> go simdata


  • Output results:


Input:
Psim2.skymodel.png
Predict:
Psim2.predict.png
Image:
Psim2.image.png
Analyze:
Psim2.analysis.png