# set simdata as current task and reset all parameters default("simdata2") # use protoplanetary disk sky model image modelimage = "input50pc_672GHz.fits" # what is the image coordinate system? # imhead("input50pc_672GHz.fits") # in the logger, we see that this image had dummy (length=1) spectral and stokes # axes, and a single 32 GHz side "channel" centered at 672 GHz: # Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units # ---------------------------------------------------------------------------------------------------- # 0 0 Direction Right Ascension SIN 257 257 18:00:00.002 255.00 -3.110000e-03 arcsec # 1 0 Direction Declination SIN 257 1 -23.00.00.002 0.00 3.110000e-03 arcsec # 2 1 Stokes Stokes 1 1 I # 3 2 Spectral Frequency 1 1 6.72e+11 0.00 3.2000000029800e+10 Hz # where is the center of the image:? # 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() # what about the brightness scale, which can be changed independently of the WCS?: # imstat("input50pc_672GHz.fits") # ... # 'max': array([ 6.52469971e-05]), # ... # that's 0.0652 mJy/pixel. # so we'll 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" antennalist = "/opt/casa/data/alma/simmos/alma.out20.cfg" # and then deconvolve the visibilities back into an image image = True vis = "$project.ms" imsize = [192, 192] # cleaning relatively deepl 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 # plotting both to the screen and the png files, and giving us lots of messages: graphics = "both" verbose = True overwrite = True