# 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='1310+323.ll');
mf.select(defpass='none',sbandav='usb');
mf.fill(msfile='1310+323.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 98apr16.ms');
shell ('cp -r 3c273.ms 98apr16.ms');

include 'ms.g';
mset:=ms(filename='98apr16.ms',readonly=F);
mset.concatenate(msfile='1310+323.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='98apr16.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('98apr16.ms');

include 'msplot.g';
# look at the phases of 1310+323 first
plot:=msplot(msfile='98apr16.ms',edit=T);  
# look at the amps of 1310+323 next
plot:=msplot(msfile='98apr16.ms',edit=T);  

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

# flag antenna 3 completely (very bad phases)
af:=autoflag('98apr16.ms'); 
af.setdata();
af.setselect(ant=3);    
af.run();
af.done();

# flag antenna 6 over 092801 to 104034 (bad phases)
af:=autoflag('98apr16.ms');
af.setdata();
af.setselect(ant=6,timerng="1998/04/16/09:28:01 1998/04/16/10:40:34");
af.run();
af.done();

# flag antenna 6 over 072700 to 084002 (bad phases)
af:=autoflag('98apr16.ms');
af.setdata();
af.setselect(ant=6,timerng="1998/04/16/07:27:00 1998/04/16/08:40:02");
af.run();
af.done();

# flag baseline 6-10 over 092801 to 113000 (bad phases)
af:=autoflag('98apr16.ms');
af.setdata();
af.setselect(baseline=[6,10],timerng="1998/04/16/09:28:01 1998/04/16/11:30:00");
af.run();
af.done();

# flag baseline 2-6 over 092801 to 113000 (bad phases)
af:=autoflag('98apr16.ms');
af.setdata();
af.setselect(baseline=[2,6],timerng="1998/04/16/09:28:01 1998/04/16/11:30:00");
af.run();
af.done();

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

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

# obtain gain solutions
include 'calibrater.g';

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

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

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

# 6. Transfer the Flux Density Scale

cal.fluxscale(tablein='98apr16.gcal',tableout='98apr16.fcal',reference='3C273',transfer=['1310+323'],refspwmap=[1,1]);
cal.plotcal(tablename='98apr16.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 [2]');
cal.setapply(type='G',t=0.0,table='98apr16.fcal');
cal.correct();

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

cal.done();

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

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

dc.summary('1310+323.split.ms');
dc.summary('ngc4826.split.ms');

# 9. Make an image of the calibrator
include 'imager.g';
imgr:=imager(filename='1310+323.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="1310+323.cln",image="1310+323.cm",residual="1310+323.resid");
imgr.done();
dc.view('1310+323.cm');
# looks good!

# 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();