* STEP5
*
* This is a collection of procedures used to set the EVPA calibration
*  for VIPS data.
*
*  1. step5setup -- to initialize variables and to give the user a chance
*                 to set key variables before calling the other
*                 procedures.
*  2. step5    -- To do everything else. Namely
*                - Load Q and U images from step 4
*                - Make POLC and POLA images from them
*                - Calculate offsets needed to match to IF3 from 
*                    peak of polc on 3C279 C4.
*                - help user make EVPA correction based on IF3 or 4
*                - fix residual R/L phase differences on strong sources
*                - load all images from step 6
*                - make POLC, POLA and RM images
* 
* 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
* ----------------
*  v12Oct2004, G. Taylor -- Wrote procedure after Chris' Step3
*

***********************************************************************
*
* First define all of the utility procedures:
*  step5vars
*  find5cals
*  step5calin
*  step5cl
*  step5init
*  loadstep5
*
***************************************************

***************************************************
* step5vars: Define variables
***************************************************

proc step5vars
   scalar k
   scalar refval
   scalar workdisk
   scalar oqcomps
   scalar mapcell
   scalar mapsize
   scalar scount
   scalar calcount
   array c4(7)
   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)
   string*16 s3(10)
finish

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

proc find5cals
   task 'getthead'
   inext 'su'
   invers 0
   keyword 'num row'
   getthead
   print ''
   j = keyvalue(1)
   print 'find5cals:',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 'find5cals: Found calibrator',cal3(calcount)
      end
  end
finish


***************************************************
* find5all: Finds all sources in the SPLAT file
***************************************************

proc find5all
   task 'getthead'
   inext 'su'
   invers 0
   keyword 'num row'
   getthead
   print ''
   j = keyvalue(1)
   print 'find5all:',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
   scount 0
   for k=1 to j
      pixxy k,4,0
      tabget
      pixxy k,2,0
      tabget
      scount = scount + 1
      s3(scount) = keystrng
      print 'find5all: Found source',s3(scount)
  end
finish

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

proc step5calin
   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


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

proc step5cl
   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

***************************************************
* step5cal1: 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 step5cal1
*
* Initialize
*
   step5calin
*
* Set up some run-specific variables
*
   inname uvroot
   docal 2
   dopol 1
   gainuse 5
   snver 3
   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 'ICLN'
   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 step5cl
*
   snver 3
   gainver 5
   gainuse 6
   step5cl
finish


***************************************************
* step5init: Initializes parameters
***************************************************

proc step5init
   if vband = 'u' 
      then dclass '_U'
      c4 280 240 0
      uvclass 'splat'
      tclass 'tbu1'
      oqname 'oq208-u'
      oqcomps 3
      mapsize 512
      mapcell 0.000125
   end
   if vband = 'c'
      then dclass '_C'
      c4 139 121 0
      uvclass 'splat'
      tclass 'tbc1'
      oqname 'oq208-c'
      oqcomps 2
      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
*
   find5cals
*
* Find all sources
*
   find5all
finish

***************************************************
* loadstep5: Loads in images produced by step 4
***************************************************

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

***************************************************
* loadcal
***************************************************

proc loadcal
*
* Initialize imlod parameters
*
   task 'imlod'
   outna ''
   outcla ''
   outseq 0
   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


***************************************************
* loadall: Loads in all images 
***************************************************

proc loadall
*
* Initialize imlod parameters
*
   task 'imlod'
   outna ''
   outcla ''
   outseq 0
   outdi indi
   ncount 1
   dotable 1
   nfiles 0
   nmaps 0
   error 1
*
* Run imlod for each of the calibrators
*
   for k=1 to scount
      infile vipsdir!!':'!!s3(k)!!'.CMP_N'; go; wait
      for j=1 to 4
         print 'loadstep2: Loading ',s3(k)
         infile vipsdir!!':'!!s3(k)!!'.QCLN_'!!vband!!char(j); go; wait
         infile vipsdir!!':'!!s3(k)!!'.UCLN_'!!vband!!char(j); go; wait
      end
   end
finish

***************************************************
* loadpair: Loads in all images in 1+2 and 3+4 pairs
***************************************************

proc loadpair
*
* Initialize imlod parameters
*
   task 'imlod'
   outna ''
   outcla ''
   outseq 0
   outdi indi
   ncount 1
   dotable 1
   nfiles 0
   nmaps 0
   error 1
*
* Run imlod for each of the calibrators
*
   for k=1 to scount
      infile vipsdir!!':'!!s3(k)!!'.CMP_N'; go; wait
      for j=5 to 6
         print 'loadstep2: Loading ',s3(k)
         infile vipsdir!!':'!!s3(k)!!'.QCLN_'!!vband!!char(j); go; wait
         infile vipsdir!!':'!!s3(k)!!'.UCLN_'!!vband!!char(j); go; wait
      end
   end
finish

***************************************************
* step5polc: Makes polc images
***************************************************

proc step5polc
*
*
* Set up some run-specific variables
*
   blc 0; trc 0
   incl 'qcln'; in2cl 'ucln'
   opcode 'polc'
   outcl opcode
   aparm 0
   innam outn
*
* Load the map for each calibrator and run comb
*
   in2di indi
   for k=1 to calcount
      innam cal3(k)
      in2na cal3(k)      
      outn innam
      print '!! Making POLC and PPOL for ',cal3(k)
      for j=1 to 4
          inseq j; in2seq j
          blc 10 10; trc 40 40; tvlod; imstat
          bparm 0 0 1
          bparm(1)=pixstd
          incl 'ucln'
          blc 10 10; trc 40 40; tvlod; imstat
          bparm(2)=pixstd
          incl 'qcln'; blc 0; trc 0
          opcode 'polc'
          outcl opcode
          go comb; wait comb
          opcode 'pola'
          outcl opcode
          go comb; wait comb
      end
   end
*
*
finish

***************************************************
* step5pall: Makes polc images
***************************************************

proc step5pall
*
*
* Set up some run-specific variables
*
   blc 0; trc 0
   incl 'qcln'; in2cl 'ucln'
   opcode 'polc'
   outcl opcode
   aparm 0
   innam outn
*
* Load the map for each calibrator and run calib
*
   in2di indi
   for k=1 to scount
      innam s3(k)
      in2na s3(k)      
      outn innam
      print '!! Making POLC and PPOL for ',s3(k)
      for j=1 to 4
          inseq j; in2seq j
          blc 10 10; trc 40 40; tvlod; imstat
          bparm 0 0 1
          bparm(1)=pixstd
          incl 'ucln'
          blc 10 10; trc 40 40; tvlod; imstat
          bparm(2)=pixstd
          incl 'qcln'; blc 0; trc 0
          opcode 'polc'
          outcl opcode
          go comb; wait comb
          opcode 'pola'
          outcl opcode
          go comb; wait comb
      end
   end
*
*
finish

***************************************************
* step5uall: Makes polc images
***************************************************

proc step5uall
*
*
* Set up some run-specific variables
*
   blc 0; trc 0
   incl 'qcln'; in2cl 'ucln'
   opcode 'polc'
   outcl opcode
   aparm 0
   innam outn
*
* Load the map for each calibrator and run calib
*
   in2di indi
   for k=1 to scount
      innam s3(k)
      in2na s3(k)      
      outn innam
      print '!! Making POLC and PPOL for ',s3(k)
      for j=5 to 8
          inseq j; in2seq j
          blc 10 10; trc 40 40; tvlod; imstat
          bparm 0 0 1
          bparm(1)=pixstd
          incl 'ucln'
          blc 10 10; trc 40 40; tvlod; imstat
          bparm(2)=pixstd
          incl 'qcln'; blc 0; trc 0
          opcode 'polc'
          outcl opcode
          go comb; wait comb
          opcode 'pola'
          outcl opcode
          go comb; wait comb
      end
   end
*
*
finish

***************************************************
* step5rm: Makes rm image and spectral index
***************************************************

proc step5rm
*
*
* Set up some run-specific variables
*
   blc 0; trc 0
   incl 'pola';
   aparm 0
*
* Load the map for each calibrator and run calib
*
   in2di indi; outdi indi
   for k=1 to scount
      innam s3(k)
      in2na s3(k)      
      outn innam
      print '!! Making RM image for ',s3(k)
      for j=1 to 4
          incl 'pola n'
          inseq j; outseq 0
          outn inna
          outcl 'pola'
          rename
      end
      incl 'pola'
      for j=1 to 8
          inseq j; in2seq j
          axinc 1; axval j 0; axref 1; axtyp 'wavsq'
          axdef
      end
      task 'mcube'
      incl 'pola'
      inseq 1
      in2seq 8
      in3seq 1
      outn innam
      outcl 'cube'
      outseq 0
      doalign 1; axref 1; ax2ref 0; npoints 8
      go mcube; wait mcube
      task 'trans'
      incl 'cube'
      outn innam
      inseq 1
      outcl 'bcue'
      transcod '312'
      blc 0; trc 0
      go trans; wait trans
      vers 'gbt'
      task 'rm'
      incl 'bcue'
      inseq 1
      blc 0; trc 0
      outn innam
      outseq 20
      scalr1 0
      aparm 4 4 0 0 0 0 0 1 5 20
      bparm 4642 5045 14937 15317
      go rm; wait rm
      outseq 21
      scalr1 1
      go rm; wait rm
      vers ''
   end
   print '!! Now blanking on total intensity '
   task 'blank'
   incl 'rotmes'; inseq 20
   in2cl 'icln'; in2seq 1
   opcode 'in2c'; doalign 1
   dparm 0 0 100 0.002
   txinc 1; tyinc 1
   for k=1 to scount
      print '!! Blanking ',s3(k)
      innam s3(k)
      in2na s3(k)      
      outn innam; outcl 'rmblan'; outseq 0
      go blank; wait blank
   end
*
   blc 0; trc 0
   incl 'qcln'; in2cl 'ucln'
   opcode 'polc'
   outcl opcode
   aparm 0
   innam outn
*
* Load the map for each calibrator and run comb
*
   in2di indi
   incl 'icln'; in2cl incl
   opcode 'spix'
   outcl 'spix'
   for k=1 to scount
      innam s3(k)
      in2na s3(k)      
      outn innam
      inseq 1; in2seq 2
      print '!! Making SPIX for ',s3(k)
      blc 10 10; trc 40 40; tvlod; imstat
      bparm 0 0 1
      bparm(1)=pixstd
      aparm(9)=3*bparm(1)
      inseq 2
      blc 10 10; trc 40 40; tvlod; imstat
      bparm(2)=pixstd
      aparm(10)=3*bparm(2)
      inseq 1; blc 0; trc 0; outseq 0
      go comb; wait comb
   end
*
*
finish

***************************************************
* step5corr1: Correct all IFs to the value of IF 3
*   Measured at a fixed pixel on 3C279
***************************************************

proc step5corr1
*
* Set up some run-specific variables
*
   inname '3c279'
   incl 'pola'
*
* determine reference value
*
   inseq 3
   pixxy c4
   imval
   refval = pixval
*
* Load the map for each IF
*
   for j=1 to 4
       inseq j
       incl 'polc'
       pixxy c4
       imval
       incl 'pola'
       imval
       clcor(j) = 2*(refval - pixval)
   end
*
* Apply the corrections
   inna uvroot
   incl 'splat'
   inseq 1
   source ''
   stokes ''
   freqid -1
   bif 0; eif 0
   timer 0
   anten 0
   subar 0
   gainver 3
   gainuse 4
   opcode 'polr'
   infil ''
   go clcor; wait clcor
*
finish

***************************************************
* step5oqim: Images OQ208
***************************************************

proc step5oqim
   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
   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

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

proc step5ccedt
   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

***************************************************
* step5lpcal: Runs lpcal
***************************************************

proc step5lpcal
   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

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

proc step5split
*
* Split out calibrated data
*
   task 'split'
   inna uvroot
   inclass uvclass
   outdi indi
   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 0
   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

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

proc step6split
*
* Split out calibrated data
*
   task 'split'
   inna uvroot
   inclass uvclass
   outdi indi
   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 6
   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

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

proc step5setup

*
* First call the setup procedure to set the
*  general variables.
*
   if (vipsdir = '') 
      then vipsdir = 'vipsdata'
   end
*
* Now print summary of STEP5 variable values
*
   print ''
   print '*** Current values of STEP5 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

***************************************************
*
* fixphases: Corrects R-L phase differences in time
*
***************************************************

proc fixphase
   step5init
   loadcal
   step5calin
   step5cal1
   step6split
   print ''
   print 'Finished with step6.  Phases fixed'
finish


***************************************************
*
* step5: Runs everything
*
***************************************************

proc step5
   step5init
   loadstep5
   step5polc
   step5corr1
   step5split
   print ''
   print 'Finished with step5.  EVPA almos calibrated'
   print 'Final correction needed by comparison with'
   print 'VLA POLCA observations then CLCOR/SPLIT/FITTP'
finish

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