# Start by trying to create decent MEM VLA only image for Feathering
imgr:=imager(filename='orion.ms');
imgr.setdata(mode="none",spwid=[1,2],fieldid=[3:11]);
imgr.setimage(nx=300,ny=300,cellx="2arcsec",celly="2arcsec",stokes='I',
	mode='mfs',fieldid=7);
imgr.weight(type='uniform',mosaic=T);
imgr.setvp(dovp=T);
imgr.setmfcontrol(scaletype='SAULT',minpb=0.07,constpb=0.4);
imgr.mem(algorithm='mfentropy',
	displayprogress=T,
	niter=5,
	sigma='0.1Jy',
	targetflux='20Jy',
	constrainflux=F,
	image='maxen1.im',
	model='maxen1.model',
	residual='maxen1.resid')
#
# Make first mask based on maxen1.im
include 'interactivemask.g'
imask:=interactivemask(refimage='maxen1.im',
	mask='mask.image');
#
# Try maxen again using mask and higher target flux
imgr:=imager(filename='orion.ms');
imgr.setdata(mode="none",spwid=[1,2],fieldid=[3:11]);
imgr.setimage(nx=300,ny=300,cellx="2arcsec",celly="2arcsec",stokes='I',
	mode='mfs',fieldid=7);
imgr.weight(type='uniform',mosaic=T);
imgr.setvp(dovp=T);
imgr.setmfcontrol(scaletype='SAULT',minpb=0.07,constpb=0.4);
imgr.mem(algorithm='mfentropy',
	displayprogress=T,
	niter=10,
	sigma='0.05Jy',
	targetflux='120Jy',
	constrainflux=F,
	image='maxen2.im',
	model='maxen2.model',
	residual='maxen2.resid',
	mask='mask.image');#
#
# Make a new mask based on maxen2.im
include 'interactivemask.g'
imask:=interactivemask(refimage='maxen3.im',
	mask='mask.image');
imask.start()
#
# Try maxen again using mask and higher target flux
imgr:=imager(filename='orion.ms');
imgr.setdata(mode="none",spwid=[1,2],fieldid=[3:11]);
imgr.setimage(nx=300,ny=300,cellx="2arcsec",celly="2arcsec",stokes='I',
	mode='mfs',spwid=[1,2],fieldid=7);
imgr.setvp(dovp=T);
imgr.weight(type="briggs", rmode="norm",robust=-1,mosaic=T);
imgr.make('maxenjmf2.model');
imgr.mem(algorithm='mfentropy',
	displayprogress=T,
	niter=45,
	sigma='4mJy',
	targetflux='190Jy',
	constrainflux=F,
	image='maxen5.im',
	model='maxen5.model',
	residual='maxen5.resid',
	mask='mask.image');
#
# Try to reproduce orion_vlamem.im
imgr:=imager(filename='orion.ms');
imgr.setdata(mode="none",spwid=[1,2],fieldid=[3:11]);
imgr.setimage(nx=300,ny=300,cellx="2arcsec",celly="2arcsec",stokes='I',
        mode='mfs',spwid=[1,2],fieldid=7);
imgr.weight(type="briggs", rmode="norm", robust=-1,mosaic=T);
imgr.setvp(dovp=T);
imgr.setmfcontrol(stoplargenegatives=-1, cyclefactor=3, cyclespeedup=500);
imgr.make('maxenjm.model');
imgr.mem(algorithm='mfentropy',
         displayprogress=T,
         niter=23,
         sigma='0.004Jy',
          targetflux='180Jy',
          constrainflux=F,
          image='maxenjm.im',
          model='maxenjm.model',
          residual='maxenjm.resid',
          mask='orion.mask6');

#
# Next try VLA only multi-scale clean
imgr:=imager(filename='orion.ms');
imgr.setdata(mode="none",spwid=[1,2],fieldid=[3:11]);
imgr.setimage(nx=300,ny=300,cellx="2arcsec",celly="2arcsec",stokes='I',
	mode='mfs',spwid=[1,2],fieldid=7);
imgr.setvp(dovp=T);
imgr.weight(type="briggs", rmode="norm",robust=-1,mosaic=T);
imgr.setmfcontrol(stoplargenegatives=-1,cyclefactor=3,cyclespeedup=500);
imgr.setscales(scalemethod='uservector',uservector=[0,3,10,30]);
imgr.make('msc3.model');
imgr.clean(algorithm='mfmultiscale',
	displayprogress=T,
	niter=12000, 
	gain=0.2,
	threshold='0.02Jy',
	image='msc3.im',
	residual='msc3.resid',
	model='msc3.model',
	mask='orion.mask6');
##### Try with different resolutions
imgr:=imager(filename='orion.ms');
imgr.setdata(mode="none",spwid=[1,2],fieldid=[3:11]);
imgr.setimage(nx=300,ny=300,cellx="2arcsec",celly="2arcsec",stokes='I',
	mode='mfs',spwid=[1,2],fieldid=7);
imgr.setvp(dovp=T);
imgr.weight(type="briggs", rmode="norm",robust=-1,mosaic=T);
imgr.setmfcontrol(stoplargenegatives=-1,cyclefactor=3,cyclespeedup=500);
imgr.setscales(scalemethod='uservector',uservector=[0,4,12,36]);
imgr.make('msc5.model');
imgr.clean(algorithm='mfmultiscale',
	displayprogress=T,
	niter=8400, 
	gain=0.2,
	threshold='0.01Jy',
	image='msc5.im',
	residual='msc5.resid',
	model='msc5.model',
	mask='orion.mask6');
# 
# Image is very similar at 8400 or 10,000 iterations.
#
# Feathering using best MSC output
include 'imager.g'
myim:=imager('orion.ms');
myim.setimage();              # setimage is a dummy to get things initialized.
myim.setvp(dovp=T);
myim.feather(image='msc_feather.im', highres='msc3.im',
	lowres='orion.gbt2.im');
myim.done();

im:=image('msc_feather.im')
im.setrestoringbeam(major='135.0arcsec', minor='75arcsec', pa='10deg')
im.sbu('Jy/beam')
im.tofits('msc_feather.fits')
#
#
# Try now using orion.gbt2.im as input model for msc 
imv:=image('msc3.im');
csys:=imv.coordsys();
img:=image('orion.gbt2.im');
imgre:=img.regrid(outfile='orion.gbt2.imre',shape=imv.shape(),csys=csys);
imv.done();
img.done();
imgre.done();
#
include 'deconvolver.g'
decon:=deconvolver(dirtyname='orion.gbt2.imre',psfname='');
decon.makegaussian(gaussianimage='gbt_gaussian.im', 
	bmaj='85arcsec',bmin='85arcsec',bpa='0deg',normalize=F);
decon.open(dirtyname='orion.gbt2.imre',psfname='gbt_gaussian.im');
decon.setscales(scalemethod='uservector',uservector=([30,60]));
decon.clean(algorithm='msclean',model='orion.gbt2.mod',niter=200,
	gain=0.3,displayprogress=T,mask='orion.mask6');
decon.done();
#
imgr:=imager(filename='orion.ms');
imgr.setdata(mode="none",spwid=[1,2],fieldid=[3:11]);
imgr.setimage(nx=300,ny=300,cellx="2arcsec",celly="2arcsec",stokes='I',
	mode='mfs',spwid=[1,2],fieldid=7);
imgr.setvp(dovp=T);
imgr.weight(type="briggs", rmode="norm",robust=-1,mosaic=T);
imgr.setmfcontrol(stoplargenegatives=-1,cyclefactor=3,cyclespeedup=500);
imgr.setscales(scalemethod='uservector',uservector=[0,3,10,30]);
imgr.clean(algorithm='mfmultiscale',
	displayprogress=T,
	niter=10000, 
	gain=0.2,
	threshold='0.01Jy',
	image='msc_gbt2.im',
	residual='msc_gbt2.resid',
	model='orion.gbt2.mod',
	mask='orion.mask6');