#! /bin/csh -f # handling standard 4 GHz polarization data # # Responsible: testing version to be used for 4GHz polarization system # Version 1.0 # testing data: 100505_04:45:58 from SMA Ken (Taco) Young, PI: Keping Qiu # The testing data set has not been available to the public yet. # Please acquire the PI's permission for use of the testing data # and the figures produced from this script. # # demo illustration figures (under construction) # fig1.ps = image of stokes I, total intensity # fig2.ps = image of stokes Q # fig3.ps = image of stokes U # fig4.ps = image of polarization vector (linear lines) # overlaid on polarization intensity (greyscaled) # #################################################################################### # #jhz 2011-02-17 received the test data from Taco #jhz 2011-02-17 received the c-code of extract baseline-based # information regarding to the polarization visibility # #jhz 2011-02-17 communicated with users for information on the 4GHz pol data and test # #jhz 2011-02-21 look at the data with mirREAD routines outside Miriad #jhz 2011-03-02 come back to look into the polarization section with mirRead # routines for 4 GHz polarization data; pol states and pol # sampling patterns are difficult to understand; the documents # available appeared to inadequate for disentangle them. # #jhz 2011-03-03 start the script (Ver0.1) for testing the 4GHz polarization # system # #jhz 2011-03-03 work 4GHz polarization codes in Miriad #jhz 2011-03-23 go through compiling #jhz 2011-03-25 smalod reads all 4GHz polarization data #jhz 2011-03-29 check the variables for polarization calibration # realized a few peculiar problems. #jhz 2011-03-30 parallactic angles for the 4GHz data are inconsistent between # the values from smalod and uvredo .... #jhz 2011-03-31 polarization labelling for the 4GHz data does not make sense. # check evector OK # check mount type OK #jhz 2011-04-04 install SMA Miriad beta for 4 GHz pol test code on Jupiter #jhz 2011-04-05 make a tarball for the 4GHz pol test Miriad to SMA Miriad # web site for testers # #test with 2GHz polarization data #jhz 2011-04-05 check the polarization sampling pattern # check the time steps on the sources # check the e-vector and mount type # check the parallactic angle coverage of # pol calibarators #jhz 2011-04-07 check gain calibrators # check the algorithm in solving for D-terms # and the error assessments #jhz 2011-04-08 worked on the pol data editing tool programs # fixed a few bugs #jhz 2011-04-11 test the calibration with the data from the GC pol # program. # ....... #jhz 2011-04-28 received 7 variable tables of the 4GHz pol test data set # 100505_04:45:58 from online software crew: # .... # Attached are seven tables. The first gives the scan number, # source name, parallactic angle (degrees), integration time # (seconds) and UT (hours) for each scan. Each of the other # 6 tables # gives the polarization state for one baseline to antenna 1 # (table2 has 1-2 data, table 3 has 1-3 data etc., there is # no table for 1-7 as 7 was not in the array). For the # last 6 tables the format is scan number, source name, # baseline name ("1-2" etc.) and polarization code. The # polarization state for any baseline can be derived from # the information in these tables. All of these tables # contain entries for every scan, including the scans # which are flagged bad on one or more baselines. # Taco # #jhz 2011-05-10 verified with Ken Young, bl.ipol, the baseline pol ids in # the Mir formatted data for 4 GHz polarization mode with # QWP on single receivers are related to circular pol states: # bl.ipol pol states # 0 Unknown # 1 RR # 2 RL # 3 LR # 4 LL # #jhz 2011-05-11 merge 2GHz script Ver1.2 with 4GHz Ver0.1 to Ver0.2 # correct pol labelling and several other minor bugs. # recompile the package and upgrate to beta 1.3.0 # install Ver 1.3.0 on Jupiter (centOS) and bigfoot (Redhat) # test smalod on both Jupiter bigfoot successfully. # update SMA Miriad web for testers downloading # #jhz 2011-05-12 add back chi2 in variables # add lineset, script variable for linetype # change stokes = i to =ii in uvflux # add smasplt # add a few script variables, nw, spltopt # add goto tag SMASPLT for 4GH data # add gcal2 and fcal2 # #jhz 2011-05-12 add goto tag BLCORR to correct for baseline errors # work out setup to correct for baseline errors. # #jhz 2011-05-13 add goto tag CHI #jhz 2011-05-18 update ephemeris for the solar objects used in smaflux #jhz 2011-05-19 update sma/mount = NASMYTH #jhz 2011-06-03 implement the imaging part #jhz 2011-06-09 changed the data path to CF where the test data located. #jhz 2011-06-13 following the suggestions from the software meeting # today, comment out the blflag and add a note for instruction. # $< echo 'Press return to continue' $< ########################################################### #enable the following line: 'goto SMATEST' #allow users to select their own visibility data #for the reduction of polarization data with #this script. ########################################################### goto SMATEST if($#argv<2) then echo "Incorrect usage: Needs Task Date Time" echo "example: sma_pol_cal.csh 100505 04:45:58" goto terminate endif ########################################## #--> Define date, time ofthe observations ########################################## set date=$1 set hour=$2 goto USERSREDUCTION SMATEST: ########################################## # for SMA software test ########################################## set date=100505 set hour=04:45:58 USERSREDUCTION: ########################################## #--> Define side band (calibration is done # separately for the USB and LSB) # nn=0 for LSB, nn=1 for USB ########################################## set sb = usb set nn = 1 set rx = rx1 set fname = "$date""_""$rx".$sb echo $fname ########################################## #--> Define the Flux, Gain, Bandpass and # Polarization Calibrators used ########################################## set fcal = titan set fcal2 = neptune set gcal = 1830+063 set gcal2 = 1751+096 set bcal = 3c273 set pcal = 3c273 set prcal = 3c454.3 ######################################### #--> Define the Source name ######################################### set font = I18360 # select source for pol mapping set source = '' ######################################### #--> Define reference antenna # average time interval ######################################### set refant = 5 set avetime = 6 ######################################### #--> Other definitions set region = quart set regdisp = 'arcsec,box(-12,-12,12,12)' set dt = flx set spltops = all set nw = 48 set lineset = 'chan,48,6,116,128' set edg = '6,6' set bc = blc # ######################################### # #goto continue # ######################################### # enable one of the next several goto # lines allow users to restart the script # from a working step. ######################################### # # #goto SMALOD #Load smadata #goto INSPECT #Inspect visdata #goto IPOLFLAG #Flag data with incorrect pol tags #goto TSYSFIX #Corrections for fringe amplitude #goto SMASPLT #Separate the ch0 data from spectral chunk data #goto FLAG #Flag bad visibilities #goto BPCAL #Bandpass calibration #goto APPLYBP #Apply bandpass solutions #goto CHI #re-calculate chi with uvredo #goto FLUX #Determine flux density scale #goto SPLIT #Split the data into single-source file #goto BLCORR #Corrections for baseline errors #goto GETFLUX #Get the flux density of the pol-calibrator #goto POLCAL #Do polarization calibration #goto GCAL #Gain calibration for target #goto AVER #Make averaging visibility to a longer time interval #goto POLMAP #Start polarization mapping # #users: please change the data path below set DATPATH = /sma/SMAusers/jzhao/data/ # ##################################### # copy the antenna file to working area ##################################### # cp ${DATPATH}/$date"_"$hour/antennas . # ##################################### # Convert raw data into miriad format ##################################### # SMALOD: echo $fname \rm -fr $fname smalod in=${DATPATH}/$date"_"$hour out=$date rxif=1 \ options=circular sideband=$nn refant=${refant} nscans=5, #goto terminate $< echo 'Press return to continue' $< # ##################################### # Inspect at ##################################### INSPECT: uvindex vis=$fname $< smacheck vis=$fname var=systemp flagval=flag range=0,1000 $< smauvplt vis=$fname axis=time,ampl device=/xs nxy=1,1 $< #goto terminate IPOLFLAG: ##################################### # Flag data with incorrect pol tags ##################################### uvflag vis=$fname select='pol(i)' flagval=flag \rm -fr $fname.flgd uvflag vis=$fname select='-pol(ll,lr,rl,rr)' flagval=f uvcal vis=$fname options=unflagged out=$fname.flgd $< #goto terminate ##################################### # Tsystem correction: ##################################### TSYSFIX: \rm -fr $fname.tsys smafix vis=$fname.flgd out=$fname.tsys device=/xw nxy=1,1 \ yrange=100,800 options=tsyscorr,dosour #goto terminate ##################################### #Separate the ch0 from spectral windows ##################################### SMASPLT: \rm -r $fname.tsys.$nw smasplt vis=$fname.tsys out=$fname.tsys.$nw options=$spltops ##################################### # Flag corrupted data: check your # visibility data to look for bad data ##################################### FLAG: ##################################### #flag edge channels ##################################### uvflag vis=$fname.tsys.$nw flagval=f edge=$edg $< ##################################### #flag bad time interval for selected source ##################################### uvflag vis=$fname.tsys.$nw flagval=f select='time(13:59,14:11),sou('$prcal')' uvflag vis=$fname.tsys.$nw flagval=f select='time(06:00,06:25),sou('$bcal')' $< ##################################### #flag bad antennas ##################################### uvflag vis=$fname select='ant(8),time(12:30:00,13:20:00)' flagval=f $< # # users can enable blflag for flagging the bad visibilities. # #blflag vis=$fname.tsys.$nw axis=time,ampl device=/xs select='sou('$bcal')' # #goto terminate ##################################### # Bandpass calibration ##################################### BPCAL: set options=averrll smamfcal select='pol(ll,rr),sou('$bcal')' refant=$refant \ interval=3000 weight=2 \ vis=$fname.tsys.$nw options=$options edge=$edg $< smagpplt vis=$fname.tsys.$nw device=/xs \ yaxis=amp,phase \ options=bandpass polyfit=5 \ nxy=1,8 #goto terminate ##################################### # Apply bandpass calibration ##################################### APPLYBP: \rm -fr $fname.tsys.$nw.bp uvaver vis=$fname.tsys.$nw out=$fname.tsys.$nw.bp options=nocal #goto terminate ################################################################# # # Define the correct quarter wave plate paramenters. # these parameters have been coded in both Miriad-CVS # and SMA beta Miriad # since Miriad-CVS 4.0.3 # #EVEC: #puthd in=$fname.tsys.$nw.bp/evector value=0.7853 # which has been coded in Miriad #puthd in=$fname.tsys.$nw.bp/mount value=4 # which has been coded in Miriad CHI: ################################ # re-calculate chi with uv-redo ################################ set chi = '' #set chi = '.chi' #\rm -fr $fname.tsys.$nw.bp$chi #uvredo vis=$fname.tsys.$nw.bp options=chi out=$fname.tsys.$nw.bp$chi # which has been coded in Miriad # #goto terminate ################################################################# # ################################################################# # Flux Calibration ################################################################# # FLUX: echo 'FLux calibrator =' $fcal smaflux vis=$fname.tsys.$nw.bp$chi select='pol(ll,rr),so('$fcal')' mirhome=$MIR #goto terminate # $< echo 'Press return to continue' $< # ################################################################# # Split visibility on source based. ################################################################# #goto terminate # SPLIT: set so=$fcal \rm -fr "$so""_""$dt".$sb uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \ out="$so""_""$dt".$sb line=$lineset set so=$fcal2 \rm -fr "$so""_""$dt".$sb uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \ out="$so""_""$dt".$sb line=$lineset set so=$gcal \rm -fr "$so""_""$dt".$sb uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \ out="$so""_""$dt".$sb line=$lineset set so=$gcal2 \rm -fr "$so""_""$dt".$sb uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \ out="$so""_""$dt".$sb line=$lineset set so=$bcal \rm -fr "$so""_""$dt".$sb uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \ out="$so""_""$dt".$sb line=$lineset set so=$pcal \rm -fr p"$so""_""$dt".$sb uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \ out=p"$so""_""$dt".$sb line=$lineset set so=$prcal \rm -fr "$so""_""$dt".$sb uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \ out="$so""_""$dt".$sb line=$lineset set so=$font \rm -fr "$so""_""$dt".$sb uvaver vis=$fname.tsys.$nw.bp$chi select='sou('$so')' \ out="$so""_""$dt".$sb line=$lineset #goto terminate BLCORR: \rm -rf p$pcal"_"$dt.$sb.$bc uvedit vis=p$pcal"_"$dt.$sb source=$pcal apfile=antennas \ smaoffset=1,-0.085,0.036,0.046,2,0.110,-0.092,0.049,3,-0.067,0.049,0.046,4,-0.088,-0.025,0.032,5,0.088,-0.032,0.064 \ out=p$pcal"_"$dt.$sb.$bc options=sma #goto terminate GETFLUX: ################################################### # Remove the phase error ################################################### selfcal vis=p$pcal"_"$dt.$sb.$bc refant=${refant} \ options=pha select='pol(ll,rr)' \ interval=${avetime} ################################################### # get the flux density for the pol calibrator ################################################### uvflux vis=p$pcal"_"$dt.$sb.$bc stokes=ii ################################################### # Averaging visibility to a long time interval ################################################### \rm -fr pp$pcal"_"$dt.$sb.$bc uvaver vis=p$pcal"_"$dt.$sb.$bc out=pp$pcal"_"$dt.$sb.$bc \ interval=${avetime} #goto terminate ################################################### # Polarization calibration ################################################### POLCAL: $< gpcal vis=pp$pcal"_"$dt.$sb.$bc refant=${refant} \ options=noxy,circular,qusolve \ interval=${avetime} flux=5.5 #goto terminate ################################################### # Gain Calibation for target. # "smamfcal options=nopassol" can also be used # instead of selfcal ################################################### GCAL: $< selfcal vis=$gcal"_"$dt.$sb refant=${refant} \ options=amp select='pol(ll,rr)' interval=${avetime} $< smagpplt vis=$gcal"_"$dt.$sb device=/xw yaxis=phas \ nxy=2,4 #goto terminate ################################################### # Applying Gain & Pol calibration to the target ################################################### gpcopy vis=$gcal"_"$dt.$sb out=$font"_"$dt.$sb mode=copy $< \rm -rf $font"_"$dt.$sb.g uvaver vis=$font"_"$dt.$sb out=$font"_"$dt.$sb.g interval=${avetime} $< gpcopy vis=pp$pcal"_"$dt.$sb.$bc out=$font"_"$dt.$sb.g options=nocal #goto terminate AVER: \rm -fr p$font"_"$dt.$sb.g uvaver vis=$font"_"$dt.$sb.g out=p$font"_"$dt.$sb.g $< #goto terminate ############################################################### # Gain calibration for a point source and apply pol calibration ############################################################### selfcal vis=$prcal"_"$dt.$sb refant=${refant} options=pha \ select='pol(ll,rr)' interval=0.5 $< \rm -fr $prcal"_"$dt.$sb.g uvaver vis=$prcal"_"$dt.$sb out=$prcal"_"$dt.$sb.g \ interval=${avetime} gpcopy vis=pp$pcal"_"$dt.$sb.$bc out=$prcal"_"$dt.$sb.g \ options=nocal $< \rm -fr p$prcal"_"$dt.$sb.g uvaver vis=$prcal"_"$dt.$sb.g out=p$prcal"_"$dt.$sb.g \ interval=${avetime} #goto terminate ################################################## echo ' Polarization Imaging' ################################################## POLMAP: ######################### # select source for mapping ######################### set source=$pcal ######################### # select pol visdata set ######################### set polvis=pp$pcal"_"$dt.$sb.$bc ######################### # setup the imaging size ######################### set imsize=512 ######################### # setup the cell size ######################### set cell=0.1 ######################### # setup clean region ######################### set clreg='arcsec,box(-3,-3,3,3)' ######################### # setup display region ######################### set disreg='arcsec,box(-8,-8,8,8)' ######################### # setup pol display region ######################### set pdisreg='arcsec,box(-4,-4,4,4)' ######################### # setup the area for statistics ######################### set statreg='arcsec,box(-10,-10,-8,-8)' set statcet='arcsec,box(-2,-2,2,2)' #goto terminate echo 'POLARIZATION IMAGING' $< echo 'Press return to continue' $< DIRTYMAP: if ($source == '') then goto POLMAP endif \rm -fr $source.$sb.i $source.$sb.q $source.$sb.u $source.$sb.v $source.$sb.b ######################################## #make sure to have a right vis file #corresponding to the polarization source #for imaging -> invert vis=???????????? ######################################## # ####################################### echo ' Make dirty images for four stokes' ####################################### # $< invert vis=$polvis beam=$source.$sb.b \ imsize=$imsize cell=$cell \ map=$source.$sb.i,$source.$sb.q,$source.$sb.u,$source.$sb.v \ options=mfs,sdb,systemp sup=0 \ stokes=i,q,u,v robust=2 line=chan,$nw,1,1 $< cgdisp in=$source.$sb.i,$source.$sb.i type=p,c \ device=dirty.map/cps \ labtyp=arcsec nxy=1,1 slev=p,10 \ levs1=2,3,4,5,6,7,8,9 \ options=wedge,mirror region=$disreg #goto terminate ####################################### echo ' Clean dirty images for four stokes' ####################################### $< CLEANMAP: if ($source == '') then goto POLMAP endif \rm -fr $source.$sb.icln $source.$sb.qcln $source.$sb.ucln $source.$sb.vcln \rm -fr $source.$sb.icmp $source.$sb.qcmp $source.$sb.ucmp $source.$sb.vcmp mfclean map=$source.$sb.i beam=$source.$sb.b out=$source.$sb.icmp \ region=$clreg gain=0.05 niters=150 $< mfclean map=$source.$sb.q beam=$source.$sb.b out=$source.$sb.qcmp \ region=$clreg gain=0.05 niters=100 $< mfclean map=$source.$sb.u beam=$source.$sb.b out=$source.$sb.ucmp \ region=$clreg gain=0.05 niters=100 $< mfclean map=$source.$sb.v beam=$source.$sb.b out=$source.$sb.vcmp \ region=$clreg gain=0.05 niters=100 #goto terminate ####################################### echo ' Restore cleaned pixel images for' echo ' four the stokes i, q, u, and v' ####################################### $< RESTOR: if ($source == '') then goto POLMAP endif $< restor map=$source.$sb.i beam=$source.$sb.b out=$source.$sb.icln \ model=$source.$sb.icmp $< restor map=$source.$sb.q beam=$source.$sb.b out=$source.$sb.qcln \ model=$source.$sb.qcmp $< restor map=$source.$sb.u beam=$source.$sb.b out=$source.$sb.ucln \ model=$source.$sb.ucmp $< restor map=$source.$sb.v beam=$source.$sb.b out=$source.$sb.vcln \ model=$source.$sb.vcmp $< #goto terminate DISPLAYSTOKES: if ($source == '') then goto POLMAP endif cgdisp in=$source.$sb.icln,$source.$sb.icln type=p,c device=fig1.ps/cps \ labtyp=arcsec nxy=1,1 slev=p,10 levs1=2,3,4,5,6,7,8,9 \ options=wedge,mirror,blacklab cols1=7 \ region=$disreg range=0.1,4,lin,2 $< cgdisp in=$source.$sb.qcln,$source.$sb.icln type=p,c device=fig2.ps/cps \ labtyp=arcsec nxy=1,1 slev=p,10 levs1=2,3,4,5,6,7,8,9 \ options=wedge,mirror,blacklab cols1=7 \ region=$disreg range=0.003,0.024,lin,2 $< cgdisp in=$source.$sb.ucln,$source.$sb.icln type=p,c device=fig3.ps/cps \ labtyp=arcsec nxy=1,1 slev=p,10 levs1=2,3,4,5,6,7,8,9 \ options=wedge,mirror,blacklab cols1=7 \ region=$disreg range=-0.003,-0.024,lin,2 $< cgdisp in=$source.$sb.vcln,$source.$sb.icln type=p,c device=v.map/cps \ labtyp=arcsec nxy=1,1 slev=p,10 levs1=2,3,4,5,6,7,8,9 \ options=wedge,mirror region=$disreg #goto terminate ######################################################## echo 'Make statistics from the images of stokes i, q, and u' ######################################################## $< STATISTICS: if ($source == '') then goto POLMAP endif imstat in=$source.$sb.icln region=$statreg $< imstat in=$source.$sb.qcln region=$statreg $< imstat in=$source.$sb.ucln region=$statreg $< imstat in=$source.$sb.vcln region=$statreg $< imstat in=$source.$sb.qcln region=$statcet $< imstat in=$source.$sb.ucln region=$statcet #goto terminate ######################################################## echo ' Make polarization intensity and vector images' echo ' from the cleaned images of stokes i, q, and u' ######################################################## $< POLVEC: if ($source == '') then goto POLMAP endif $< set irms=0.01 set qurms=0.005 \rm -fr $source.$sb.pi $source.$sb.pie \rm -fr $source.$sb.pm $source.$sb.pme \rm -fr $source.$sb.pa $source.$sb.pae $< impol in=$source.$sb.qcln,$source.$sb.ucln,$source.$sb.icln \ poli=$source.$sb.pi,$source.$sb.pie \ polm=$source.$sb.pm,$source.$sb.pme \ pa=$source.$sb.pa,$source.$sb.pae \ sigma=$qurms,$irms sncut=3,10 $< cgdisp in=$source.$sb.pi,$source.$sb.pm,$source.$sb.pa \ type=pix,amp,ang device=fig4.ps/cps \ labtyp=arcsec options=mirror,blacklab \ vecfac=2.0,3.0 cols1=0.5 \ range=0,0,lin,1 \ lines=5,8,8 region=$pdisreg terminate: echo 'end' exit