#JEH> Continuation of test log from aips++ data reduction
#JEH> of NGC 4038 D-array HI spectral line data.
#JEH> Log 3. 12/13-14/01
#
#JEH> First, we need to resolve why BP soln for Ant14 is so
#JEH> bad. (shape is good, but level is ~half of others)
#JEH> Full glish commands for reduction up to this point are
#JEH> in n4038D.g (appended at end of this log).

TEST 0 (TESTING_LOGS/n4038D.g):
Ran 'n4038D.g' with ant 14,10 flagged as in AIPS, using
both 3c286 and 1130-148 as bp cals.
In plotting bcal, ant14 shows up at amplitude of ~1.4, 
while all other antennae are around ~2.4
Header reads "Fields = [1,2]"

TEST 1 (n4038Da.g):
Ran 'n4038D.g' with ant 14,10 flagged as in AIPS, using
only 1130-148 as bp cal.
Exact same result. 

TEST 2 (n4038Db.g):
Same as test1, but flagged ant 14 for full run.
In plotting bcal, all antenna are at the
same level, but at ~2.4, not 1.0. ???

TEST 3 (n4038Dc.g):
As test2, but only using 1130-148 as phase cal.
Exact same result as TEST 2.

TEST 4 (n4038Dd.g):
delete ant14, 10 for the full run.
Use both 3c286 and 1130-148 as bp cals.
Exact same result as TEST 2.

#JEH> So is the bp not supposed to be normalized to 1???
#JEH> Emailed George & Athol.

#JEH> Note that amplitude of solutions in n4038D.gcalout
#JEH> is ~0.23, so 1/10 of what is plotted from bcal.
#JEH> Actually, in aips I self-cal'ed so I didn't have to
#JEH> delete Ant 14. Lets try selfcal latter. For now, 
#JEH> I guess we'll stick with the script that flags
#JEH> ant14 for the whole run, and uses both for bpass.
#JEH> Copied into n4038D.g.
#JEH> For now, lets see what kind of fluxes we get out 
#JEH> when imaging. So first, I have to apply the calibration:

cal:=calibrater('n4038D.ms');
cal.setdata(msselect='FIELD_ID IN [1:3]');
cal.setapply(type='G',t=0.0,table='n4038D.gcalout',select='FIELD_ID<=2');
#JEH> Why field_id==2???
cal.setapply(type='B',t=0.0,table='n4038D.bcal');
cal.correct();
cal.done();

#JEH> Lets first make a map of "ch 0" (ch 11-107)

include 'imager.g'
imagr:=imager('n4038D.ms');
imagr.setdata(fieldid=3, mode='channel', nchan=96, start=11, step=1,
                spwid=1);

imagr.setimage(nx=256, ny=256, cellx='10arcsec', celly='10arcsec', 
                 stokes='I', fieldid=3, mode='channel', nchan=1, 
                 start=11, step=96, spwid=1);
imagr.weight(robust=-1);
imagr.make('n4038Dch0.model');
imagr.clean(algorithm='clark', niter=3000, gain=0.1, threshold='0.0Jy',
               model='n4038Dch0.model', image='n4038Dch0.restored',
               residual='n4038Dch0.residual');
imagr.restore(model='n4038Dch0.model', 
                 image='n4038Dch0.restored2');
imagr.done();

im:=image('n4038Dch0.restored2');
im.view();   

#JEH> Looks good. Looks like AIPS map.
#JEH> Lets compare fluxes between AIPS map and AIPS++:
#JEH>              aips++        AIPS
#JEH> noise:      0.45 mJy/b     0.63 mJy/b
#JEH> n4038 peak: 18.8 mJy/b    122.6 mJy/b
#JEH> n4038 flux: 91.5 mJy      589.4 mJy
#JEH> sw peak:    11.5 mJy/b     66.3 mJy/b
#JEH> sw flux:    11.6 mJy       70.8 mJy
#JEH> se peak:    20.0 mJy/b    113.3 mJy/b
#JEH> se flux:    22.0 mJy      114.9 mJy
#JEH> east peak:   5.9 mJy/b     32.3 mJy/b
#JEH> east flux:   5.2 mJy       31.6 mJy
#JEH> 
#JEH> Something is amiss. Here is calibration script,
#JEH> if anyone wants to check it:

#=======================================================

####################################################
####################################################
## n4038D.g
####################################################
####################################################
# Version of n4038D script that flags as in aips.
# No; flags ant 14 for whole run. I'll not do then
# when I try selfcal.
#
#include 'ms.g'
#fitsfilename :='/export/data_1/jhibbard/DATA/n4038/AH582/AH582.960603.LINE';
#myms:=fitstoms(fitsfile=fitsfilename, msfile='n4038D.ms');
#myms.done();

#########################################
##
## Flagging
##
######################################### 

include 'flagger.g'
fg:=flagger(msfile="/export/data_1/jhibbard/aips++/n4038D.ms" );
fg.setchan(1:127);
fg.setflagmode('flag');
fg.filter(column='DATA', operation='range', comparison='Amplitude', range=['0.00001Jy','5Jy']);
#
# THIS DOESN'T WORK - SOMEONE NEEDS TO TELL ME HOW TO FLAG A FRIGGIN
# BASELINE!
#fg.query('ANTENNA1==8 and ANTENNA2==23');
#
fg.setantennas([4]);
fg.setids(fieldid=[1]);
fg.state();
fg.timerange('1996/06/03/00:57:30','1996/06/03/01:02:00');
#
fg.setantennas([18]);
fg.setids(fieldid=[3]);
fg.state();
fg.timerange('1996/06/03/00:59:00','1996/06/03/02:56:00');
#
fg.setantennas([10]);
fg.setids(fieldid=[1]);
fg.state();
fg.timerange('1996/06/03/00:56:00','1996/06/03/01:12:00');
#
fg.setantennas([14]);
fg.setids(fieldid=[1,2,3]);
fg.state();
#fg.timerange('1996/06/03/00:56:00','1996/06/03/02:41:30');
fg.timerange('1996/06/03/00:56:00','1996/06/03/04:23:30');
#
fg.done();
#

#include 'msplot.g'
#mymsplot:=msplot('n4038D.ms', edit=F);
#mymsplot:=msplot('n4038D.ms', edit=T, flagfile="n4038Dflag01.flag");
#

########################################
##
## Calibration
##
######################################### 

include 'imager.g';
im:=imager('n4038D.ms');
im.setjy(fieldid=1);
im.done();
#
include 'calibrater.g';
cal:=calibrater('n4038D.ms');
#cal.setdata(mode='channel',nchan=1,start=11,step=96,msselect='FIELD_ID<=2');
cal.setdata(msselect='FIELD_ID<=2');
cal.setsolve(type='G',t=300.0,preavg=300.0,refant=24,table='n4038D.gcal');
cal.solve();
cal.plotcal(tablename='n4038D.gcal');
cal.fluxscale(tablein='n4038D.gcal', tableout='n4038D.gcalout',
      reference='1331+305', transfer=['1130-148']);
cal.plotcal(tablename='n4038D.gcalout');
cal.done();

##IT IS REQUIRED TO CLOSE THE TOOL AND RECREATE A NEW ONE.
cal:=calibrater('n4038D.ms');
cal.setdata(mode='channel', nchan=127, start=1, step=1, msselect='FIELD_ID<=2');
cal.setapply(type='G',t=0.0,table='n4038D.gcalout');
cal.setsolve(type='B',t=86400.0,preavg=86400.0,refant=24,table='n4038D.bcal');
cal.solve();
cal.plotcal(tablename='n4038D.bcal');
cal.done();

##IT IS REQUIRED TO CLOSE THE TOOL AND RECREATE A NEW ONE.
cal:=calibrater('n4038D.ms');
cal.setdata(msselect='FIELD_ID IN [1:3]');
cal.setapply(type='G',t=0.0,table='n4038D.gcalout',select='FIELD_ID==2');
#JEH> Why field_id==2???
cal.setapply(type='B',t=0.0,table='n4038D.bcal');
cal.correct();
cal.done();

#########################################
##
## MAPPING
##
######################################### 

#JEH> Lets first make a map of "ch 0" (ch 11-107)
include 'imager.g'
imagr:=imager('n4038D.ms');
imagr.setdata(fieldid=3, mode='channel', nchan=96, start=11, step=1,
                spwid=1);
imagr.setimage(nx=256, ny=256, cellx='10arcsec', celly='10arcsec', 
                 stokes='I', fieldid=3, mode='channel', nchan=1, 
                 start=11, step=96, spwid=1);
imagr.weight(robust=-1);
imagr.make('n4038Dch0.model');
imagr.clean(algorithm='clark', niter=3000, gain=0.1, threshold='0.0Jy',
               model='n4038Dch0.model', image='n4038Dch0.restored',
               residual='n4038Dch0.residual');
imagr.restore(model='n4038Dch0.model', 
                 image='n4038Dch0.restored2');
imagr.done();
im:=image('n4038Dch0.restored2');
im.view();