#
# SCRIPTS FOR AIPS++ ALMA TST2.0
# (SINGLE-DISH INTERFEROMETER DATA COMBINATION METHODS)
#
# James Di Francesco, National Research Council of Canada  
#

# 1. "FEATHERING" TEST - using full channel maps not zero-moment maps

include 'imager.g';                     # - open the imaging tool

# set up primary beam file
imff:=imagefromfits(outfile='n4826_12mchan.beam',infile='NGC4826.primary_gain.fits'); # 
                                        # - convert the beam FITS file into an image
imff.summary();                         # - check the header axes, need to add extra axis
new_image:=imff.adddegaxes(stokes=T,outfile='n4826_12mchan.2.beam'); # 
                                        # - add a fourth header axes of polarization to 
                                        #   match other datasets
csys:=new_image.coordsys();             # - retrieve the new image's coordinate system
csys.reorder(order=[1,3,2]);            # - reorder the coordinate systems (RA&Dec=1)
new_image.setcoordsys(csys=csys);       # - put the coordinate system back in the image
new_image.done();                       # - close out the temporary tool
newim:=image('n4826_12mchan.2.beam');   # - open the new image
newim.summary();                        # - verify the new image axes

# set up single-dish data file
imff:=imagefromfits(outfile='n4826_12mchan.im',infile='NGC4826.12motf.chan.fits'); #
                                        # - convert the beam FITS file into an image
imff.summary();                         # - check the header axes, no extra axis needed!

# combine the data by "feathering"
img:=imager('n4826_both.ms'); 
img.setvp(dovp=T);
img.feather(image='n4826.fthr.im',highres='n4826_bima.im',lowres='n4826_12mchan.im',lowpsf='n4826_12mchan.2.beam');                   # - do the feathering, seemed to work OK
img.done();                             # - close the imaging tool

# examine the feathered data
include 'viewer.g';                     # - open the viewer tool
dv.gui;                                 # - fire up the data viewer

# data viewer statistics of BIMA-only data (n4826_bima.im):
#    Peak flux (all channels) = 1.889 Jy/beam at 12:56:45.388,+21:40:51.100
#    Flux Density on channel 7 = 3.084 Jy (assuming beam is present)
#    RMS (full plane, no mask) = 90.65 mJy/beam
#    Mean of RMS region = 1.65 mJy/beam

# data viewer statistics of feathered data (n4826.fthr.im):
#    Peak flux (all channels) = 1.980 Jy/beam at 12:56:45.388,+21:40:51.100
#    Flux Density on channel 7 = 10.05 Jy (assuming beam is present)
#    RMS (full plane, no mask) = 92.88 mJy/beam
#    Mean of RMS region = 3.55 mJy/beam

include 'imager.g';                     # - open the imaging tool
img:=image('n4826.fthr.im');            # - load in the feathered image
img.stats();                            # - determine bulk statistics
img.done();                             # - close the imaging tool 

# img.stats statistics of feathered data (n4826.fthr.im):
#    Peak flux (all channels) = 1.980 Jy/beam at 12:56:45.388,+21:40:51.100
#    Total flux Density = 132.4 Jy (assuming correct beam is present)
#    RMS (full plane, no mask) = 92.88 mJy/beam
#    Mean of RMS region = 3.55 mJy/beam

# Comments:
# This combination method seems adequate but I don't seem much difference between the 
# interferometer-only images and the combined images.  (This also seem true from the 
# combination shown on page 16 of the "detailed data information and test help" file.)
# It may be that there simply isn't much extended flux for this compact object that was
# resolved out, but the 12m flux is not terribly strong either.  Perhaps I am biased
# because I know that in low S/N situations the linear combinations performed by the
# BIMA SONG group may be safest route to pursue.

# The feathering option doesn't have any capability of "tweaking", namely altering some 
# parameters to vary (and improve) the combination.  It is what it is.  For example, is 
# there a way to scale the flux of one dataset so that it best matches the flux of the
# other in some region of overlap in the spatial frequency distributions of both datasets?
# (Note: this capability is actually described on pg. 101 of the Cookbook, but doesn't 
# say how the scale factor can be determined in AIPS++.)

# Also, I note the statistics I find from the provided BIMA-only image differ somewhat
# from those listed in the "detailed data information and test help" file...

# Notes:
# The above test was performed using the primary beam file provided with the test, and 
# not generating a voltage pattern as described in the Cookbook.  If time permits, I'll
# re-do the above with the vpmanager

# 2. "DECONVOLUTION" TEST - using full channel maps not zero-moment maps

#     a) produce a model by "directly converting" the SD image cube

include 'imager.g';                      # - open the imaging tool

# define the imaging parameters
im:=imager(filename='n4826_both.ms');    # - get the BIMA Measurement Set
im.setdata(fieldid=[1:7],spwid=[1,3]);   # - select the fields and spectral windows
im.setimage(nx=256,ny=256,cellx='1arcsec',celly='1arcsec',stokes='I',mode='channel',nchan=30,start=47,step=4,fieldid=1,spwid=[1:3]);   # - define the imaging parameters
im.setmfcontrol(scaletype="NONE",minpb=0.1); # - set no scaling, level at which PB is applied
im.setvp(dovp=T);                        # - do PB correction

# convert the single-dish image to a model
im.makemodelfromsd(sdimage='n4826_12mchan.im',modelimage='n4826_joint1',maskimage='n4826.mask');                                       # - make starting model image from SD image
im.clean(algorithm='mfclark',model='n4826_joint1',gain=0.2,niter=500); 
                                         # - do the CLEAN of the BIMA data with the model 

im.done();                               # - close the imaging tool

img:=image('n4826_joint1.restored');     # - grab the restored image
img.stats();                             # - provide image stats...

# img.stats statistics of joint deconvolved data (n4826_joint1.restored, full cube):
#    Peak flux (all channels) = 2.339 Jy/beam at 12:56:51.629,+21:40:04.090 (at the edge)
#    Total flux density = 167.5 Jy (assuming correct beam is present)
#    RMS (full, no mask) = 219.3 mJy/beam
#    Mean of RMS region = 4.50 mJy/beam

# data viewer statistics of joint deconvolved data (n4826_joint1.restored, inner box, ch.7):
#    Peak flux (channel 7) =  997 mJy/beam
#    Flux density on channel 7 = 3.205 Jy (assuming correct beam is present)
#    RMS (inner box, channel 7) = 141.8 mJy/beam
#    mean of rms region = 9.98 mJy/beam

#     b) produce a model by deconvolving the SD image cube

include 'imager.g';                      # - open the image tool

# regrid the SD image cube
ima:=image('n4826_joint1.restored');     # - grab the restored image cube from before
csys:=ima.coordsys();                    # - retrieve details of the mosaic image
csys.summary();                          # - get a summary of the coordinates
img:=image(infile='n4826_12mchan.im');   # - grab the SD image cube
imb:=img.regrid(outfile='n4826.12motf.chregrid.im',shape=[256,256,1,30],csys=csys); #
                                         # - regrid the SD image cube
ima.done();                              # - close tool
imb.done();                              # - close tool
img.done();                              # - close tool

# deconvolve the SD image cube
include 'deconvolver.g';                 # - open the deconvolver tool
decon:=deconvolver(dirtyname='n4826.12motf.chregrid.im',psfname=''); # 
                                         # - enter info for tool, note warning fixed below
decon.makegaussian(gaussianimage='n12m_gaussian.im',bmaj='77.1arcsec',bmin='77.1arcsec',bpa='134deg',normalize=F);                     # - make a Gaussian PSF
decon.open(dirtyname='n4826.12motf.chregrid.im',psfname='n12m_gaussian.im'); #
                                         # - enter info for tool
decon.setscales(scalemethod='uservector',uservector=[30,60]); #
                                         # - define scales to clean on
decon.clean(algorithm='msclean',model='n4826_joint2',niter=100,gain=0.3); #
                                         # - deconvolve the SD data to generate a model
decon.done();                            # - close the deconvolver tool

# deconvolve the interferometer data using the model
im:=imager(filename='n4826_both.ms');    # - open the imager tool with synthesis data
im.setdata(fieldid=[1:7],spwid=[1:3]);   # - select appropriate fields & spec. windows
im.setimage(nx=256,ny=256,cellx='1arcsec',celly='1arcsec',stokes='I',mode='channel',nchan=30,start=47,step=4,fieldid=1,spwid=[1:3]);   # - define image characteristics
im.setmfcontrol(scaletype="NONE",minpb=0.1); # 
                                         # - set multi-field processing parameters
im.setvp(dovp=T);                        # - do PB correction when imaging
im.clean(algorithm='mfclark',model='n4826_joint2',gain=0.2,niter=500,mask='n4826.mask'); #
                                         # - deconvolve BIMA data
im.done();                               # - close the imaging tool

img:=image(infile='n4826_joint2.restored'); #
                                         # - open the restored file for statistics
img.stats();                             # - get the statistics

# img.stats statistics of joint deconvolved data (n4826_joint2.restored, full cube):
#    Peak flux (all channels) = 2.621 Jy/beam at 12:56:51.272,+21:42:17.091 (at the edge)
#    Total flux density = 68.36 Jy (assuming correct beam is present)
#    RMS (full, no mask) = 223.7 mJy/beam
#    Mean of RMS region = 1.84 mJy/beam
# (matches with data viewer statistics of full cube)

# data viewer statistics of joint deconvolved data (n4826_joint2.restored, inner box, ch.7):
#    Peak flux (channel 7) =  992 mJy/beam
#    Flux density on channel 7 = 6.325 Jy (assuming correct beam is present)
#    RMS (inner box, channel 7) = 144.2 mJy/beam
#    mean of rms region = 24.5 mJy/beam

# Note: I'm not sure if the deconvolution succeeded, the command-line returned "F" after
# running the im.clean job.
# Note: restored data look OK (similar to before)

# Comments:
# I'm finding differences by factors of ~3 in the fluxes between the various combination 
# methods.  Am I doing something wrong, or is this expected?  (The latter doesn't appear
# expected from the examples in the "detailed description" file.)  I can't think of any
# other parameters to vary in the above, beyond gains and niters in the deconvolutions.
# The spectra for all 4 datasets, BIMA-only, Method 1, Method 2a, and Method 2b, at the
# same locations look all very similar to each other in shape and intensity, however,
# suggesting that the statistics (and grey scaling) have been skewed by outliers...