Difference between revisions of "PPdisk simdata (CASA 3.1)"

From CASA Guides
Jump to: navigation, search
(Protoplanetary disk)
Line 6: Line 6:
  
 
*Simdata2 version of script: [[File:Ppdisk.simdata2.txt]]
 
*Simdata2 version of script: [[File:Ppdisk.simdata2.txt]]
 +
 +
---------------
 +
<br>
 +
*Explanation of the script:
 +
<br>
 +
 +
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  
 
*To run the script, type  
Line 14: Line 84:
  
  
 +
--------
 +
--------
 
*Output results:
 
*Output results:
  

Revision as of 18:16, 8 June 2010

Simulating Observations in CASA

Protoplanetary disk



  • 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:
Psim2.skymodel.png
Predict:
Psim2.predict.png
Image:
Psim2.image.png
Analyze:
Psim2.analysis.png