#! /bin/csh -f # handling standard 2 GHz polarization data: # Demonstrating to make polarization images of # a target source from loading SMA data, calibration # to imaging in Miriad. # # Responsible: Josep Miquel Girart (CSIC-IEEC) & Jun-Hui Zhao (SAO) # Version 2.0 # testing data: 070728_02:39:09 from SMA archive # # demo illustration figures # fig1.ps = image of polarization vector (linear lines) # overlaid on polarization intensity (greyscaled) # and continuum emission (contours) for g3141 # ################################################################## ################################################################# # # The correct quarter wave plate paramenters have been coded in # both Miriad-CVS and SMA beta Miriad since Miriad-CVS 4.0.3 # ################################################################# echo ' ' echo '************************************************************' echo 'Acknowledgement:' echo 'The reduction procedure for SMA polarization data written in' echo 'the following script was based on the communications with' echo 'Josep Miquel Girart, Meredith Hughes, Dan Marrone, Ram Rao,' echo 'Bob Sault and Mel Wright.' echo ' ' echo 'For Further Information on The SMA Polarimeter' echo 'Please read the article written by Dan Marrone and Ram Rao' echo 'http://adsabs.harvard.edu/abs/2008SPIE.7020E..60M' echo 'and the Ph.D. Thesis of Dan Marrone.' echo '************************************************************' echo ' ' ################################################################## # #jhz 2011-04-04 start the script (Ver1.0) for reduction of 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-20 contact SMA users for scripts of polarization data # reductions. Get positive responses. # ....... # #jhz 2011-05-04 receive a pol-data reduction script from Josep Miquel Girart # #jhz 2011-05-04 merge the two scripts to Ver 1.1 # change the number of argv to 2 # make an e2e script # check variable setups # comment out puthd for setup of QWP evector and type of mount # comment out uvredo for recalculate parallatic angles # which have been coded in both Miriad-CVS and SMA beta Miriad # since MIRIAD CVS 4.0.3 # #jhz 2011-05-05 correct fcal input # correct pcal input # correct fname setup # define variable dt # comment out foreach loop which was failed # split the individual sources and calibrators separately # change variable name ant to refant # add uvflag to flag pol(i) # add variable prcal for reference point pol source # remove interval=5 in the uvaver in splitting the sources # before gain calibration # change to line=ch,24,7,114,128 # add the variable avetime. # #jhz 2011-05-06 reset the interval 30 to 3000 in smamfcal in # solving for bandpass # add selfcal to calibrate the gains before solving for # D-terms. # add uvflux to calculate the flux density # of pol-calibrator that is used in GPCAL. # test the polarization calibration with a reference # point source. # receive a part of polarization imaging. # #jhz 2011-05-09 add the instruction lines for polimage # work out imaging for four stokes. # add the statistics for stokes images # derive the rms for i and sqrt(q*u) which # are used in impol. # work out imaging for pol intensity # and pol vector. # create four postscript figures for illustration. # make Ver1.2 # #jhz 2011-05-10 send to pol experts/users for review and comments # #jmg 2011-05-27 Major rearrange. The calibration steps follow a more # detailed pathway. Added menus to select different options # for the different steps. # #jhz 2011-06-02 add DATPATH #jhz 2011-06-02 test this script and made a few minor corrections # delete "exit" before goto SMALOD in the MAIN MANU # make hard copy fig1.ps for polarization vector # map of the target source. #jhz 2011-06-03 change the script name sma2GHz1Rx_pol_src.csh # upgrade to 2.0 #jhz 2011-06-09 changed the data path to CF where the test data located. echo '' echo 'Press return to continue' $< ########################################## #--> Define observational parameters: # date (yymmdd) # time (hh:mm:ss) # side band (nn=0 for LSB, nn=1 for USB)** # reciever: rx1 for 345 Ghz # # ** the calibration is done separately # for the USB and LSB ########################################## set dt = 070728 set hour = 02:39:09 set sb = usb set nn = 1 set rx = rx1 set fname = "$dt""_""$rx".$sb echo $fname ########################################## #--> Define the Flux, Gain, Bandpass and # Polarization Calibrators used ########################################## set fcal = uranus set gcal = 1749+096 set bcal = 3c454.3 set pcal = 3c454.3 set prcal = 3c273 ######################################### #--> Define the Source name ######################################### set font = g3141 ######################################### #--> Define reference antenna # average time interval ######################################### set refant = 2 set avetime = 5 ######################################### ########################################## #--> Define mappimg parameters: # polvis # source # imsize Image size in pixels # cell Cellsize in arcsecond # rob Robust weighting parameter # clreg Cleaning region. it can be a ascii file with many boxes or # polygons defined (e.g. clreg=@filename) # clcut1 Cutoff for cleaning Stokes I (typically 2-3 * sigma) # clcut2 Cutoff for cleaning Stokes Q & U # disreg Region for displaying the Stokes images # statreg Region for statistical analysis (e.g. determination of sigma) # irms rms noise of Stokes I # qurms rms noise of Stokes Q and U # ########################################## set source = $font set polvis = $source"_"$dt.$sb.tsys.bp.flx.gcal.pol set imsize = 256 set cell = 0.15 set rob = 0.5 set clreg = 'arcsec,box(-3,-3,3,3)' set clcut1 = 0.10 set clcut2 = 0.02 set disreg = 'arcsec,box(-8,-8,8,8)' set statreg = 'box(20,20,80,230)' set irms = 0.050 set qurms = 0.008 ################################################################## # For more detailed and complete information on standard SMA # calibration, # please check: # http://www.cfa.harvard.edu/sma/miriad/miriad4sma/stepbystep.html ################################################################## # script based on the pol cal done of NGC 1333 IRAS 4A & G31.41+0.31: # Girart, Rao, Marrone 2006, Science, 313, 812 # Girart, et al. 2009, Science, 324, 1408 echo;echo echo "=====================================================================" echo "WARNING: This is a simplified script. It is recommended to check all" echo "different calibration steps using the appropiate visualization tools" echo "for the visibility data (some are already included). " echo "=====================================================================" echo;echo MENU: echo '';echo '' echo '==========================' echo '|| MAIN MENU ||' echo '==========================' echo '' echo '1. Load SMA data into miriad' echo '2. Calibration' echo '3. Mapping' echo '' echo ' Select one of the following options' set sel=$< switch ($sel) case 1: echo 'Loading data into miriad' goto SMALOD breaksw case 2: echo 'Standard Calibration' goto MENUcal breaksw case 3: echo 'Polarization map' goto MENUmap breaksw endsw goto MENU ##################################### # Convert raw data into miriad format ##################################### # SMALOD: # #users: please change the data path below # set DATPATH = /sma/SMAusers/jzhao/data # echo $fname \rm -fr $fname smalod in=${DATPATH}/$dt"_"$hour out=$dt rxif=1 \ options=circular sideband=$nn refant=${refant} nscans=5, echo 'Done. Press return to continue' $< goto MENU ##################################### ##################################### # # CALIBRATION PROCEDURES # ##################################### ##################################### MENUcal: echo;echo echo '-------------------------------' echo '=> SMA CAlibration procedure <=' echo '-------------------------------' echo;echo echo 'select (number) for the following options' echo ' 1. Inspect visibility data. Tsystem correction' echo ' 2. Flag bad data' echo ' 3. Bandpass and flux Calibration' echo ' 4. Gain Calibration for the source' echo ' 5. Polarization Calibration' echo ' 6. Go back to main menu' echo echo 'Select number: '; set sel=$<;echo MENUcal2: echo; echo; echo 'Options:' switch ($sel) case 1: echo 'INSPECT Create a file with a summary of the observations' echo ' Check for high Tsys and flag them' echo ' Plot visibility amplitude of all sources' echo 'TSYSFIX Plot Tsys. Tsys corrections to the visibilies' echo; echo 'Select option (INSPECT, TSYSFIX):' set lab1=$< set out=tsys breaksw case 2: echo 'FLAG1 Flag data with incorrect pol tags' echo ' Flag data with Tsys and flag them' echo 'FLAG2 (optional in case data is corrupted)' echo ' Flag interatively bad data' echo ' If data is corrupted in a given interval, or for a given ' echo ' antenna, baseline, etc uncomment uvflag and define properly' echo ' select and other parameters' echo 'APPLYF Apply flags to the data' echo; echo 'Select option (FLAG1, FLAG2, APPLYF):' set lab1=$< set lab=tsys set out=tsys.flgd breaksw case 3: echo 'BPCAL Bandpass calibration ' echo ' Visualize solutions and data' echo 'APPLYBP Apply the bandpass calibration' echo 'FLUX Flux calibration' echo 'SPLIT Split data' echo; echo 'Select option (BPCAL, APPLYBP, FLUX, SPLIT):' set lab1=$< set lab=tsys.flgd set out=tsys.flgd.bp breaksw case 4: echo 'GCAL Gain calibration' echo ' Visualize solutions and data' echo 'FLAG3 Flag data (only if necessary)' echo 'APPLYG Apply gain calibration to the target' echo; echo 'Select option (GCAL, FLAG3, APPLYG):' set lab1=$< set lab=tsys.bp.flx set out=tsys.bp.flx.gcal breaksw case 5: echo 'GCAL2 First selfcalibrate pol.calibrator' echo ' Visualize solutions and data' echo 'POLAVE Average data in 5 minutes interval for polcal' echo 'GETFLUX Get the flux density for the pol calibrator' echo 'POLCAL Derive the leakage solutions ' echo 'APPLYPOL Apply leakages to the target' echo; echo 'Select option (GCAL2, POLAVE, GETFLUX, POLCAL, APPLYPOL):' set lab1=$< set lab=tsys.bp.flx set out=tsys.bp.flx.gcal set out2=tsys.bp.flx.gcal.pol breaksw case 6: set lab1=MENU breaksw endsw goto $lab1 ##################################### # Inspect Data ##################################### INSPECT: uvindex vis=$fname >& Summary"_""$dt"".""$rx" echo 'Generated file called: Summary_'$dt'.'$rx echo 'Press return to continue'; $< smacheck vis=$fname var=systemp flagval=flag range=0,1500 echo 'Press return to continue'; $< smauvplt vis=$fname axis=time,ampl device=/xs nxy=1,1 echo 'Press return to continue'; $< ##################################### # Tsystem correction: ##################################### TSYSFIX: \rm -fr $fname.$out smafix vis=$fname out=$fname.$out device=/xw nxy=1,1 \ yrange=100,800 options=tsyscorr,dosour echo 'Press return to continue'; $< goto MENUcal ##################################### # Flag data with incorrect pol tags # Flag ##################################### FLAG1: \rm -fr $fname.flgd uvflag vis=$fname.$lab select='-pol(ll,lr,rl,rr)' flagval=f uvflag vis=$fname.$lab flagval=f edge=7 echo 'Press return to continue'; $< goto MENUcal2 ##################################### # Flag corrupted data: check your # visibility data to look for bad data ##################################### FLAG2: #uvflag vis=$fname.$lab flagval=f select='time(13:59,14:11),sou('$bcal')' foreach src ($bcal $gcal) echo; echo 'Clik h on the screen for help' blflag vis=$fname.$lab axis=time,ampl device=/xs select='sou('$src')' end goto MENUcal2 ##################################### # Apply flags to data. create a new file ##################################### APPLYF: \rm -fr $fname.$out uvcal vis=$fname.$lab options=unflagged out=$fname.$out echo 'Press return to continue'; $< goto MENUcal ##################################### # Bandpass calibration ##################################### BPCAL: set options=averrll smamfcal select='pol(ll,rr),sou('$bcal')' refant=$refant \ interval=3000 weight=2 \ vis=$fname.$lab options=$options edge=7 echo 'Press return to continue. Check up the bandpass solutions' $< smagpplt vis=$fname.$lab select='pol(ll,rr),sou('$bcal')' \ device=/xs yaxis=phase,amp \ nxy=2,4 options=bandpass polyfit=5 echo 'Press return to continue'; $< echo echo 'Compare the raw and calibrated spectra (amplitude and phase).' echo 'First raw data' echo smauvspec vis=$fname.$lab select='pol(ll,rr),sou('$bcal')' \ device=1/xs interval=9999 axis=freq,both nxy=4,6 options=nopass echo 'Press return to continue. Now calibrated data' $< echo smauvspec vis=$fname.$lab select='pol(ll,rr),sou('$bcal')' \ device=2/xs interval=9999 axis=freq,both nxy=4,6 echo 'Press return to continue'; $< goto MENUcal2 ##################################### # Apply bandpass calibration ##################################### APPLYBP: \rm -fr $fname.$out uvaver vis=$fname.$lab out=$fname.$out options=nocal echo 'Press return to continue'; $< goto MENUcal2 ################################################################# # Flux Calibration ################################################################# FLUX: echo $fcal smaflux vis=$fname.$out select='pol(ll,rr),so('$fcal')' mirhome=$MIR echo 'Press return to continue'; $< goto MENUcal2 ################################################################# # Split visibility on source based. ################################################################# SPLIT: foreach so ($fcal $gcal $pcal $bcal $prcal $font) \rm -fr $so'_'$dt.$sb.tsys.bp.flx uvaver vis=$fname.$out select='sou('$so')' \ out=$so'_'$dt.$sb.tsys.bp.flx line=ch,24,7,114,128 end echo 'Press return to continue'; $< goto MENUcal ################################################### # Gain Calibation for target. ################################################### GCAL: set vis=$gcal'_'$dt.$sb.$lab smamfcal vis=$vis refant=${refant} \ options=nopassol select='pol(ll,rr)' interval=5 echo 'Plot amp & phase solutions. Press return to continue'; $< smagpplt vis=$vis device=1/xs yaxis=amp,pha nxy=2,4 echo 'Plot calibrated phase. Press return to continue'; $< smauvplt vis=$vis device=1/xs axis=tim,pha nxy=2,3 echo 'Plot calibrated amp. Press return to continue'; $< smauvplt vis=$vis device=1/xs axis=tim,pha nxy=2,3 echo 'Press return to continue'; $< goto MENUcal2 ##################################### # Flag corrupted data: check your # visibility data to look for bad data ##################################### FLAG3: set vis=$gcal'_'$dt.$sb.$lab echo; echo 'Clik h on the screen for help' blflag vis=$vis axis=time,ampl device=/xs options=noba blflag vis=$vis axis=time,phas device=/xs options=noba echo 'Press return to continue'; $< goto MENUcal2 ################################################### # Applying Gain calibration to the target ################################################### APPLYG: gpcopy vis=$gcal"_"$dt.$sb.$lab out=$font"_"$dt.$sb.$lab mode=copy \rm -rf $font"_"$dt.$sb.$out uvaver vis=$font"_"$dt.$sb.$lab out=$font"_"$dt.$sb.$out interval=${avetime} echo 'Press return to continue'; $< goto MENUcal ################################################### # Remove the phase error of pol calibrator ################################################### GCAL2: set vis=$pcal'_'$dt.$sb.$lab selfcal vis=$vis refant=${refant} \ options=pha select='pol(ll,rr)' \ interval=3 echo 'Plot phase solutions. Press return to continue'; $< smagpplt vis=$vis device=/xs yaxis=phas nxy=2,3 echo 'Plot calibrated phase. Press return to continue'; $< smauvplt vis=$vis device=/xs axis=time,phas nxy=2,3 echo 'Press return to continue'; $< goto MENUcal2 ################################################### # Averaging visibility to a long time interval ################################################### POLAVE: \rm -fr $pcal'_'$dt.$sb.$out uvaver vis=$pcal'_'$dt.$sb.$lab out=$pcal'_'$dt.$sb.$out \ interval=${avetime} echo 'Press return to continue'; $< goto MENUcal2 ################################################### # get the flux density for the pol calibrator ################################################### GETFLUX: uvflux vis=$pcal"_"$dt.$sb.$out stokes=i echo 'Press return to continue'; $< goto MENUcal2 ################################################### # Polarization calibration ################################################### POLCAL: set vis=$pcal'_'$dt.$sb.$out echo 'Delete bad points' blflag vis=$vis axis=real,imag options=nobase \ device=/xs echo 'Press return to continue'; $< gpcal vis=$vis refant=${refant} \ options=noxy,circular,qusolve \ interval=${avetime} flux=14.47 echo '' echo 'Pol calibration done. First visualize uncalibrated data. ' echo 'Press return to continue'; $< smauvplt vis=$vis device=1/xs axis=real,imag options=nocal,nopol,noba echo 'Visualize calibrated data. Press return to continue'; $< smauvplt vis=$vis device=2/xs axis=real,imag options=noba echo 'Press return to continue'; $< goto MENUcal2 ################################################### # Applying Pol calibration to the target ################################################### APPLYPOL: gpcopy vis=$pcal'_'$dt.$sb.$out out=$font'_'$dt.$sb.$out options=nocal \rm -rf $font'_'$dt.$sb.$out2 uvaver vis=$font'_'$dt.$sb.$out out=$font'_'$dt.$sb.$out2 interval=${avetime} echo 'Press return to continue'; $< goto MENUcal ##################################### ##################################### # # POLARIZATION IMAGING PROCEDURES # ##################################### ##################################### MENUmap: echo;echo;echo echo '-------------------------------' echo '=> SMA Polarization imaging <=' echo '-------------------------------' echo;echo echo; echo; echo 'Options:' echo 'DMAP Generate dirty map for all the Stokes' echo ' Display the dirty maps' echo 'CLEAN Make cleaned images all Stokes parameters' echo 'DISPLAY Display the cleaned maps for all Stokes parameters' echo 'STAT Make statistics from the Stokes images' echo 'POLVEC Compute polariation intensity and vector images' echo ' Display polarization vectors and intensity.' echo; echo 'Select option (DMAP, CLEAN, DISPLAY, STAT, POLVEC)' set lab1=$< goto $lab1 DMAP: ######################################## echo ' Make dirty images for four stokes' echo ' Press return to continue'; $< ####################################### \rm -fr $source.$sb.i $source.$sb.q $source.$sb.u $source.$sb.v $source.$sb.b 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,systemp \ stokes=i,q,u,v robust=$rob echo 'Press return to continue'; $< set dev=dirty_map.ps/cps set dev=/xs cgdisp in=$source.$sb.i,$source.$sb.i type=p,c \ device=$dev \ labtyp=arcsec nxy=1,1 slev=p,10 \ levs1=1,2,3,4,5,6,7,8,9 \ options=wedge,mirror region=$disreg echo 'Press return to continue'; $< goto MENUmap CLEAN: ####################################### echo ' Clean dirty images for four stokes' echo ' Press return to continue'; $< ####################################### \rm -fr $source.$sb.icmp clean map=$source.$sb.i beam=$source.$sb.b out=$source.$sb.icmp \ region=$clreg gain=0.1 niters=500 cutoff=$clcut1 foreach stk (q u v) \rm -fr $source.$sb.$stk'cmp' clean map=$source.$sb.$stk beam=$source.$sb.b out=$source.$sb.$stk'cmp' \ region=$clreg gain=0.05 niters=100 cutoff=$clcut2 end ####################################### echo ' Restore cleaned pixel images for' echo ' four the stokes i, q, u, and v' echo ' Press return to continue'; $< ####################################### foreach stk (i q u v) \rm -fr $source.$sb.$stk'cln' restor map=$source.$sb.$stk beam=$source.$sb.b out=$source.$sb.$stk'cln' \ model=$source.$sb.$stk'cmp' end echo 'Press return to continue'; $< goto MENUmap DISPLAY: ####################################### echo 'Display Stokes I image ' ####################################### set dev=$font'_'i.ps/cps set dev=/xs cgdisp in=$source.$sb.icln,$source.$sb.icln type=p,c \ device=$dev \ labtyp=arcsec nxy=1,1 slev=p,10 levs1=1,2,3,4,5,6,7,8,9 \ options=wedge,mirror,blacklab cols1=7 \ region=$disreg range=0.14,14,lin,2 echo 'Press return to continue'; $< ####################################### echo 'Display Stokes Q, U and V images' ####################################### foreach stk (q u v) set dev==$font'_'$stk.ps/cps set dev=/xs echo 'Display Stokes '$stk echo 'Press return to continue'; $< cgdisp in=$source.$sb.$stk'cln',$source.$sb.icln type=c,p device=$dev \ labtyp=arcsec nxy=1,1 slev=a,$qurms \ levs1=-10,-8,-6,-4,-2,2,4,6,8,10 \ options=wedge,mirror,blacklab cols1=7 \ region=$disreg range=-0.08,14,lin,-2 end echo 'Press return to continue'; $< goto MENUmap STAT: ######################################################## echo 'Make statistics from the images' ######################################################## foreach stk (i q u v) echo 'Stokes '$stk'. Press return to continue'; $< imstat in=$source.$sb.$stk'cln' region=$statreg end goto MENUmap POLVEC: ######################################################## echo ' Make polarization intensity and vector images' echo ' from the cleaned images of stokes i, q, and u' ######################################################## \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 echo 'Press return to continue'; $< set dev=$font'_'pol.ps/cps set dev=/xs cgdisp in=$source.$sb.icln,$source.$sb.pi,$source.$sb.pm,$source.$sb.pa \ type=c,pix,amp,ang device=$dev \ labtyp=arcsec options=mirror,blacklab,rot90 \ vecfac=2.0,3.0 cols1=0.5 \ range=0,0,lin,1 slev=p,10 levs1=-0.5,0.5,1,2,3,4,5,6,7,8,9 \ lines=5,8,8 region=$disreg $< cgdisp in=$source.$sb.icln,$source.$sb.pi,$source.$sb.pm,$source.$sb.pa \ type=c,pix,amp,ang device=fig1.ps/cps \ labtyp=arcsec options=mirror,blacklab,rot90 \ vecfac=2.0,3.0 cols1=0.5 \ range=0,0,lin,1 slev=p,10 levs1=-0.5,0.5,1,2,3,4,5,6,7,8,9 \ lines=5,8,8 region=$disreg echo 'Press return to continue'; $< #–––––––––– exit