* STEP3 
*
* This is a collection of procedures used to do the polarization calibration
*  for VIPS data.
*
*  1. step3setup -- to initialize variables and to give the user a chance
*                 to set key variables before calling the other
*                 procedures.
*  2. step3    -- To do everything else.
*
* NB: All of the supporting utility procedures are at the beginning of
*     this file, while the main procedures (listed above) are at
*     the end.
*
* REVISION HISTORY
* ----------------
*  v04Aug2004, C. Fassnacht -- First working version.
*  v09Oct2004, G. Taylor -- Tweak to get 15 GHz going
*

***********************************************************************
*
* First define all of the utility procedures:
*  step3vars
*  find3cals
*  step3calin
*  step3cl
*  step3init
*  loadstep2
*  step3cal1
*  step3cal2
*  step3oqim
*  step3ccedt
*  step3lpcal
*  step3split
*
***************************************************

***************************************************
* step3vars: Define variables
***************************************************

proc step3vars
   scalar k
   scalar workdisk
   scalar oqcomps
   scalar mapcell
   scalar mapsize
   scalar calcount
   array oqcc(4,10)
   array oqcl(4,50)
   string*1  vband
   string*6  uvclass
   string*6  dclass
   string*6  tclass
   string*12 uvroot
   string*12 fcalname
   string*12 oqname
   string*12 tmpstrng
   string*12 rootname
   string*12 fullname
   string*6  classname
   string*6  tmpclass
   string*16 vipsdir
   string*16 cal3(10)
finish

***************************************************
* find3cals: Finds calibrators in the SPLAT file
***************************************************

proc find3cals
   task 'getthead'
   inext 'su'
   invers 0
   keyword 'num row'
   getthead
   print ''
   j = keyvalue(1)
   print 'find3cals:',j,' total sources in table'
   task 'tabget'
* NB: Column 2 (second argument in pixxy) is the source name
* NB: Column 4 (second argument in pixxy) is the calcode
   calcount 0
   for k=1 to j
      pixxy k,4,0
      tabget
      if keystrng <> ''
         then pixxy k,2,0
         tabget
         calcount = calcount + 1
         cal3(calcount) = keystrng
         print 'find3cals: Found calibrator',cal3(calcount)
      end
  end
finish

***************************************************
* step3calinit: Initializes variables for CALIB task
*          NB: Variables that are NOT set by this
*              procedure: docal, gainuse, solmode, 
*              solint, refant.
***************************************************

proc step3calin
   task 'calib'
   qual    -1
   selband  -1
   selfreq  -1
   freqid   -1
   timerang 0,0
   bchan 0; echan 0
   antenn   0
   dofit    0
   antuse   0
   subarray 0; docrt -1
   uvrange  0,0
   wtuv     0
   weightit 0
   dopol    -1
   doband   -1
   bpver    -1
   flagver  0
   smooth   0
   solsub   0
   solmin   0
   aparm    0
   soltyp '  '
   solcon   0
   snver    0
   minamper 0
   minphser 0
   cparm    0
   snver    0
   antwt 0; gainerr 0
   clro
*
* Put in parameters for a point-source model
*  as a default.
*
   clr2n
   cmodel ''; smodel 0
   ncomp 0; flux 0; nmaps 0
   cmethod 'dft'
*
* Save the inputs.
*
   tput calib
finish

***************************************************
* step3cl: Runs clcal
*  NB: snver, gainuse, refant, and gainver are NOT 
*   set in this procedure.  Neither is the input file.
***************************************************

proc step3cl
   task 'clcal'
   source ''
   soucode ''
   calsour ''
   qual -1
   calcode ''
   timer 0
   subarray 0
   antenna 0
   selband -1
   selfreq -1
   freqid -1
   opcode 'cali'
   interpol '2pt'
   cutoff 0
   samptype ''
   bparm 0
   icut 0.1
   doblank 0
   dobtween 0
   smotype ''
   inver 0
*
   go; wait
finish

***************************************************
* step3init: Initializes parameters
***************************************************

proc step3init
   if vband = 'u' 
      then dclass 'dclnu'
      uvclass 'splat'
      tclass 'tbu1'
      oqname 'oq208-u'
      oqcomps 3
      mapsize 512
      mapcell 0.000125
      oqcc  0.001, -7.50E-04, -1.13E-03, 0.002
      oqcc(1,2) = 2.50E-04
      oqcc(2,2) = -2.75E-03
      oqcc(3,2) =-1.88E-03
      oqcc(4,2) = -8.75E-04
      oqcc(1,3) = -4.25E-03
      oqcc(2,3) = -5.88E-03
      oqcc(3,3) = -6.50E-03
      oqcc(4,3) = -2.50E-03
      oqcl  248 251 265 273 254 235 271 250 290 210 308 237
   end
   if vband = 'c'
      then dclass 'dclnc'
      uvclass 'splat'
      tclass 'tbc1'
      oqname 'oq208-c'
      oqcomps 2
      oqcc 0.00325,-6.50E-03,-3.00E-03,0.0065,-3.50E-03,-9.25E-03,-0.01025,0.0015
      oqcl 115,103,140,155,142,92,169,135
      mapsize 256
      mapcell 0.00025
   end
*
* Load the proper uv file
*
   indi workdisk
   outdi workdisk
   in2di workdisk
   inname uvroot
   inclass uvclass
   inseq 1
*
* Find the calibrators
*
   find3cals
finish

***************************************************
* loadstep2: Loads in images produced by step 2
***************************************************

proc loadstep2
*
* Initialize imlod parameters
*
   task 'imlod'
   outna ''
   outcla dclass
   outseq 1
   outdi outdi
   ncount 1
   dotable 1
   nfiles 0
   nmaps 0
   error 1
*
* Run imlod for each of the calibrators
*
   for k=1 to calcount
      print 'loadstep2: Loading ',cal3(k)
      infile vipsdir!!':'!!cal3(k)!!'.cmp_n'; go; wait
   end
finish

***************************************************
* step3cal1: Runs calib using the maps made of
*   each calibrator by difmap in step 2.  After
*   running calib, run clcal to apply the calibration.
*   The selfcal is phase only.
***************************************************

proc step3cal1
*
* Initialize
*
   step3calin
*
* Set up some run-specific variables
*
   inname uvroot
   docal 2
   gainuse 1
   snver 1
   solmode 'p'
   aparm(7) 2
   solint 20/60
*
* Load the map for each calibrator and run calib
*
   print '!! Running phase only calibration'
   in2di indi
   in2class dclass
   in2seq 1
   outcl  in2clas
   for k=1 to calcount
      calsour cal3(k),''
      in2na cal3(k)
      print '!! Running CALIB on',cal3(k)
      go; wait
   end
*
* Apply the solutions by running step3cl
*
   snver 1
   gainver 1
   gainuse 2
   step3cl
finish

***************************************************
* step3cal2: Re-runs calib, this time with amplitude
*   plus phase selfcal.
*   After running calib, run clcal.
***************************************************

proc step3cal2
*
* Initialize
*
   step3calin
*
* Set up some run-specific variables
*
   inname uvroot
   docal 2
   gainuse 2
   snver 2
   solmode 'a&p'
   solint 4
*
* Load the map for each calibrator and run calib
*
   print '!! Running amp + phase calibration'
   calsour '3c279',''
   in2di indi
   in2class dclass
   in2seq 1
   outcl  in2clas
   for k=1 to calcount
      calsour cal3(k),''
      in2na cal3(k)
      print '!! Running CALIB on',cal3(k)
      go; wait
   end
*
* Apply the solutions by running step3cl
*
   snver 2
   gainver 2
   gainuse 3
   step3cl
finish

***************************************************
* step3oqim: Images OQ208
***************************************************

proc step3oqim
   task 'imagr'
   source 'oq208',''
   calcode ''
   timer 0
   docal 2
   gainuse 3
   dopol -1
   blver -1
   doband -1
   smooth 0
   stokes ''
   bchan 1; echan 0
   channel 0
   nchav 1
   bif 1; eif 0
   outnam oqname
   outdi indi
   outver 0
   clr2n
   cellsize mapcell, mapcell
   imsize mapsize mapsize
   nfield 1
   do3dimag -1
   fldsize mapsize,mapsize,0
   rashift 0; decshift 0
   uvtaper 0
   uvrange 0
   guard 0
   rotate 0
   zerosp 0
   uvwtfn 'na'
   uvsize 0
   robust 0
   uvbox 0
   uvbxfn 1
   xtype 5; ytype 5
   niter 2000
   allokay 0
   nboxes 0
   clbox 0
   boxfile ''
   oboxfile ''
   gain 0.1
   flux 0.0005
   minpatch 256
   bmaj 0; bmin 0; bpa 0
   overlap 0
   phat 0
   factor 0
   cmethod ''
   imagrprm 25,0; imagrprm(10)=1
   ngauss 0
   wgauss 0
   fgauss 0
   maxpix 0
   dotv 0
   go; wait
finish

***************************************************
* step3ccedt: runs ccedt on the newly created
*  OQ208 map
***************************************************

proc step3ccedt
   task 'ccedt'
   inna oqname
   inclass 'icl001'
   inseq 0
   inver 1
   outver 2
   bcount 1; ecount 0
   cutoff 0
   boxfile ''
   nboxes oqcomps
   clbox oqcl
   nccbox -1 * oqcomps
   ccbox oqcc
   cparm 0
   go; wait
finish

***************************************************
* step3lpcal: Runs lpcal
***************************************************

proc step3lpcal
   task 'lpcal'
   inna uvroot
   inclass uvclass
   inseq 1
   calsour 'oq208',''
   timer 0
   selband -1
   selfreq -1
   freqid -1
   bchan 1; echan 0
   bif 1; eif 0
   antennas 0
   uvrange 0
   subarray 0
   docal 2
   gainuse 3
   blver -1
   flagver 0
   doband -1
   in2na oqname
   in2class 'icl001'
   in2seq 0
   in2disk indi
   in2vers 2
   ncomp oqcomps, 0
   nmaps 1
   solint 4
   prtlev 0
   cparm 0
   go; wait
finish

***************************************************
* step3split: splits out all the data and writes
*  to disk
***************************************************

proc step3split
*
* Split out calibrated data
*
   task 'split'
   inna uvroot
   inclass uvclass
   inseq 1
   sour ''
   calcode ''
   timer 0
   stokes ''
   selband -1
   selfreq -1
   freqid -1
   bif 1; eif 0
   bchan 1; echan 0
   subarray 0
   docal 2
   gainuse 3
   dopol 1
   blver -1
   flagver 0
   doband -1
   smooth 0
   outclass tclass
   outseq 0
   douvcomp 1
   aparm 2,0
   go; wait
*
* Set up fittp to write split files out to disk
*
   task 'fittp'
   doall -1
   inclass tclass
   intype 'uv'
   outtape 0
   dostokes -1
   donewtab 1
   format 0
   blocking 10
   catno 0
   tput fittp
*
* Get the number of rows in the SU table
*
   inna uvroot
   inclass uvclass
   inseq 1
   task 'getthead'
   inext 'su'
   invers 0
   keyword 'num row'
   getthead
*
* Now loop through the SU table getting the source name for each row.
* For each calibrator found, write it out to disk
*
   print ''
   j = keyvalue(1)
   print 'Found ',j,' sources'
   task 'tabget'
   for k=1 to j
      inna uvroot; inclass uvclass
      pixxy k,2,0
      tabget
      tget fittp
      inna keystrng
      outfil = vipsdir!!':'!!keystrng!!'.'!!tclass;  go fittp; wait fittp
      print outfil
  end
finish

**********************************************************************
*
* step3setup:  Prompt the user to define appropriate variables.
*
***********************************************************************

proc step3setup

*
* First call the setup procedure to set the
*  general variables.
*
   if (vipsdir = '') 
      then vipsdir = 'vipsdata'
   end
*
* Now print summary of STEP3 variable values
*
   print ''
   print '*** Current values of STEP3 variables ***'
   print ''
   print 'workdisk = ',workdisk
   print 'uvroot = ',uvroot
   print 'vband  = ',vband
   print 'refant = ',refant
   print 'vipsdir= ',vipsdir
   print ''
   print ''
   print 'Examples:'
   print '  workdisk: 1       (value for indisk, outdisk, in2disk)'
   print '  uvroot: bf072a-c  (root name for uv files)'
   print '  vband:  u         (observing band: u or c)'
   print '  vipsdir: vipsdata (directory with data)'
   print '   NB: this environment variable must be set before starting'
   print '     aips'
   print ''
   print 'Please examine all variables (scroll up if needed)'
   print '*** If changes are desired, reset values now ***'
   print ''
*   print 'For a description of these variables, type monhelp'
   print ''

finish

***************************************************
*
* step3: Runs everything
*
***************************************************

proc step3
   step3init
   loadstep2
   step3cal1
   step3cal2
   step3oqim
   step3ccedt
   step3lpcal
   step3split
   print ''
   print 'Finished with step3.  Enjoy your polarized images'
finish

***************************************************
*
* This is the STEP3 pipeline for the VIPS program.
*
* To run this requires 2 sub-steps:
*   1. type step3setup to display variables for
*      step3.  Change the variables as needed.
*      Type step3setup as many times as required
*      to confirm settings.
*   2. type step3
*
* NOTE: Does not include EVPA calibration yet.
*
***************************************************