# set simdata as current task and reset all parameters
default("simdata") 
# use ppdisk sky model image  
modelimage         =  "ppdisk672_GHz_50pc.fits"

# what is the image coordinate system?
# imhead("ppdisk672_GHz_50pc.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     0.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

# in simdata, we'll us the coordinate system as it is:
ignorecoord        =  False

# what about the brightness scale, which can be changed independently of the 
# WCS?:

# imstat("ppdisk672_GHz_50pc.fits")
# ...
#  'max': array([  6.52469971e-05]),
# ...
# that's 0.0652 mJy/pixel.  

# so we'll leave the brightness of the image as it is also:
inbright           =  "unchanged"

# we want to observe over an 8GHz bandwidth 672-680 GHz:
startfreq          =  "672.0GHz"
chanwidth          =  "8.0GHz"


# we'll call the project "psim"
project            =  "psim" 

# and observe for 20 min, 10sec integrations
totaltime          =  "1200s" 
integration        =  "10s" 

# this parameter doesn't matter much since we're not going to do a mosaic:
scanlength         =  5

# we'll point at the source:
direction          =  ['J2000 18h00m00.0s -23d00m00s']

# the input image is only 0.8arcsec in size, so to do a single pointing, 
# we can make pointingspacing large and/or mosaicsize small:
pointingspacing    =  "10arcsec"
mosaicsize = "1arcsec"


# I'll look at a graphical layout of the setup:
checkinputs        =  "yes"

# on my machine this is where the configuration files are located:
#antennalist        =  "/Applications/CASA.app/Contents/data/alma/simmos/alma.out20.cfg"
antennalist = "/users/rindebet/casa/svn/data/alma/simmos/alma.out20.cfg"


# I chose out20 to have a synthesized beam of about XXX"

# let's add thermal noise:
noise_thermal      =  True
# and observe under marginal high-freq conditions:
noise_mode         =  "tsys-atm"
# with tsys-atm, we only set PWV and ground temp - simdata will construct 
# an atmospheric model using the ALMA site characteristics:
user_pwv           =  1.0
t_ground           =  269.0


# now I'll set up the output / clean:  
# I want to make a synthesis image with small enough pixels to resolve the beam
cell               =  "0.004arcsec" 
imsize             =  [192, 192]

# we want to clean relatively deeply:
niter              =  10000
threshold          =  "0.01mJy"   # mJy/bm
psfmode            =  "clark"
weighting          =  "natural"


# and see the output displayed graphically.
fidelity           =  False
display            =  True
verbose            =  True

# now type go :)