# set simdata as current task and reset all parameters
default("simdata") 

# 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            =  "psim"

# we'll leave the sky model the way it is 

ignorecoord        =  False
modelimage           =  "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. 

totaltime          =  "1200s"
direction          =  "J2000 18h00m00.031s -22d59m59.6s"
mosaicsize            =  "0.76arcsec"

# the default pointingspacing is fine - we'll only fit one pointing in the small mapsize

antennalist        =  "/opt/casa/data/alma/simmos/alma.out20.cfg"


# let's not do noise this time:

noise_thermal = False

# cleaning relatively deeply using the same size cell as the input
cell = '.00311arcsec'
imsize             =  [192, 192]
niter              =  10000
threshold          =  "1e-7Jy"
weighting          =  "natural"


# we'd like to calculate a difference and fidelity image, and see some diagnostics:

fidelity            =  True

# plotting both to the screen and the png files, and giving us lots of messages:
checkinputs = "yes"
verbose            =  True