########################################
## spectline_calib.g
##
## CALIBRATION OF SPECTRAL LINE DATACUBE
##
## JHIBBARD 5/14/02
## called by specline_doall.g
######################################### 
##
## REQUIRES "BASENAME", PHASECAL, PHASECAL FIELDID 
## TO BE DEFINED FROM WITHOUT

## User Specified Parameters

FLUXCAL			:= '1331+305';
FLUXCAL_FIELDID		:= 1;
BPASSCAL		:= PHASECAL;
BPASSCAL_FIELDID	:= PHASECAL_FIELDID;
BPASSCAL		:= FLUXCAL;
BPASSCAL_FIELDID	:= FLUXCAL_FIELDID;
SOURCE_FIELDID		:= 3;
                                   ## The following specifies the line
			           ## channels used to calculate the GAIN
				   ## ie, the range for "channel 0"	
BCHAN			:= 10;     # first good channel
ECHAN			:= 118;    # last good channel

REFANT			:= 24;     # reference antennae. Choose one near the
                                   # center of the array and that is present
                                   # in all calibrator scans.

SOLINT_GCAL		:= 300.0;   #averaging time for gain solutions
SOLINT_BCAL		:= 186400.0; #averaging time for bpass solutions
                                    #MADE LARGER THAN 1 DAY FOR
				    #OBSERVATIONS THAT CROSS MIDNIGHT


## Set up names for calibration files

GCAL_TABLENAME		:= spaste(BASENAME, '.gcal');
GCALOUT_TABLENAME	:= spaste(BASENAME, '.gcalout');
BCAL_TABLENAME		:= spaste(BASENAME, '.bcal');

######################################### 

print 'Starting Calibration';

## First, set the theoretical flux for the flux calibrator:

include 'imager.g';

print 'Setting Flux Sacale for ',FLUXCAL;
#imagr:=imager(spaste(BASENAME,'.ms'), compress=T);
# could the compress have been fouling up msplot??

imagr:=imager(spaste(BASENAME,'.ms'));
imagr.setjy(fieldid=FLUXCAL_FIELDID);
imagr.done();

#######################
## GAIN CALIBRATION
#######################
##
## Solve for time dependant gains for flux and phase calibrators.
## (equivalent to CALIB and CLCAL in AIPS)


include 'calibrater.g';
cal:=calibrater(spaste(BASENAME,'.ms'));

print 'Setting up inputs for Gain calibration';

# Restrict data used for Gain calibration to channels with
# a decent response. Usually, the inner 75% of the band is used
# for this so-called "pseudo-continuum" or "channel 0"

cal.setdata(mode='channel',
	    start=BCHAN,
            nchan=(ECHAN-BCHAN+1),
            msselect=spaste('FIELD_ID==',FLUXCAL_FIELDID,
                          ' OR FIELD_ID==',PHASECAL_FIELDID)); 

## Solve for the relative gain variations, putting results into
## GCAL_TABLENAME .
## (equivalent to CALIB in AIPS).

cal.setsolve(type='G',
             t=SOLINT_GCAL,
             refant=REFANT,
             table=GCAL_TABLENAME);

cal.state();

print 'Solving for relative gains';
cal.solve();

#######################
## BANDPASS CALIBRATION
#######################

# Note: now we use all channels

cal.setdata(msselect=spaste('FIELD_ID==',BPASSCAL_FIELDID)); 

## Apply the relative gain table (GCAL_TABLENAME)
## before solving for the B Jones matrix, so that the resulting 
## bandpass is normalized to 1. 
## (equivalent to BPASS in AIPS).

cal.reset();                  # undoes previous setsolve

cal.setapply(type='G',
             t=0.0,
             table=GCAL_TABLENAME);

cal.setsolve(type='B',
             t=SOLINT_BCAL,
             refant=REFANT,
             table=BCAL_TABLENAME);

cal.state();

print 'Solving for bandpass using calibrator', BPASSCAL;
cal.solve();

# cal.plotcal(tablename=BCAL_TABLENAME, multiplot=T);

## I don't like that plotting too much, so I use the script 
## specline_plotbp.g instead.

include 'specline_plotbp.g';

##########################
## SET ABSOLUTE FLUX SCALE
##########################

## Scale relative gains in GCAL_TABLENAME based on known and observed 
## flux of flux calibrator and observed flux of phase calibrator.
## Resulting table GCALOUT_TABLENAME has absolute gain calibration.
## (equivalent to GETJY in AIPS).

print 'Bootstrapping flux for ',PHASECAL;
cal.fluxscale(tablein=GCAL_TABLENAME, 
              tableout=GCALOUT_TABLENAME,
              reference=FLUXCAL, 
              transfer=[PHASECAL]);

#FYI, got a flux of 5.335+/-0.002 for 1130-148

cal.plotcal(tablename=GCALOUT_TABLENAME);

#######################
## CORRECT DATA
#######################
## Finally, apply the properly scaled gain table (GCALOUT_TABLENAME) 
## and the relative bandpass table, and write the corrected data.
## (equivalent to CLCAL in AIPS).

cal.setdata();    # selects all fields, channels

cal.reset();      # undoes previous setapply, setsolve

cal.setapply(type='G',
             t=0.0,
	     select=spaste('FIELD_ID==',PHASECAL_FIELDID),
             table=GCALOUT_TABLENAME);

cal.setapply(type='B',
             t=0.0,
	     select=spaste('FIELD_ID==',BPASSCAL_FIELDID),
             table=BCAL_TABLENAME);

cal.state();

print 'Writing bandpass and gain corrected data';
cal.correct();    # interpolates G solutions from PHASECAL to SOURCE
                  # (stored in GCALOUT_TABLENAME), interpolates B 
                  # solutions from BPASSCAL to SOURCE, applies to 
                  # raw data and writes into CORRECT column of ms
cal.done();

print 'Calibration done';

# You might want to run msplot, specifying "corrected" under the
# dataselection
include 'msplot.g';
mymsplot:=msplot(spaste(BASENAME,'.ms'), edit=T);

## THIS IS SCREWED UP!!! 
## BEFORE: I CANNOT SEE MY CALIBRATED DATA IN MSPLOT WHEN I RUN THIS SCRIPT! 
## NEW: REFUSES TO RUN.
## I HAVE TO GET OUT OF AIPS++
## AND BACK IN AGAIN. SOME PROBLEM WITH UPDATING TABLE ROWS.
## GM SAYS THIS IS DUE TO FLAGGER NOT PROPERLY CLOSING A TABLE

## THERE IS A POINT AT 0,0 ON ALL POINTS!