# AIPS++ path used to reduce BIMA data of NGC 4826

# 1. Fill BIMA data from disk into an AIPS++ measurement set

# i) fill phase calibrator data
include 'mirfiller.g';
mf := mirfiller(mirfile='1230+123.ll');
mf.select(defpass='none',sbandav='usb');
mf.fill(msfile='1230+123.ms',verbose=T);

# ii) fill passband/flux calibrator data
mf := mirfiller(mirfile='3c273');
mf.select(defpass='none',sbandav='usb');
mf.fill(msfile='3c273.ms',verbose=T);

# iii) fill target data
mf := mirfiller(mirfile='ngc4826.ll');
mf.select(defpass='none',splwin=[5,6,7,8]);
mf.fill(msfile='ngc4826.ms',verbose=T);

mf.done();

#) iv) concatenate individual source files
shell ('rm -r 98apr22.ms');
shell ('cp -r 3c273.ms 98apr22.ms');

include 'ms.g';
mset:=ms(filename='98apr22.ms',readonly=F);
mset.concatenate(msfile='1230+123.ms',freqtol='10MHz');
mset.concatenate(msfile='ngc4826.ms',freqtol='10MHz');
mset.summary();

mset.done();

# 2. Fix the Flux Density of the Primary Flux Calibrator

# set the flux density of the flux calibrator 
# Note: 23 Jy from "private communication"
include 'imager.g';
imgr:=imager(filename='98apr22.ms',compress=F)
imgr.setjy(fieldid=1,fluxdensity=23.0);
imgr.setdata(fieldid=1);
imgr.plotvis(type='model');
imgr.done();

# 2. Examine and Flag MS data

include 'catalog.g';
dc.summary('98apr22.ms');

include 'msplot.g';
# look at the amps,phases of 1230+123 first
plot:=msplot(msfile='98apr22.ms',edit=T);  

# look at the amps,phases of 3c273 next
plot:=msplot(msfile='98apr22.ms',edit=T);  

# flag blocks of data identified from msplot
include 'autoflag.g';   

# flag antenna 9 over 100629 to 110000 (bad phases)
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(ant=9,timerng="1998/04/23/10:06:29 1998/04/23/11:00:00");
af.run();
af.done();

# flag antenna 10 over 084600 to 110000 (bad phases)
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(ant=[6,10],timerng="1998/04/23/08:46:00 1998/04/23/11:00:00");
af.run();
af.done();

# flag baselines 3-7 over 084600 to 110000 (bad phases)
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(baseline=[3,7],timerng="1998/04/23/08:46:00 1998/04/23/11:00:00");
af.run();
af.done();

# flag antenna 6 over 000001 to 071720 (bad phases)
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(ant=6,timerng="1998/04/23/00:00:01 1998/04/23/07:17:20");
af.run();
af.done();

# flag baseline 3-5 (terrible phases!)
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(baseline=[3,5]);
af.run();
af.done();

# flag baseline 5-9 (terrible phases!)
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(baseline=[5,9]);
af.run();
af.done();

# flag baseline 1-8 (bad phase over one interval)
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(baseline=[1,8],timerng="1998/04/23/08:05:32 1998/04/23/09:20:57");
af.run();
af.done();

# flag baseline 1-3 (bad phase over one interval)
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(baseline=[1,3],timerng="1998/04/23/08:05:32 1998/04/23/09:20:57");
af.run();
af.done();

# flag baseline 4-9 (bad phase over one interval)
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(baseline=[4,9],timerng="1998/04/23/05:22:13 1998/04/23/06:38:43");
af.run();
af.done();

# flag antenna 8 over some time interval 
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(ant=8,timerng="1998/04/23/04:36:00 1998/04/23/04:57:00");
af.run();
af.done();

# flag antenna 6 over some time interval 
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(ant=6,timerng="1998/04/23/08:46:00 1998/04/23/11:00:00");
af.run();
af.done();

# flag end channels: 1,2 and 63,64
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(chan=[1,2,63,64],spwid=[2,3,4,5],field=[3,4,5,6,7,8,9]);
af.run();
af.done();

# flag target data on antenna 7 between 060302 and 063844 (big jump in
# phase solution)
af:=autoflag('98apr22.ms');
af.setdata();
af.setselect(ant=7,timerng="1998/04/23/06:03:02 1998/04/23/06:38:44");
af.run();
af.done();

# 4. Solve for the Amp/Phase Gains from the Phase Calibrator

# obtain gain solutions
include 'calibrater.g';

cal:=calibrater(filename='98apr22.ms');
cal.reset();
cal.setdata(msselect='FIELD_ID IN [1,2] && SPECTRAL_WINDOW_ID IN [1]');
cal.setsolve(type='G',t=0,refant=4,table='98apr22.gcal');
cal.state();
cal.solve();

# examine gain solutions
cal.plotcal(tablename='98apr22.gcal',plottype='AMP',spwids=1,multiplot=T);
cal.plotcal(tablename='98apr22.gcal',plottype='PHASE',spwids=1,multiplot=T);

# 5. Solve for Bandpass Gains
# not needed for (these) BIMA data!

# 6. Transfer the Flux Density Scale

cal.fluxscale(tablein='98apr22.gcal',tableout='98apr22.fcal',reference='3C273',transfer=['1230+123'],refspwmap=[1]);
cal.plotcal(tablename='98apr22.fcal',plottype='AMP',multiplot=F); 

# 7. Apply the calibrations to the data
# first to the phase calibrator
cal.reset();
cal.setdata(msselect='FIELD_ID IN [2] && SPECTRAL_WINDOW_ID IN [1]');
cal.setapply(type='G',t=0.0,table='98apr22.fcal');
cal.correct();

# then to the target fields
cal.reset();
cal.setdata(msselect='FIELD_ID IN [3:9] && SPECTRAL_WINDOW_ID IN [2:5]');
cal.setapply(type='G',t=0.0,table='98apr22.fcal',spwmap=[1,1,1,1,1]);
cal.correct();

cal.done();

plot:=msplot(msfile='98apr22.ms',edit=T);

# 8. Split the calibrated data into individual files
mset:=ms(filename="98apr22.ms",readonly=F);
mset.split(outputms="1230+123.split.ms",fieldids=2,spwids=1,nchan=1,start=1,step=1,whichcol="CORRECTED_DATA");
mset.split(outputms="ngc4826.split.ms",fieldids=[3,4,5,6,7,8,9],spwids=[2,3,4,5],nchan=64,start=1,step=1,whichcol="CORRECTED_DATA");
mset.done();

dc.summary('1230+123.split.ms');
dc.summary('ngc4826.split.ms');

# 9. Make an image of the calibrator
include 'imager.g';
imgr:=imager(filename='1230+123.split.ms');
imgr.setdata(fieldid=[1],spwid=[1]);
imgr.setimage(nx=256,ny=256,cellx='1arcsec',celly='1arcsec',stokes='I');
imgr.clean(algorithm="clark",niter=100,gain=0.1,threshold='0.01Jy',model="1230+123.cln",image="1230+123.cm",residual="1230+123.resid");
imgr.done();
dc.view('1230+123.cm');
# looks alright, appearance of second peak to the west a bit worrisome!

# 10. Make an image of the target fields
imgr:=imager(filename='ngc4826.split.ms');
imgr.setdata(fieldid=[1:7],spwid=[1:3],nchan=120,start=47,step=1);
imgr.setimage(nx=256,ny=256,cellx='1arcsec',celly='1arcsec',stokes='I',mode='channel',nchan=30,start=47,step=4,fieldid=1,spwid=[1:3]);

imgr.weight(mode='natural',mosaic=T); # this step did not work
imgr.setvp(dovp=T);
imgr.setoptions(padding=1.0,ftmachine='mosaic');
imgr.clean(algorithm='mfclark',niter=300,gain=0.1,model='ngc4826.cln',residual='ngc4826.resid',image='ngc4826.cm');

# 11. Get the zeroth and first moments
include 'image.g';
im:=image(infile='ngc4826.cm');
im_mom0:=im.moments(outfile='ngc4826.mom0.cm',moments=0,axis=4,mask='indexin(4,[5:26])',includepix=[0.070,1000.0]);
im_mom1:=im.moments(outfile='ngc4826.mom1.cm',moments=1,axis=4,mask='indexin(4,[5:26])',includepix=[0.28,1000.0]);
im_mom0.done();
im_mom1.done();
im.done();

# 12. Check out the moment maps!
include 'viewer.g';
dv.gui();