*
* CALAX - set of procedures to semi-automatically reduce polarization
*         calibration data for the VLA for 3C48 datasets.
*
* Version STM-07-02-01 for running on SANDROCK USER 1553 Disks 1 and 2
* THIS ASSUMES THE TST VERSION OF AIPS! NOT COMPATIBLE WITH NEW OR EARLIER.
*
* NOTE: The session-specific variables are defined in procedure SETIT.
* This script should be useable for any run if these are changed
* appropriately.
*
* Steps to run, with data in SLOT I:
*
* 1. SETIT -- define the session specific variables (IMPORTANT!!!!)
*
*    Before starting be sure to quack off the first 10 sec (to zap the
*    pointing scans and bad data at start), and examine the data in
*    TVFLG for obvious problems (use AMPL DIFF and PHA DIFF with 3.3sec
*    records).
*
*    If you read in the data with 1min default CL entries:
*    Delete default CL 1 and run 
*    REINDXR(I) -- INDXR with CPARM(3) 10/60 to make CL 1 with 10sec entries
*    Also after UVCOP, use this.
*
* 2. CALDATA(I) -- gain calibration, using a Gaussian model, produces CL 2 and
*    CL 3 which is the one to be applied in following steps.  You can
*    examine the output in file *.CALIB to check the process.
* 4. LISTCAL(I) -- examine calibrators with LISTR OPTYP='MATX'
*
* Remember, if you need to redo things, proc RESETIT(I) will reset
*
Proc Setit
* Procedure to define global pseudoverbs, arrays and strings
* as well as to define the calibrator names and number for later use
*
  string*16 abscal, absnam, srccal
  string*12 absmod, predir
  string*8 abstag
  string*6 iclass, bclass
  string*2 prefix
  string*1 fqband
  scalar abcell, absref, numbox, calimode, sjparm, clniter 
  scalar doifavg, moddis, mapdis
  scalar gflux, gbmaj, gbmin, gbpa
  array absflx(2), absant(50), notant(50), absuvr(2), useuvr(2)
  array absbox(4), rmsbox(4)
*  
* Set up some class suffixes for the appropriate version of AIPS
*
  type 'Setting up for AIPS TST'
  iclass 'ICL001'; bclass 'IBM001'; 
*
* Calibrate weights (DOCALIB=2)
*
  calimode 2;
*  
* NOTE: These values are specific to this session
*
  type 'Using setup values for A-config X-Band 3C48 data'
*
  abstag '-22DEC00'
  predir 'PWD:'; type 'PREDIR =',predir
  prefix 'AX'; type 'PREFIX =',prefix
  fqband 'X'; type 'FQBAND =',fqband
  abscal '0137+331'; type 'ABSCAL =',abscal; 
  absnam '3C48'; type 'ABSNAM =',absnam
  absmod ''; type 'ABSMOD =',absmod
  moddis 0; mapdis 2; doifavg 1;
* setjy calc 1999.2
  absflx 0, 0; sjparm 0; 
* Gaussian model parameters
  gflux 3.1543; gbmaj 0.316; gbmin 0.0964; gbpa 19.1
* open recommended 0-40 klam to 0-60 klam and stations out to 24 (no model)
*  absuvr 0 60; type 'ABSUVR =',absuvr
*  absant 2 3 6 7 8 10 12 14 15 16 21 22 23 24 28 0
*  notant -2 -3 -6 -7 -8 -10 -12 -14 -15 -16 -21 -22 -23 -24 -28 0
* use 0-300 klam for Gaussian model for amplitudes
  absuvr 0 300; type 'ABSUVR =',absuvr
  absant 0
  type 'ABSANT =',absant
  notant 0
* use 0-400 klam for Gaussian model for phases
  useuvr 0 400; type 'USEUVR =',useuvr
*
  absref 8; type 'ABSREF =',absref
  numbox 0; type 'NUMBOX =',numbox
  absbox 0, 0, 0, 0; type 'ABSBOX =',absbox
  rmsbox 32, 32, 96, 96; type 'RMSBOX =',rmsbox
  abcell 0.05; type 'ABCELL =',abcell
  srccal ' '
  clniter 300
*
Finish
*
*
Proc Resetit(i)
* Reset dataset I for redo-ing of all procedures
*
  getn i
  task 'setjy'
  optyp 'rejy'
  source abscal,''
  bif 1; eif 2; zerosp 0
  freqid 1
  go setjy; wait setjy
*
  inext 'sn';inver -1;extd
*
  inext 'cl';inver 0; extd
*
  imh
  type 'There should now be no SN tables and 1 CL table'
  type 'May have to repeat EXTD if CL tables remain'
*
Finish
*
*
Proc Reindxr(i)
* Create new CL table with 10 sec entries
* Run after deleting CL 1 or running UVCOP
*
  getn i
  task 'indxr'
  infile ''; prtlev 0
  cparm 0, 0, 10/60, 0
  go indxr; wait indxr
  imh
  type 'There should now be 1 CL table with 10sec entries'
*
Finish
*
*
Proc Mapem
*
* Procedure to run IMAGR with stokes I
*
  tget imagr
  if doifavg > 0 then outn sour(1);
    outn sour(1)
    bif 1; eif 2
    outseq 1
    clr2n; in2di indi; outdi mapdis
    niter clniter
    stokes 'i'; outcl iclass
    wait imagr; go imagr
  else outn sour(1);
    for j=1 to 2;
      bif j; eif j
      outseq (j)
      clr2n; in2di indi; outdi mapdis
      niter clniter
      stokes 'i'; outcl iclass
      wait imagr; go imagr
    end
  end
  wait imagr
  type 'finished MAPEM'
  Finish
*
*
Proc caldata(i)
*
* Specialized procedure to calibrate data, possibly with ICLN model.  
* First do short (typically 10s) phase solutions followed by longer
* (~1min) full gain solutions.
*
  clrtemp
  prtas ''; clrmsg
*
  getn i
  task 'setjy'
  optyp ''
  aparm 0
  if absflx(1) = 0.0 then optyp 'calc'; aparm(2) sjparm; end
  source abscal,''
  bif 1; eif 1; zerosp absflx(1) 0
  freqid 1
  go setjy; wait setjy
  bif 2; eif 2; zerosp absflx(2) 0
  freqid 1
  go setjy; wait setjy
*
* First pass - phase-only
*
  task 'calib'
  if absmod <> '' then in2n absmod;
* Use CC 1 model attached to ABSMOD on Disk 1
    type 'Phase calibration using Clean model - all data'
    in2cl iclass; in2seq 1; in2d moddis; cmodel 'COMP'; cmethod 'dft'
    invers 1; flux 0; ncomp 0; nmaps 1; smodel 0;
    uvr 0 0; anten 0; basel 0
  else in2n '';
    if gflux > 0.0 then in2cl ''; 
      type 'Phase calibration using Gaussian model - restricted data'
      in2seq 0; in2d 0; cmodel ''; cmethod 'dft'
      invers 0; flux 0; ncomp 0; nmaps 0; 
      smodel gflux, 0, 0, 1, gbmaj, gbmin, gbpa
      uvr useuvr; anten 0; basel 0
    else in2cl ''; 
      type 'Phase calibration using Point model - restricted data'
      in2seq 0; in2d 0; cmodel ''; cmethod ''
      invers 0; flux 0; ncomp 0; nmaps 0; smodel 0;
      uvr useuvr; anten 0; basel 0
    end
  end
  calsour abscal,''
  flagver 0; gainuse 0; docalib calimode;
  refant absref
  aparm 4 0
* 10 second solutions
  solmode 'p'; solint 10/60
  cparm 0 0 10 10 0
  snver 1
  go calib; wait calib
*
* Apply each source to itself
*
  task 'clcal'
  type 'Applying Self Calibration for ',abscal
  source abscal,''
  calsour source
  opcode ' '
  interp 'simp'
  gainver 1
  gainuse 2
  go clcal; wait clcal
*
* Second pass - full gain (amplitude) solutions on longer timescales
*
  task 'calib'
  if absmod <> '' then in2n absmod;
* Use CC 1 model attached to ABSMOD on Disk 1
    type 'Amplitude calibration using Clean model - all data'
    in2cl iclass; in2seq 1; in2d moddis; cmodel 'COMP'; cmethod 'dft'
    invers 1; flux 0; ncomp 0; nmaps 1; smodel 0;
    uvr 0 0; anten 0; basel 0
  else in2n '';
    if gflux > 0.0 then in2cl ''; 
      type 'Amplitude calibration using Gaussian model - restricted data'
      in2seq 0; in2d 0; cmodel ''; cmethod 'dft'
      invers 0; flux 0; ncomp 0; nmaps 0; 
      smodel gflux, 0, 0, 1, gbmaj, gbmin, gbpa
      uvr absuvr; anten absant; basel anten
    else in2cl ''; 
      type 'Amplitude calibration using Point model - restricted data'
      in2seq 0; in2d 0; cmodel ''; cmethod ''
      invers 0; flux 0; ncomp 0; nmaps 0; smodel 0;
      uvr absuvr; anten absant; basel anten
    end
  end
  calsour abscal,''
  flagver 0; gainuse 0; docalib calimode
  refant absref
  aparm 4 0
* 30s gain solutions
  solmode 'a&p'; solint 0.5
  cparm 0 0 10 10 0
  snver 2
  go calib; wait calib
*
* Apply each source to itself
*
  task 'clcal'
  source abscal,''
  if absmod <> '' then calsour source;
    type 'Applying Self Calibration for ',abscal
    opcode ' '
    anten 0
    interp 'simp'
    gainver 2
    gainuse 3
    go clcal; wait clcal
  else calsour source;
    type 'Applying Self Calibration for ',abscal
    opcode ' '
    anten absant
    interp 'simp'
    gainver 2
    gainuse 3
    go clcal; wait clcal
    anten 0
  end
*
  prtas ''
  outprint predir!!prefix!!abstag!!'.CALIB'
  docrt -1
  prtmsg
  docrt 1
*
  type ' Completed CALDATA '
*
Finish
* 
*
Proc listcal(i)
*
* use LISTR 'GAIN' to output amplitude solutions
*
  task 'listr'
  getn i
  sources ''
  optyp 'gain'
  docalib 1; gainuse 0; dopol -1
  docrt -1
  antenn 0; basel 0
*
  bif 1; eif 1; stokes 'r'
  outpr predir!!prefix!!abstag!!'-1R.AMPL'
  inext 'sn';inver 2
  dparm 0 0 5 0
  go listr; wait listr
  inext 'cl';inver 0
  dparm 1 0 4 0
  go listr; wait listr
*
  bif 1; eif 1; stokes 'l'
  outpr predir!!prefix!!abstag!!'-1L.AMPL'
  inext 'sn';inver 2
  dparm 0 0 5 0
  go listr; wait listr
  inext 'cl';inver 0
  dparm 1 0 4 0
  go listr; wait listr
*
  bif 2; eif 2; stokes 'r'
  outpr predir!!prefix!!abstag!!'-2R.AMPL'
  inext 'sn';inver 2
  dparm 0 0 5 0
  go listr; wait listr
  inext 'cl';inver 0
  dparm 1 0 4 0
  go listr; wait listr
*
  bif 2; eif 2; stokes 'l'
  outpr predir!!prefix!!abstag!!'-2L.AMPL'
  inext 'sn';inver 2
  dparm 0 0 5 0
  go listr; wait listr
  inext 'cl';inver 0
  dparm 1 0 4 0
  go listr; wait listr
*
  docrt 1; outpr ''; inext ''; inver 0  
*
* use LISTR 'MATX' to examine gain calibration
*
  task 'listr'
  getn i
  sources ''
  stokes 'half'
  optyp 'matx'
  docalib 1; gainuse 0; dopol -1
  dparm 5,0
  docrt 1
*
  bif 1; eif 1; go listr; wait listr
*
  bif 2; eif 2; go listr; wait listr
*
  bif 1; eif 2; antenn 0; basel 0
*
  type ' Completed LISTCAL '
*
Finish
*
* 
Proc outrms 
*
* short proc to use IMSTAT for current image to output the off-source 
* statistics for I,Q,U images in IF 1 and 2 for the current source name 
* SOURCE(1).
*
* Expects BLC,TRC for off-source window stored in RMSBOX
  clrtemp
  inna source(1)
*
* IF 1 data, INSEQ=1
*
  if doifavg > 0 then indisk mapdis;
*
    type 'IF AVG On-source image statistics for source '!!inna
    blc absbox(1), absbox(2), 0; trc absbox(3), absbox(4), 0
    incl iclass; inseq 1; imstat; 
    type inna!!' IF=AVG I Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy '
    type inna!!' IF=AVG I Max: at ('!!CHAR(pixxy(1))!!','!!CHAR(pixxy(2))!!')'
*
    type 'IF AVG Off-source image statistics for source '!!inna
    blc rmsbox(1), rmsbox(2), 0; trc rmsbox(3), rmsbox(4), 0
    incl iclass; imstat; 
    type inna!!' IF=AVG I RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy '
*
  else indisk mapdis;
*
    type 'IF 1 On-source image statistics for source '!!inna
    blc absbox(1), absbox(2), 0; trc absbox(3), absbox(4), 0
    incl iclass; inseq 1; imstat; 
    type inna!!' IF=1 I Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy '
    type inna!!' IF=1 I Max: at ('!!CHAR(pixxy(1))!!','!!CHAR(pixxy(2))!!')'
*
    type 'IF 1 Off-source image statistics for source '!!inna
    blc rmsbox(1), rmsbox(2), 0; trc rmsbox(3), rmsbox(4), 0
    incl iclass; imstat; 
    type inna!!' IF=1 I RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy '
*
* IF 2 data, INSEQ=2
*
    type 'IF 2 On-source image statistics for source '!!inna
    blc absbox(1), absbox(2), 0; trc absbox(3), absbox(4), 0
    incl iclass; inseq 2; imstat; 
    type inna!!' IF=2 I Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy '
    type inna!!' IF=2 I Max: at ('!!CHAR(pixxy(1))!!','!!CHAR(pixxy(2))!!')'
*
    type 'IF 2 Off-source image statistics for source '!!inna
    blc rmsbox(1), rmsbox(2), 0; trc rmsbox(3), rmsbox(4), 0
    incl iclass; imstat; 
    type inna!!' IF=2 I RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy '
  end
*
  Finish
*
*
Proc outmap(i)
*
* Maps source
*
  clrtemp
  prtas ''; clrmsg
  task 'imagr'
  indi 1
  getn i
  docalib 1; gainuse 0; dopol -1
*** This has to change for different configs and bands
  imsize 256; cell abcell
*** No need to clean the entire map either
  nbox numbox; clbox absbox, 0
***
  uvwtfn ' '
  cmethod ''; flagver 0; bchan 1; echan 0; subar 0; clr2n; nfield 1
  bmaj 0; bmin 0; bpa 0; cparm 0
  dotv true
  baddis 2 0
*
* Do ABSCAL first
  source abscal,''
  uvra 0 0
  zerosp 0
  outver 0
  tput imagr
  mapem
  outrms
  outprint predir!!source(1)!!'-'!!prefix!!abstag!!'.TXT'
  docrt -1; prtmsg; docrt 1; outprint ''
*
  type ' Completed Mapping '
*
Finish