Date - 2000.02.14
Tester - S.T. Myers (AOC)
Platform - Linux (kernow)
Version - Weekly

Note: bugs are denoted by lines with prefix:
>>>BUG
and queries are denoted by lines with prefix:
>>>QUERY

-------------------------------------------------------------------------------

GOAL: continue automapping construction

1. Repeat previous steps to get first clean model on data

#
# Run aips++ on 1608 test data for automapping emulation
#
# Read data into measurements set
#
include 'readfits.g'
#
# Read in functions to use
#
include 'makedirty.g'
include 'imagestats.g'
include 'cleanresid.g'
#
# Make big low-res dirty map and find peak
#
stat_rec := makedirty(prefix='1608',qual='.bg',cell_siz='0.1arcsec',imag_siz=1024)
#
# Check if there is something above 6 sigma
#
pk_flx := stat_rec.flux
rms_flx := stat_rec.rms
if ( pk_flx.value > 6.0*rms_flx.value ) 
   {
   # Start new imager and make small hi-res field and clean it
   imgr := imager('1608.ms')
   phdir := stat_rec.direction
   resid_rec := cleanresid(ref imgr,prefix='1608',qual='.sm',cell_siz='0.05arcsec',imag_siz=512,phase_dir=phdir)
   #
   pk_flx := resid_rec.flux
   rms_flx := resid_rec.rms
   # 
   # Iterate on this position while there is signal
   #
   while ( pk_flx.value > 6.0*rms_flx.value ) 
      {
      resid_rec := cleanresid(ref imgr,prefix='1608',qual='.sm',cell_siz='0.05arcsec',imag_siz=512,phase_dir=phdir)
      #
      pk_flx := resid_rec.flux
      rms_flx := resid_rec.rms
      }
   #
   # Done with this field
   #
   imgr.close()
   imgr.done()
}

   This completed successfully.

2. Creat selfcal function

#
# Function to do phase-only selfcal on data
# Works on existing imager object imgr, and
# makes calibrater object from existing clean model
# (else uses point source)
# Returns glish record of last residual map stats
#
function selfcalphase(ref imgr,prefix,qual)
{	
   # Make calibrater object
   ms_name := spaste(prefix,".ms")
   cl:=calibrater(ms_name)
   # Some local parameters
   pre_name := spaste(prefix,qual)
   cln_name := spaste(pre_name,".clean")
   # check if clean model already exists, else use a 1 Jy point source
   # put the model visibilities in the MODEL_DATA column
   if ( tableexists(cln_name) )
      {
      print 'Using existing clean model table ',cln_name
      imgr.ft(model=cln_name)
      } 
   else
      {
      print 'Using 1 Jy point source'
      cln_name := ""
      imgr.setjy(fluxdensity="1 Jy")
      }
   # Make calibration table name
   tab_name := spaste(pre_name,".phasecal")
   # Set selfcal parameters for phase-only gain solution, then solve
   cl.setsolve (type="T", t=10.0, dophase=T, table=tab_name, append=F);
   cl.solve()
   # Apply solution to data
   cl.setapply (type="T", t=10.0, table=tab_name, select="");
   cl.correct()
   # Done with calibrater
   cl.close()
   # now check through the calibration table to get some stats
   # return glish record of phase cal stats
   pha_rms := 0.0
   return [ rms=dq.quantity(pha_rms,'deg') ]
}

   Currently returns dummy rms until I figure out the table format to
   extract stats.

include 'selfcalphase.g'
imgr := imager('1608.ms')
phstat := selfcalphase(imgr,'1608','.sm')

   Output: errors
   Caught an exception! Event type=run exception=(/aips++/weekly/code/trial/implement/MeasurementComponents/DOcalibrater.cc : 445) Failed AlwaysAssertExit retval
   Caught exception: Error in select expression: Parse error at or near 'ð'
Scanned so far: select from $1 where 

   It is not clear to me what the problem is from this cryptic
   message!  I think it is in the cl.setapply (with the null select
   expression.  Try by removing select="" parameter. Should be able to
   just pick up where it crashed

cl:=calibrater('1608.ms')
tab_name := '1608.sm.phasecal'
cl.setapply (type="T", t=10.0, table=tab_name);
cl.correct()

   This seems to have worked.
>>>BUG: It should be able to handle a null expression for a parameter
        and treat it as if it wasnt there (eg. default).
   Bug() it. 

3. Now look at table 1608.sm.phasecal with GUI table module.
   Construct table tool.
   Try graphical browser .browse - unclear just what is in the FIT
   columns.  The GAIN cells contain the phase solutions, GAIN[1,1,1]
   is sufficient.

   Try colnames to return column names.  
>>>BUG: DOESNT WORK FROM GUI!
   Bug() it. Try 

mytable.colnames().

   Output: TIME TIME_EXTRA_PREC INTERVAL ANTENNA1 FEED1 FIELD_ID
   ARRAY_ID OBSERVATION_ID SCAN_NUMBER CORRELATOR_ID PULSAR_BIN
   PULSAR_GATE_ID  FREQ_GROUP FIELD_NAME SOURCE_NAME GAIN REF_ANT
   REF_FEED REF_RECEPTOR REF_FREQUENCY REF_DIRECTION CAL_DESC_ID
   CAL_HISTORY_ID TOTAL_SOLUTION_OK TOTAL_FIT TOTAL_FIT_WEIGHT
   SOLUTION_OK FIT FIT_WEIGHT 

   Try to grab gains using getcolslice: GAIN contains [2,2,2] cells

gains := mytable.getcolslice(columnname='GAIN',blc=[1,1,1],trc=[1,1,1])
gains::shape

   Output: [1 1 1 224]

   Problem.  Only 112 entries have valid data in 1,1,1.
   Should now be able to extract

cgain := gains[1,1,1,1:112]

   Now to change to phases in degrees.  
   There doesnt seem to be a built-in glish fuction to extract
   amp,phase from real,imag though it would seem to be an obvious
   thing to do!  I guess just use the arg(x) function

phases := arg( cgain )*180.0/pi
include 'mathematics.g'
pharms := stddev( phases )
print pharms

   Output: 12.8047733

   To automate this, we extract the SOLUTION_OK column which is
   boolean [2,N]

gains := mytable.getcolslice(columnname='GAIN',blc=[1,1,1],trc=[1,1,2])
gains::shape

   Output: [1 1 2 224]

goods := mytable.getcol(columnname='SOLUTION_OK')
goods::shape

   Output: [2 224]

   Use the mask to turn the good gains into a vector

ggain := gains[goods]
print len(ggain)

   Output: 188

phases := arg( ggain )*180.0/pi
pharms := stddev( phases )
print pharms

   Output: 15.3763501

   Get min,max phase correction also

print min( phases ), max( phases )

   Output: -48.0400245 59.9661395

   Incorporate this into selfcalphase.g

#
# Function to do phase-only selfcal on data
# Works on existing imager object imgr, and
# makes calibrater object from existing clean model
# (else uses point source)
# Returns glish record of last residual map stats
#
include 'mathematics.g'
function selfcalphase(ref imgr,prefix,qual)
{	
   # Make calibrater object
   ms_name := spaste(prefix,".ms")
   cl:=calibrater(ms_name)
   # Some local parameters
   pre_name := spaste(prefix,qual)
   cln_name := spaste(pre_name,".clean")
   # check if clean model already exists, else use a 1 Jy point source
   # put the model visibilities in the MODEL_DATA column
   if ( tableexists(cln_name) )
      {
      print 'Using existing clean model table ',cln_name
      imgr.ft(model=cln_name)
      } 
   else
      {
      print 'Using 1 Jy point source'
      cln_name := ""
      imgr.setjy(fluxdensity="1 Jy")
      }
   # Make calibration table name
   tab_name := spaste(pre_name,".phasecal")
   # Set selfcal parameters for phase-only gain solution, then solve
   cl.setsolve (type="T", t=10.0, dophase=T, table=tab_name, append=F);
   cl.solve()
   # Apply solution to data
   cl.setapply (type="T", t=10.0, table=tab_name);
   cl.correct()
   # Done with calibrater
   cl.close()
   # open the table
   gtable := table(tab_name)
   # now check through the calibration table to extract the good gains
   gains := gtable.getcolslice(columnname='GAIN',blc=[1,1,1],trc=[1,1,2])
   goods := gtable.getcol(columnname='SOLUTION_OK')
   # done with table
   gtable.done()
   # mask out good gain sols into vector, do stats
   ggain := gains[goods]
   ngain := len(ggain)
   phases := arg( ggain )*180.0/pi
   pharms := stddev( phases )
   phamin := min( phases )
   phamax := max( phases )
   # return glish record of phase cal stats
   return [ rms=dq.quantity(pharms,'deg'), min=dq.quantity(phamin,'deg'), max=dq.quantity(phamax,'deg'), number=ngain ]
}