* * CALCX - set of procedures to semi-automatically reduce polarization * calibration data for the VLA for dataset: * 09-30-2001 1100-1400 LST SMyers TESTT/AC524 program * * Version STM-10-03-01 for running on SANDROCK USER 1553 Disks 1 and 2 * Include VPol * THIS ASSUMES THE 31DEC01 VERSION OF AIPS OR LATER! * * NOTE: The session-specific variables are defined in procedure SETIT. * This script should be useable for any run if these are changed * appropriately. * * This script assumes the UV data and ICLN models are on Disk 1, and * Disk 2 is available for temporary map storage. * * Because the primary calibrator is 1331+305 (3C286) or 0137+331 (3C84) * which are resolved, a CC model is preferred. However, the following * procedure will sort of work with a point model by skipping step 2. * * Calibrators are broken into lists: * * abscal = primary flux and EVPA calibrator * polcal = polarization leakage calibrator(s) * * srclist = all calibrators except abscal * | * strongcal = bright calibrators from callist (shorter solutions) * phacal = normal calibrators from callist * weacal = weak calbirators from callist (longer solution intervals) * * 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 * * 2. Normally, at this point, you would split off the data for the cal * and map it in SCMAP or DIFMAP to produce a ICLN with CC=1. For * D-configuration (and probably C) a point model is sufficient. * * 3. CALDATA(I) -- gain calibration, using the model, produces CL 3 and * CL 4 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' * 5. CALPOL(I) -- runs PCAL and does a LISTR OPTYP='MATX' on cross-pols * for calibrator (RL and LR). * 6. Examine LISTR output files *.PCAL and *.POLC * 7. OUTPOL(I) -- runs IMAGR and then IMSTAT/IMVAL on the maps * 8. Examine LISTR output files *.TXT and produce summary by, for example * grep 'IF=' *.TXT * The lines at the end of the .TXT files contain the angles we want, * and can be extracted by * grep 'PHASE' *.TXT * * Remember, if you need to redo things, proc RESETIT(I) will reset to * the beginning, and to redo steps 5 or 7 just requires blasting the listr * output files. * Proc Setit * Procedure to define global pseudoverbs, arrays and strings * as well as to define the calibrator names and number for later use * scalar angle, angle1, angle2, polang1, polang2, corang1, corang2 scalar phdiff, flxtot1, flxtot2, fpol1, fpol2, spol1, spol2 * string*16 abscal, absnam, srccal, absnot, srclist(30) string*16 phacal(30), polcal(30), weacal(30), strongcal(30) string*12 absmod, predir string*6 iclass, bclass, uclass, qclass, vclass string*2 prefix string*1 fqband scalar abcell, absref, absang, numbox, calimode, sjparm scalar numcal, numsrcs, numweak, numstrong array absflx(2), absant(50), notant(50), absuvr(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'; uclass 'UCL001'; qclass 'QCL001'; vclass 'VCL001'; * * Calibrate weights (DOCALIB=2) * calimode 2; * * NOTE: These values are specific to this session * type 'Using setup values for DnC-config X-Band 09-39-2001 TESTT/AC524 run' * predir 'PWD:'; type 'PREDIR =',predir prefix 'CX'; type 'PREFIX =',prefix fqband 'X'; type 'FQBAND =',fqband abscal '1331+305'; type 'ABSCAL =',abscal; absnot '' absnam '3C286'; type 'ABSNAM =',absnam absmod '1331+305-AX'; type 'ABSMOD =',absmod absang 66; type 'ABSANG =',absang * setjy calc flux 5.1915,5.1702 (1995.2) * for point model flux reduced 1% absflx 0; sjparm 2; type 'ABSFLX =',absflx * use recommended 50-300 klam and all stations absuvr 50 300; type 'ABSUVR =',absuvr absant 0 notant 0 type 'ABSANT =',absant * absref 2; type 'ABSREF =',absref numbox 0; type 'NUMBOX =',numbox absbox 0, 0, 0, 0; type 'ABSBOX =',absbox rmsbox 32, 32, 96, 96; type 'RMSBOX =',absbox abcell 0.5; type 'ABCELL =',abcell * numsrcs 9; type 'NUMSRCS =',numsrcs srclist '0713+438','0854+201','0927+390','1146+399','1159+292','1229+020',''; srclist(7) '1256-057'; srclist(8) '1310+323'; srclist(9) '1337-129'; numcal 5 phacal '0854+201','0927+390','1159+292','1310+323','1337-129',''; type 'PHACAL =',phacal numstrong 2 strongcal '1229+020','1256-057',''; type 'STRONGCAL =',strongcal numweak 2 weacal '0713+438','1146+399',''; type 'WEACAL =',weacal polcal '1159+292',''; type 'POLCAL =',polcal srccal ' ' * 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 * 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 * Specialized procedure to run IMAGR with stokes I,Q, and U * sets, and map them. * tget imagr for j=1 to 2; outn sour(1) bif j; eif j outseq (j+6) clr2n; in2di indi; outdi 2 niter 200 stokes 'i'; outcl iclass wait imagr; go imagr niter 200 stokes 'q'; outcl qclass wait imagr; go imagr stokes 'u'; outcl uclass wait imagr; go imagr stokes 'v'; outcl vclass wait imagr; go imagr end wait imagr * now make polc and pola maps for j=1 to 2; indi 2; in2di 2; outdi 2 innam outn blc 0; trc 0 inseq (j+6); outseq inseq incl qclass in2name innam in2seq inseq in2cl uclass opcode 'polc' outcl opcode aparm 0 blc 70 30; trc 140 94; tvlod; imstat bparm 0 0 1 bparm(1)=pixstd incl uclass blc 70 30; trc 140 94; tvlod; imstat bparm(2)=pixstd incl qclass; blc 0; trc 0 go comb; wait comb opcode 'pola' outcl opcode go comb; wait comb end 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 * * First TACOP CL1 to CL 2 * task 'tacop' indi 1; getn i inext 'cl' outdi indi geto i inver 1 outver 2; ncount 1 keyword ''; keyval 0 go tacop; wait tacop inext ''; inver 0; outver 0; clroname * 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 1; cmodel 'comp'; cmethod 'dft'; invers 1; flux 0; ncomp 0; nmaps 1; uvr 0 0; anten 0; basel 0 else in2n ''; type 'Phase calibration using Point model - all data' in2cl ''; in2seq 0; in2d 0; cmodel ''; cmethod ''; invers 0; flux 0; ncomp 0; nmaps 0; smodel 0; uvr 0 0; anten 0; basel 0 end calsour abscal,'' flagver 0; gainuse 0; docalib calimode; refant absref aparm 4 0 cparm 0 0 10 10 0 * * Main calibrator at 20 sec * solmode 'p'; solint 20/60 snver 1 go calib; wait calib * * Normal sources at 20 sec * calsour phacal uvr 0 0; anten 0; basel 0 clr2n; invers 0; flux 0; ncomp 0; nmaps 0 solmode 'p'; solint 20/60 snver 1 go calib; wait calib * * Strong sources at 10 sec * if numstrong > 0 then calsour strongcal; calsour strongcal uvr 0 0; anten 0; basel 0 clr2n; invers 0; flux 0; ncomp 0; nmaps 0 solmode 'p'; solint 10/60 snver 1 go calib; wait calib end * * Weaker sources at 30 sec * if numweak > 0 then calsour weacal; uvr 0 0; anten 0; basel 0 clr2n; invers 0; flux 0; ncomp 0; nmaps 0 solmode 'p'; solint 30/60 snver 1 go calib; wait calib end * * Apply each source to itself * task 'clcal' type 'Applying Self Calibration for ',abscal source abscal,'' calsour source opcode ' ' interp 'simp' gainver 2 gainuse 3 go clcal; wait clcal * for i = 1 to numsrcs; srccal = srclist(i) type 'Applying Self Calibration for ',srccal task 'clcal' source srccal,'' calsour source opcode ' ' interp 'simp' gainver 2 gainuse 3 go clcal; wait clcal end * * 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 1; cmodel 'comp'; cmethod 'dft'; invers 1; flux 0; ncomp 0; nmaps 1; uvr 0 0; anten 0; basel 0 else in2n ''; type 'Amplitude calibration using Point model - restricted data' in2cl ''; in2seq 0; in2d 0; cmodel ''; cmethod ''; invers 0; flux 0; ncomp 0; nmaps 0; smodel 0; uvr absuvr; anten absant; basel anten end calsour abscal,'' flagver 0; gainuse 0; docalib calimode refant absref aparm 4 0 cparm 0 0 10 10 0 * * One minute gain solutions * solmode 'a&p'; solint 1 snver 2 go calib; wait calib * * Normal calibrators 1m average * calsour phacal uvr 0 0; anten 0; basel 0 clr2n; invers 0; flux 0; ncomp 0; nmaps 0 solmode 'a&p'; solint 1 snver 2 go calib; wait calib * * Strong calibrators with 30s solutions * if numstrong > 0 then calsour strongcal; uvr 0 0; anten 0; basel 0 clr2n; invers 0; flux 0; ncomp 0; nmaps 0 solmode 'a&p'; solint 30/60 snver 2 go calib; wait calib end * * Weak calibrators scan average only * if numweak > 0 then calsour weacal; uvr 0 0; anten 0; basel 0 clr2n; invers 0; flux 0; ncomp 0; nmaps 0 solmode 'a&p'; solint 0 snver 2 go calib; wait calib end * * Get fluxes from amplitude solutions * task 'getjy' source srclist calsour abscal,'' anten 0 bif 1; eif 2 go getjy * * 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 3 gainuse 4 go clcal; wait clcal else calsour source; type 'Applying Self Calibration for ',abscal opcode ' ' anten absant interp 'simp' gainver 3 gainuse 4 go clcal; wait clcal calsour absnot,'' type 'Applying External Calibration for ',abscal anten notant go clcal; wait clcal anten 0 end * for i = 1 to numsrcs; srccal = srclist(i) type 'Applying Self Calibration for ',srccal task 'clcal' source srccal,'' calsour source opcode ' ' anten 0 interp 'simp' gainver 3 gainuse 4 go clcal; wait clcal end * prtas '' outprint predir!!prefix!!'.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' inext 'sn';inver 2 docalib 1; gainuse 0; dopol -1 dparm 0 0 5 0 docrt -1 antenn 0; basel 0 * bif 1; eif 1; stokes 'r' outpr predir!!prefix!!'-1R.AMPL' go listr; wait listr * bif 1; eif 1; stokes 'l' outpr predir!!prefix!!'-1L.AMPL' go listr; wait listr * bif 2; eif 2; stokes 'r' outpr predir!!prefix!!'-2R.AMPL' go listr; wait listr * bif 2; eif 2; stokes 'l' outpr predir!!prefix!!'-2L.AMPL' 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 calpol(i) clrtemp * indisk 1 getn i task 'pcal' calsour polcal bif 1; eif 2 docalib 1; gainuse 0 solint 2 refant absref anten 0 soltyp ''; bparm 0; cparm 0 prtas ''; clrmsg go pcal; wait pcal prtas 'pcal' outprint predir!!prefix!!'.PCAL' docrt -1 prtms docrt 1 * task 'listr' sour abscal,'' uvr absuvr; anten absant; basel anten docalib 1; gainuse 0 dopol 1 optyp 'matx' stokes 'polc' dparm 5 0 docrt -1 outprint predir!!prefix!!'-'!!absnam!!'.POLC' bif 1 go listr; wait listr bif 2 go listr; wait listr docrt 1 uvr 0 0; anten 0; basel anten; stokes ' '; optyp ' ' * clrtemp getn i task 'listr' sour srclist uvra 0 0 docalib 1; gainuse 0 dopol 1 optyp 'matx' stokes 'polc' dparm 5 0 docrt -1 outprint predir!!prefix!!'-ALLCAL.POLC' bif 1 go listr; wait listr bif 2 go listr; wait listr docrt 1 * type ' Completed CALPOL ' * 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) blc rmsbox(1), rmsbox(2), 0; trc rmsbox(3), rmsbox(4), 0 * * IF 1 data, INSEQ=7 * inseq 7; type 'IF 1 Off-source image statistics for source '!!inna incl iclass; imstat; type inna!!' IF=1 I RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy ' incl qclass; imstat; type inna!!' IF=1 Q RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy ' incl uclass; imstat; type inna!!' IF=1 U RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy ' incl vclass; imstat; type inna!!' IF=1 V RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy ' * * IF 2 data, INSEQ=8 * inseq 8; type 'IF 2 Off-source image statistics for source '!!inna incl iclass; imstat; type inna!!' IF=2 I RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy ' incl qclass; imstat; type inna!!' IF=2 Q RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy ' incl uclass; imstat; type inna!!' IF=2 U RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy ' incl vclass; imstat; type inna!!' IF=2 V RMS: FLUX = '!!CHAR(pixstd*1000.0)!!' milliJy ' * Finish * * Proc setang (angle) * * Set vector angle into -180 to 180 deg range of validity * if angle < -180.0 then angle = angle + 360.0; end if angle > 180.0 then angle = angle - 360.0; end * Finish * * Proc outstat (angle1,angle2) * * short proc to use IMSTAT for current image to output the polarization * angles for IF 1 and 2 for the current source name SOURCE(1) * Should have already set BLC,TRC for window (or 0,0 for whole map) * * NOTE: Polarization angles output by COMB are 1/2 the R/L phase difference * and thus run from -90 to +90 degrees * clrtemp inna source(1); * * IF 1 data, INSEQ=7 * type 'IF 1 On-source image statistics for source '!!inna incl iclass; inseq 7; 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))!!')' flxtot1 = pixval*1000.0 * incl qclass; imval; type inna!!' IF=1 Q Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy ' incl uclass; imval; type inna!!' IF=1 U Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy ' incl vclass; imval; type inna!!' IF=1 V Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy ' incl 'POLC'; imval; type inna!!' IF=1 P Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy ' fpol1 = pixval*1000.0 incl 'POLC N'; imval; type inna!!' IF=1 P Max: Noise = '!!CHAR(pixval*1000.0)!!' milliJy ' spol1 = pixval*1000.0 incl 'POLA'; imval; type inna!!' IF=1 P Max: Angle = '!!CHAR(pixval)!!' Deg ' angle1 = pixval incl 'POLA N'; imval; type inna!!' IF=1 P Max: Error = '!!CHAR(pixval)!!' Deg ' * type inna!!' IF=1 Poln Flux = '!!CHAR(fpol1)!!' +/- '!!CHAR(spol1)!!' mJy ' type inna!!' IF=1 Frac Poln = '!!CHAR(100.0*fpol1/flxtot1)!!' % ' * * IF 2 data, INSEQ=8 * inna source(1); type 'IF 2 On-source image statistics for source '!!inna incl iclass; inseq 8; 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))!!')' flxtot2 = pixval*1000.0 * incl qclass; imval; type inna!!' IF=2 Q Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy ' incl uclass; imval; type inna!!' IF=2 U Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy ' incl vclass; imval; type inna!!' IF=2 V Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy ' incl 'POLC'; imval; type inna!!' IF=2 P Max: Flux = '!!CHAR(pixval*1000.0)!!' milliJy ' fpol2 = pixval*1000.0 incl 'POLC N'; imval; type inna!!' IF=2 P Max: Noise = '!!CHAR(pixval*1000.0)!!' milliJy ' spol2 = pixval*1000.0 incl 'POLA'; imval; type inna!!' IF=2 P Max: Angle = '!!CHAR(pixval)!!' Deg ' angle2 = pixval * incl 'POLA N'; imval; type inna!!' IF=2 P Max: Error = '!!CHAR(pixval)!!' Deg ' * type inna!!' IF=2 Poln Flux = '!!CHAR(fpol2)!!' +/- '!!CHAR(spol2)!!' mJy ' type inna!!' IF=2 Frac Poln = '!!CHAR(100.0*fpol2/flxtot2)!!' % ' * * Compute phase diff between IF 1 and 2 phdiff = 2.0*(angle1 - angle2); setang( phdiff ) type inna!!' IF=1-2 RL Phase Difference = '!!CHAR(phdiff)!!' Deg ' * Finish * * Proc outpol(i) * * Maps sources, calculates polarizations and angles and writes them out * clrtemp 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 false baddis 2 0 * * Do ABSCAL first to establish angle setting source abscal,'' uvra 0 0 zerosp 0 outver 0 tput imagr mapem indi 2 prtas ''; clrmsg outrms blc absbox(1), absbox(2), 0; trc absbox(3), absbox(4), 0 outstat(angle1, angle2) * These angles are 1/2*RL phase diff * At X band the RL phase of ABSCAL is ABSANG deg corang1 = absang - 2.0*angle1; setang( corang1 ) corang2 = absang - 2.0*angle2; setang( corang2 ) * These correction factors are the instrumental R/L phase difference * correction factors ANGLE(TRUE) - ANGLE(OBS) * so we can just add these to the observed R/L angles for the sources type source(1)!!' IF=1 RL Phase Angle = '!!CHAR(absang)!!' Deg ' type source(1)!!' IF=2 RL Phase Angle = '!!CHAR(absang)!!' Deg ' type source(1)!!' IF=1 RL Phase Correction = '!!CHAR(corang1)!!' Deg ' type source(1)!!' IF=2 RL Phase Correction = '!!CHAR(corang2)!!' Deg ' outprint predir!!source(1)!!'-'!!prefix!!'.TXT' docrt -1; prtmsg; docrt 1; outprint '' for j=1 to 24; getn j; zap; end * * Loop over other calibrators * for i = 1 to numsrcs; * clrtemp tget imagr srccal = srclist(i) * type 'Imaging data from source ',srccal * source srccal,'' indi 1 tput imagr mapem indi 2 prtas ''; clrmsg outrms blc absbox(1), absbox(2), 0; trc absbox(3), absbox(4), 0 outstat(angle1, angle2) polang1 = 2.0*angle1 + corang1; setang( polang1 ) polang2 = 2.0*angle2 + corang2; setang( polang2 ) type source(1)!!' IF=1 RL Phase Angle = '!!CHAR(polang1)!!' Deg ' type source(1)!!' IF=2 RL Phase Angle = '!!CHAR(polang2)!!' Deg ' outprint predir!!source(1)!!'-'!!prefix!!'.TXT' docrt -1; prtmsg; docrt 1; outprint '' for j=1 to 24; getn j; zap; end * end * type ' Completed Polarization Calibration Run ' type ' Results can be found in *'!!prefix!!'.TXT files ' * Finish