This is a set of commands to take an observed, calibrated measurement
set and inject some models to test whether they are recoverable by
imaging.  

Delete false attempts or change base name if you go wrong and try again

Identify fully-calibrated parent data set to provide noise

Get the synthesis scripts

http://svnhq16.hq.eso.org/p16/pub/arc/software/casarelated/simulations/gensynthcube.py
http://svnhq16.hq.eso.org/p16/pub/arc/software/casarelated/simulations/

start casa

# set some variables to avoid having to set parameters multiple times
thisconfig = "alma;5.0arcsec"  # 5.0arcsec is ~the synthesized beam of the parent data set

mycell = '1.0arcsec'  # pixel size
mim = 216             # image size in pixels 
base = 'test'

####################################################################
# Make the parent  measurement set.
# This assumes that you can squash all to a single channel

# EXAMPLE using the calibrated data set from the NGC3256 casa guide (start of the imaging)
# split out come continuum-only channels
split(vis='ngc3256_line_target.ms',outputvis=base+'.ms',datacolumn='data',spw='2:31~120',width=90)

# OPTIONAL to remove signal from parent data set - not needed if pure noise anyway
clean(vis=base+'.ms',imagename=base+'orig.clean',interactive=T,imsize=mim,cell=mycell)
uvsub(vis=base+'.ms')
# END OF OPTIONAL BIT

####################################################################
# Make the model to inject

# Alternatively, you could use an observed image as long as it has
# clean components corresponding only and exactly to the desired
# model, in position and flux density, with no significant extra noise
# (i.e. well-masked during clean).

execfile('/home/amsr/scratch/ALMA/CASANOTES/TEST/gensynthcube.py')

# enter parameters below to correspond to test.ms and to the fake
#  sources you want to inject

myimdict = {
    'num_x_pixels'        : 50,       #[int] Number of x pixels in image
    'num_y_pixels'        : 50,       #[int] Number of y pixels in image
    'num_freq_pixels'     : 1,         #[int] Number of freq channels in image
    'freq_first_chan_ghz' : 99.6314,     #[float] Freq of first channel, in GHz
    'channel_width_ghz'   : 1.40625, #[float] Channel width, in GHz
    'restfreq_ghz'        :  99.6314,    #[float] rest frequency, in GHz
    'pixsize'             : 1.0,       #[float] pixel size of final cube in arcsec
    'RA'                  : 156.965,  #[float] right ascension of the image in degrees
    'DEC'                 : -43.905       #[float] declination of the image in degrees
    }
#------------------------------------------------------------
#Defining Gaussian properties
# - for N Gaussians, each array will contain N elements
# Enter the following into CASA

mygaussiandict = {
    'title': 'gaussian',
    'number'    : 4,                           #[int] 
    'x_center'         : np.array([10.0,20.,30.,40.]),    #[float]
    'y_center'         : np.array([20.,10.,40.,30.]),    #[float]
    'z_center'         : np.array([45.,45.,45.,45.]),           #[float]
    'x_fwhm'           : np.array([2.0,1.0,1.,2.]),         #[float]
    'y_fwhm'           : np.array([1.,2.,1.,2.]),         #[float]
    'z_fwhm'           : np.array([1.0,1.,1.,1.]),         #[float]
    'arot'             : np.array([0.0,45.,90.,135.]),         #[float]
    'integ_flux_first_plane' : np.array([0.1,0.01, 0.001, 0.01]),   #[float] set to zero if no continuum needed
    'spectral_index'    : np.array([0.0,0.0,0.,0.]),         #[float]
    'line_peak_flux'    : np.array([0.0,0.0,0.,0.])          #[float] set to zero if no line needed
    }
# Other models are possible, see gensynthcube-test.py

# Make the fake cube (I can only get it to work for a single channel)
gensynthcube(base+'cube', myimdict, mygaussiandict)


# Use sim_observe to produce a model ms
sim_observe(project=base+'',skymodel=base+'cube.im',graphics='none',antennalist=thisconfig)

# Clean the model, setting masks to select the model regions (this
# could also be done non-interactively by deriving a mask based on
# mygaussiandict inputs).

mms=thisconfig.replace(';','_') 
clean(vis=base+'/'+base+'.'+mms+'.ms',imagename=base+'cc.clean',interactive=T,imsize=mim,cell=mycell)

####################################################################
# Put it all together

# insert model into measurement set
ft(vis=base+'.ms',model=base+'cc.clean.model')

# add model visibilities to data
uvsub(vis=base+'.ms',reverse=T)

# Finally, clean data set with model injected
clean(vis=base+'.ms',imagename=base+'model.clean',interactive=T,imsize=mim,cell=mycell)