|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ########################################################################## 00002 # # 00003 # 2008 NRAO Synthesis Summer School # 00004 # VLA Continuum Polarimetry Tutorial for CASA # 00005 # This script reads VLA archive files and stores them as a CASA # 00006 # measurement set, and then processes it. # 00007 # # 00008 # Original version GMoellenbrok 2008-06-08 Patch 2 summer school # 00009 # Modified/merge STM 2008-07-31 Clean up # 00010 # Updated version STM 2008-08-07 finishing touches # 00011 # Updated version GMoellenbrock2009-11-20 Add evlabands=F # 00012 # # 00013 # Based on 2008 Synthesis Imaging Workshop script # 00014 # # 00015 # Script Notes: # 00016 # o Uses the VLA export files AT166_1 to 3, which should be in the # 00017 # working directory that casapy was started in # 00018 # # 00019 # o This script has some interactive commands, such as with # 00020 # plotxy and the viewer. If scriptmode=True, then this script # 00021 # will stop and require a carriage-return to continue at these # 00022 # points. # 00023 # # 00024 # o Sometimes cut-and-paste of a series of lines from this script # 00025 # into the casapy terminal will get garbled (usually a single # 00026 # dropped character). In this case, try fewer lines, like groups # 00027 # of 4-6. # 00028 # # 00029 # # 00030 ########################################################################## 00031 00032 import time 00033 import os 00034 import pickle 00035 00036 # 00037 # This is a script that reduces B- and C-configuration VLA continuum 00038 # polarimetry data at 5 GHz on 3C129 and calibrators, and images the 00039 # dual-config 3C129 in full polarization. 00040 # 00041 # Observations: AT166 00042 # B-array: 1994Jul25 00043 # C-array: 1994Nov03 00044 # 00045 # Sources: Science Target: 3C129, a radio galaxy 00046 # Calibrator: 0420+417, observed altenately with 3C129 00047 # Calibrator: 0518+165 (3C138), for flux density/poln p.a. 00048 # Calibrator: 0134+329 (3C48), for flux density/poln p.a. 00049 # 00050 # Obs Modes: Two 50 MHz continuum (single chan) sub-bands at 4585.1 & 4885.1 MHz 00051 # Full polarization: RR,RL,LR,LL 00052 # 00053 00054 # In the following script, the B array and C array data are separately 00055 # filled and reduced. The calibration parameters are fairly standard, 00056 # and some variations may be suggested. In general, after each calibration 00057 # opearation, a data or calibration plotting command is included, again 00058 # with suggestions on possible variations. 00059 # 00060 # This script is not written to be run automatically; doing so will likely 00061 # yield a less-than-optimal result. It is intended that students will 00062 # cut-and-paste this script one task at a time to progress through the 00063 # reduction. 00064 # The methods are run in the function-call style: task(param1=x,param2=y) 00065 # To see the range of parameters for a function, use 'inp <task>'. After 00066 # running a task (in function-call style), type 'tget <task>', and 00067 # 'inp <task>' to review the parameters that were used. 00068 00069 # if you are running it and want it to stop during interactive parts. 00070 00071 scriptmode = F 00072 00073 # Enable benchmarking? 00074 benchmarking = True 00075 00076 #===================================================================== 00077 # 00078 # Set up some useful variables 00079 00080 pathname=os.environ.get('CASAPATH').split()[0] 00081 00082 prefix='3c129.tutorial' 00083 00084 # Sets a shorthand for fixed input script/regression files 00085 scriptprefix='3c129_tutorial_regression' 00086 00087 # A few parameters used repeatedly below, are defined here: 00088 # (Be sure to reset these if you exit CASA and start it up again) 00089 msnameB='at166B.ms'; 00090 msnameC='at166C.ms'; 00091 00092 # 00093 #===================================================================== 00094 # Clean up old versions of files to be created in this script 00095 os.system('rm -rf '+prefix+'.* at166B.* at166C.* 3C129BC.*') 00096 00097 # Start timing 00098 if benchmarking: 00099 startTime=time.time() 00100 startProc=time.clock() 00101 00102 # 00103 #===================================================================== 00104 # IMPORT, FLAG, CALIBRATE, IMAGE THE B-CONFIGURATION DATA 00105 #===================================================================== 00106 # 00107 print "" 00108 print "Loading B-config data..." 00109 print "" 00110 00111 #===================================================================== 00112 # Fill B-config data at C-band (5 GHz) 00113 print "--Import (Bconfig)--" 00114 importvla(archivefiles=['AT166_1', 'AT166_2'],vis=msnameB,bandname='C',evlabands=F); 00115 00116 #===================================================================== 00117 # List a summary of the dataset in the logger 00118 print "--Listobs--" 00119 listobs(vis=msnameB); 00120 00121 #--> Note scan sequence, fields, and spectral window information, etc. 00122 00123 #===================================================================== 00124 # set flux density calibrator total intensity models 00125 # NB: these sources are resolved, so we use model images provided 00126 # by the VLA (copy them from the data repository: data/VLA/CalModels/ 00127 # or point to the location on your system) 00128 # 00129 # Location of Cal Models 00130 # e.g. for MacOSX 00131 #fluxcaldir = '/opt/casa/data/nrao/VLA/CalModels/' 00132 # or standard distro 00133 fluxcaldir = pathname + '/data/nrao/VLA/CalModels/' 00134 # or in place 00135 #fluxcaldir = './' 00136 00137 # NB: By default, the model for 0420+417 is a 1 Jy point source 00138 print "--Setjy--" 00139 # Should say: 00140 ## 0518+165 spwid= 0 [I=3.688, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00141 ## 0518+165 spwid= 1 [I=3.862, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00142 setjy(vis=msnameB,field='0518+165',modimage=fluxcaldir+'3C138_C.im',scalebychan=False,standard='Perley-Taylor 99') 00143 ## 0134+329 spwid= 0 [I=5.405, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00144 ## 0134+329 spwid= 0 [I=5.739, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00145 setjy(vis=msnameB,field='0134+329',modimage=fluxcaldir+'3C48_C.im',scalebychan=False,standard='Perley-Taylor 99') 00146 00147 #===================================================================== 00148 # Plot data and interactively edit 00149 # One spw at a time, calibrators only, RR, LL only 00150 print "--Plotxy--" 00151 if scriptmode: 00152 plotxy(vis=msnameB,spw='0',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL'); 00153 print 'Plotting 0420+417,0518+165,0134+329 SPW 0' 00154 print " Note VA01 seems to be bad (low) on one integration of one scan" 00155 print " Mark, Locate, and Flag this" 00156 user_check=raw_input('hit Return to continue script\n') 00157 00158 plotxy(vis=msnameB,spw='1',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL'); 00159 print 'Plotting 0420+417,0518+165,0134+329 SPW 1' 00160 print " Mark, Locate, and Flag bad VA01 integration in this SPW 1 also" 00161 user_check=raw_input('hit Return to continue script\n') 00162 else: 00163 plotxy(vis=msnameB,spw='0',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL',interactive=F,figfile='at166B.plotxy.initial.spw0.png'); 00164 plotxy(vis=msnameB,spw='1',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL',interactive=F,figfile='at166B.plotxy.initial.spw1.png'); 00165 print "--Flagdata--" 00166 tflagdata(vis=msnameB,antenna='VA01',timerange='1994/07/25/14:21:10.0~14:21:20.0',mode='manual') 00167 00168 #--> Variations: plot on different axes (e.g., phase vs. time, real vs imag, amp/phase vs. uvdist) 00169 #--> plot model data for calibrators 00170 #--> plot 3C129 data 00171 00172 #===================================================================== 00173 # Solve for gains on calibrators 00174 # NB: reference phases to VA12; pre-apply parallactic angle correction 00175 print "--Gaincal--" 00176 gaincal(vis=msnameB,caltable='at166B.gcal',field='0420+417,0518+165,0134+329',refant='VA12',parang=T); 00177 00178 #--> (22/22 good solutions) 00179 00180 #===================================================================== 00181 # Examine solutions: 00182 print "--Plotcal--" 00183 # Use subplot panels to stack amp, phase, snr above each other 00184 if scriptmode: 00185 plotcal(caltable='at166B.gcal',yaxis='amp',subplot=311); 00186 plotcal(caltable='at166B.gcal',yaxis='phase',subplot=312); 00187 plotcal(caltable='at166B.gcal',yaxis='snr',subplot=313); 00188 user_check=raw_input('hit Return to continue script\n') 00189 else: 00190 plotcal(caltable='at166B.gcal',yaxis='amp',subplot=211,showgui=F,figfile=''); 00191 plotcal(caltable='at166B.gcal',yaxis='phase',subplot=212,showgui=F,figfile='at166B.gcal.plotcal.png'); 00192 00193 #--> Variations: plot solutions per antenna using iteration='antenna' 00194 00195 #===================================================================== 00196 # Scale gain solution from 0420+417 according to f.d. calibrators 00197 print "--Fluxscale--" 00198 fluxscale(vis=msnameB,caltable='at166B.gcal',fluxtable='at166B.fcal',reference='0518+165,0134+329'); 00199 00200 # --> 0420: 1.441/1.443 Jy (broadly consistent with VLA Cal man values at L & X) 00201 00202 #===================================================================== 00203 # Examine solutions: 00204 print "--Plotcal--" 00205 if scriptmode: 00206 plotcal(caltable='at166B.fcal',yaxis='amp'); 00207 user_check=raw_input('hit Return to continue script\n') 00208 else: 00209 plotcal(caltable='at166B.fcal',yaxis='amp',showgui=F,figfile='at166B.fcal.plotcal.png'); 00210 00211 #--> Note that gain amps are ~constant now 00212 00213 #===================================================================== 00214 # Solve for instrumental polarization on 0420+417 (and also for source poln) 00215 # NB: IMPORTANT: use gcal---not fcal---here because model is _still_ 1.0 Jy 00216 print "--Polcal (D)--" 00217 polcal(vis=msnameB,caltable='at166B.dcal',field='0420+417',refant='VA12',gaintable='at166B.gcal',gainfield='0420+417'); 00218 00219 # --> 2/2 good solutions; 00220 # --> 0420: 0.027Jy @ 61.2deg / 0.024Jy @ -83.0 Jy 00221 00222 #===================================================================== 00223 # Examine solutions: 00224 print "--Plotcal--" 00225 if scriptmode: 00226 plotcal(caltable='at166B.dcal',xaxis='antenna',yaxis='amp'); 00227 user_check=raw_input('hit Return to continue script\n') 00228 else: 00229 plotcal(caltable='at166B.dcal',xaxis='antenna',yaxis='amp',showgui=F,figfile='at166B.dcal.plotcal.png'); 00230 00231 # You can plot imag vs. real also 00232 # plotcal(caltable='at166B.dcal',xaxis='real',yaxis='imag',plotrange=[-0.05,0.05,-0.05,0.05]); 00233 00234 #===================================================================== 00235 # Set full polarization model for 0518+165 (pol is 11.1% @ -11 deg [RL = -22]) 00236 # NB: neglecting source structure here) 00237 print "--Setjy (X)--" 00238 setjy(vis=msnameB,field='0518+165',spw='0',scalebychan=False,fluxdensity=[3.688, 0.380, -0.153, 0.0]); 00239 setjy(vis=msnameB,field='0518+165',spw='1',scalebychan=False,fluxdensity=[3.862, 0.397, -0.161, 0.0]); 00240 00241 #===================================================================== 00242 # Solve for polarization position angle on 0518+165 00243 print "--Polcal (X)--" 00244 polcal(vis='at166B.ms',caltable='at166B.xcal',field='0518+165',refant='VA12',poltype='X',gaintable=['at166B.fcal', 'at166B.dcal'],gainfield=['0518+165','0420+417']); 00245 00246 # --> 2/2 good solutions 00247 # --> 0518: 2.7deg / 41.8deg (I need to check these numbers) 00248 00249 #===================================================================== 00250 # apply all calibration... 00251 # (NB: different fields are selected from each caltable, depending on selected data fields) 00252 # (NB: parang=T is set to rotate polarization p.a. frame from antennas to sky 00253 print "--Applycal--" 00254 print " apply calibration to 0420+417,3C129" 00255 applycal(vis=msnameB,field='0420+417,3C129',gaintable=['at166B.fcal','at166B.dcal','at166B.xcal'],gainfield=['0420+417','0420+417','0518+165'],parang=T); 00256 00257 print " apply calibration to 0518+165" 00258 applycal(vis=msnameB,field='0518+165',gaintable=['at166B.fcal','at166B.dcal','at166B.xcal'],gainfield=['0518+165','0420+417','0518+165'],parang=T); 00259 00260 print " apply calibration to 0134+329" 00261 applycal(vis=msnameB,field='0134+329',gaintable=['at166B.fcal','at166B.dcal','at166B.xcal'],gainfield=['0134+329','0420+417','0518+165'],parang=T); 00262 00263 #===================================================================== 00264 # Examine (edit?) calibrated data (calibrators) 00265 print "--Plotxy (corrected)--" 00266 if scriptmode: 00267 plotxy(vis=msnameB,datacolumn='corrected',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL'); 00268 user_check=raw_input('hit Return to continue script\n') 00269 else: 00270 plotxy(vis=msnameB,datacolumn='corrected',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL',interactive=F,figfile='at166B.plotxy.final.png'); 00271 00272 #--> Variations: plot on different axes (e.g., phase vs. time, real vs imag) 00273 #--> overlay model data for calibrators 00274 #--> plot 3C129 data 00275 # For example: 00276 # Examine cross-hand data (real vs. imag) 00277 # plotxy(vis=msnameB,xaxis='real',yaxis='imag',datacolumn='corrected',selectdata=True,correlation='RL,LR',plotrange=[-0.5,0.5,-0.5,0.5],field = '0518+165'); 00278 # plotxy(vis=msnameB,xaxis='real',yaxis='imag',datacolumn='corrected',selectdata=True,correlation='RL,LR',plotrange=[-0.5,0.5,-0.5,0.5],field = '0134+329'); 00279 # plotxy(vis=msnameB,xaxis='real',yaxis='imag',datacolumn='corrected',selectdata=True,correlation='RL,LR',plotrange=[-0.5,0.5,-0.5,0.5],field = '0420+417'); 00280 00281 #--> NB: RL and LR signal are complex conjugates of each other (Q+iU & Q-iU) 00282 00283 #===================================================================== 00284 # do some simple imaging of each source 00285 print "--Clean--" 00286 # 3C129: 00287 # We will do a simple image-plane Hogbom clean with psfmode='hogbom' and imagermode='' 00288 # This will clean the IQUV planes consecutively 00289 # 00290 # Non-interactive (no clean boxes) for now 00291 # 00292 clean(vis=msnameB,imagename='at166B.3c129',field='3C129',psfmode='hogbom',imagermode='',niter=2500,imsize=[2048,2048],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F); 00293 00294 print " Clean image is at166B.3c129.image" 00295 00296 # You can do a a Cotton-Schwab clean with psfmode='clark' and imagermode='csclean' 00297 # You can try a threshold also. 00298 00299 # You can clean the calibrators also: 00300 # 0518: 00301 # clean(vis=msnameB,imagename='at166B.0518',field='0518+165',psfmode='hogbom',niter=500,imsize=[512,512],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F); 00302 00303 # 0134: 00304 # clean(vis=msnameB,imagename='at166B.0134',field='0134+329',psfmode='hogbom',niter=500,imsize=[512,512],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F); 00305 00306 # 0420: 00307 # clean(vis=msnameB,imagename='at166B.0420',field='0420+417',psfmode='hogbom',niter=500,imsize=[512,512],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F); 00308 00309 #===================================================================== 00310 # Examine images in viewer... 00311 if scriptmode: 00312 print "--Viewer--" 00313 viewer('at166B.3c129.image') 00314 user_check=raw_input('When done, close viewer and hit Return to continue script\n') 00315 00316 # Questions: Do calibrator polarizations come out right? 00317 # Is calibrator structure as expected? 00318 # Is 3C129 image any good? (cf C-config imaging below) 00319 00320 # Variations: clean interactively (clean boxes) 00321 # try different weighting 00322 00323 #===================================================================== 00324 # IMPORT, FLAG, CALIBRATE, IMAGE THE C-CONFIGURATION DATA 00325 #===================================================================== 00326 # NOW, reduce C-config data in much the same way... 00327 print "" 00328 print "Now reducing C-config data..." 00329 print "" 00330 00331 # Fill C-config data at C-band (5 GHz) 00332 print "--Import (Cconfig)--" 00333 importvla(archivefiles='AT166_3',vis=msnameC,bandname='C',evlabands=F); 00334 00335 #===================================================================== 00336 # List a summary of the dataset in the logger 00337 print "--Listobs--" 00338 listobs(vis=msnameC); 00339 00340 #===================================================================== 00341 # set flux density calibrator total intensity models 00342 # NB: these sources are resolved, so we use model images provided 00343 # by the VLA, with location set in B-config part of script 00344 # NB: By default, the model for 0420+417 is a 1 Jy point source 00345 print "--Setjy--" 00346 setjy(vis=msnameC,field='0518+165',modimage=fluxcaldir+'3C138_C.im',scalebychan=False,standard='Perley-Taylor 99'); 00347 setjy(vis=msnameC,field='0134+329',modimage=fluxcaldir+'3C48_C.im',scalebychan=False,standard='Perley-Taylor 99'); 00348 00349 #===================================================================== 00350 # Plot data and interactively edit 00351 # One spw at a time, calibrators only, RR, LL only 00352 print "--Plotxy--" 00353 if scriptmode: 00354 plotxy(vis=msnameC,spw='0',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL'); 00355 print 'Plotting 0420+417,0518+165,0134+329 SPW 0' 00356 user_check=raw_input('hit Return to continue script\n') 00357 00358 plotxy(vis=msnameC,spw='1',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL'); 00359 print 'Plotting 0420+417,0518+165,0134+329 SPW 1' 00360 user_check=raw_input('hit Return to continue script\n') 00361 else: 00362 plotxy(vis=msnameC,spw='0',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL',interactive=F,figfile='at166C.plotxy.initial.spw0.png'); 00363 plotxy(vis=msnameC,spw='1',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL',interactive=F,figfile='at166C.plotxy.initial.spw1.png'); 00364 00365 #--> Variations: plot on different axes (e.g., phase vs. time, real vs imag, amp/phase vs. uvdist) 00366 #--> plot model data for calibrators (datacolumn='model') 00367 #--> plot 3C129 data 00368 #--> plot individual baselines (iteration='baseline') 00369 00370 # Solve for gains on calibrators 00371 # NB: reference phases to VA12; pre-apply parallactic angle correction 00372 print "--Gaincal--" 00373 gaincal(vis=msnameC,caltable='at166C.gcal',field='0420+417,0518+165,0134+329',refant='VA12',parang=T,solint=500); 00374 00375 #--> (12/14 good solutions) NB: no solutions for 0134+329 (too little data?) 00376 00377 #===================================================================== 00378 # Examine solutions: 00379 print "--Plotcal--" 00380 if scriptmode: 00381 plotcal(caltable='at166C.gcal',yaxis='amp'); 00382 user_check=raw_input('hit Return to continue script\n') 00383 00384 plotcal(caltable='at166C.gcal',yaxis='phase'); 00385 user_check=raw_input('hit Return to continue script\n') 00386 else: 00387 plotcal(caltable='at166C.gcal',yaxis='amp',subplot=211,showgui=F,figfile=''); 00388 plotcal(caltable='at166C.gcal',yaxis='phase',subplot=212,showgui=F,figfile='at166C.gcal.plotcal.png'); 00389 00390 #--> Variations: plot solutions per antenna, spw, etc. 00391 00392 #===================================================================== 00393 # Scale gain solution from 0420+417 according to f.d. calibrators 00394 print "--Fluxscale--" 00395 fluxscale(vis=msnameC,caltable='at166C.gcal',fluxtable='at166C.fcal',reference='0518+165'); 00396 # --> 0420: 1.282/1.292 Jy (broadly consistent with VLA Cal man values at L & X) 00397 00398 #===================================================================== 00399 # Examine solutions: 00400 print "--Plotcal--" 00401 if scriptmode: 00402 plotcal(caltable='at166C.fcal',yaxis='amp'); 00403 user_check=raw_input('hit Return to continue script\n') 00404 else: 00405 plotcal(caltable='at166C.fcal',yaxis='amp',showgui=F,figfile='at166C.fcal.plotcal.png'); 00406 00407 #--> gain amplitudes now ~constant 00408 00409 # Solve for instrumental polarization on 0420+417 (and also for source poln) 00410 # NB: IMPORTANT: use gcal---not fcal---here because model is _still_ 1.0 Jy 00411 print "--Polcal (D)--" 00412 polcal(vis=msnameC,caltable='at166C.dcal',field='0420+417',refant='VA12',gaintable='at166C.gcal',gainfield='0420+417'); 00413 # --> 2/2 good solutions; 00414 # --> 0420: 0.035Jy @ -41.4deg / 0.033Jy @ 18.1deg 00415 00416 #===================================================================== 00417 # Examine solutions: 00418 print "--Plotcal--" 00419 if scriptmode: 00420 plotcal(caltable='at166C.dcal',xaxis='antenna',yaxis='amp'); 00421 user_check=raw_input('hit Return to continue script\n') 00422 else: 00423 plotcal(caltable='at166C.dcal',xaxis='antenna',yaxis='amp',showgui=F,figfile='at166C.dcal.plotcal.png'); 00424 00425 # You can plot imag vs. real also 00426 # plotcal(caltable='at166C.dcal',xaxis='real',yaxis='imag',plotrange=[-0.05,0.05,-0.05,0.05]); 00427 00428 #===================================================================== 00429 # Set full polarization model for 0518+165 (pol is 11.1% @ -11 deg [RL = -22]) 00430 # NB: neglecting source structure here) 00431 print "--Setjy (X)--" 00432 setjy(vis=msnameC,field='0518+165',spw='0',scalebychan=False,fluxdensity=[3.688, 0.380, -0.153, 0.0]); 00433 setjy(vis=msnameC,field='0518+165',spw='1',scalebychan=False,fluxdensity=[3.862, 0.397, -0.161, 0.0]); 00434 00435 #===================================================================== 00436 # Solve for polarization position angle on 0518+165 00437 print "--Polcal (X)--" 00438 polcal(vis='at166C.ms',caltable='at166C.xcal',field='0518+165',refant='VA12',poltype='X',gaintable=['at166C.fcal', 'at166C.dcal'],gainfield=['0518+165','0420+417']); 00439 # --> 2/2 good solutions 00440 # --> 0518: 77.0deg / -42.0deg 00441 00442 #===================================================================== 00443 # apply all calibration... 00444 # (NB: different fields are selected from each caltable, depending on selected data fields) 00445 # (NB: parang=T is set to rotate polarization p.a. frame from antennas to sky 00446 #...to 0420 & 3C129 00447 print "--Applycal--" 00448 print " apply calibration to 0420+417,3C129" 00449 applycal(vis=msnameC,field='0420+417,3C129',gaintable=['at166C.fcal','at166C.dcal','at166C.xcal'],gainfield=['0420+417','0420+417','0518+165'],parang=T); 00450 00451 print " apply calibration to 0518+165" 00452 applycal(vis=msnameC,field='0518+165',gaintable=['at166C.fcal','at166C.dcal','at166C.xcal'],gainfield=['0518+165','0420+417','0518+165'],parang=T); 00453 00454 #===================================================================== 00455 # Examine (edit?) calibrated data (calibrators) 00456 print "--Plotxy (corrected)--" 00457 if scriptmode: 00458 plotxy(vis=msnameC,datacolumn='corrected',field='0420+417,0518+165',selectdata=T,correlation='RR,LL'); 00459 user_check=raw_input('hit Return to continue script\n') 00460 else: 00461 plotxy(vis=msnameC,datacolumn='corrected',field='0420+417,0518+165',selectdata=T,correlation='RR,LL',interactive=F,figfile='at166C.plotxy.final.png'); 00462 00463 #--> Variations: plot on different axes (e.g., phase vs. time, real vs imag) 00464 #--> overlay model data for calibrators 00465 #--> plot 3C129 data 00466 # For example: 00467 # Examine cross-hand data (real vs. imag) 00468 # plotxy(vis=msnameC,xaxis='real',yaxis='imag',datacolumn='corrected',selectdata=True,correlation='RL,LR',plotrange=[-0.5,0.5,-0.5,0.5],field = '0518+165'); 00469 # plotxy(vis=msnameC,xaxis='real',yaxis='imag',datacolumn='corrected',selectdata=True,correlation='RL,LR',plotrange=[-0.5,0.5,-0.5,0.5],field = '0420+417'); 00470 00471 #--> NB: RL and LR signal are complex conjugates of each other (Q+iU & Q-iU) 00472 00473 #===================================================================== 00474 # do some simple imaging of each source 00475 print "--Clean--" 00476 # 3C129: 00477 # We will do a simple image-plane Hogbom clean with psfmode='hogbom' and imagermode='' 00478 # This will clean the IQUV planes consecutively 00479 # 00480 # Non-interactive (no clean boxes) for now 00481 # 00482 clean(vis=msnameC,imagename='at166C.3c129',field='3C129',psfmode='hogbom',imagermode='',niter=2500,imsize=[2048,2048],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F); 00483 00484 print " Clean image is at166C.3c129.image" 00485 00486 # You can do a a Cotton-Schwab clean with psfmode='clark' and imagermode='csclean' 00487 # You can try a threshold also. 00488 00489 # You can also clean the calibrators: 00490 # 0518: 00491 # clean(vis=msnameC,imagename='at166C.0518',field='0518+165',psfmode='hogbom',niter=500,imsize=[512,512],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F); 00492 00493 # 0420: 00494 # clean(vis=msnameC,imagename='at166C.0420',field='0420+417',psfmode='hogbom',niter=500,imsize=[512,512],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F); 00495 00496 #===================================================================== 00497 # Examine images in viewer... 00498 if scriptmode: 00499 print "--Viewer--" 00500 viewer('at166C.3c129.image') 00501 print " Notice better representation of large-scale emission than in B-config" 00502 user_check=raw_input('hit Return to continue script\n') 00503 00504 # Questions: Do calibrator polarizations come out right? 00505 # Is calibrator structure as expected? 00506 # Is 3C129 image good? (cf B-config imaging above) 00507 00508 # Variations: clean interactively (clean boxes) 00509 # try different weighting 00510 # in viewer, blink between B and C 3C129 images (in all polarizations) 00511 00512 00513 #===================================================================== 00514 # IMAGING OF COMBINED B+C CONFIGURATIONS 00515 #===================================================================== 00516 print "" 00517 print "Combining B and C config data..." 00518 print "" 00519 00520 # The next steps extract the 3C129 data from the above datasets, and combine 00521 # them to permit a dual-config imaging run 00522 00523 #===================================================================== 00524 # split out 3C129 00525 print "--Split (3C129)--" 00526 split(vis=msnameB,outputvis='at166B.3C129.ms',field='3C129'); 00527 split(vis=msnameC,outputvis='at166C.3C129.ms',field='3C129'); 00528 00529 #===================================================================== 00530 # make one MS so we can image the combined config 00531 print "--Concat (B+C config)--" 00532 concat(vis=['at166B.3C129.ms','at166C.3C129.ms'],concatvis='3C129BC.ms'); 00533 00534 #===================================================================== 00535 print "--Listobs--" 00536 listobs(vis='3C129BC.ms'); 00537 00538 #===================================================================== 00539 # Clean the image 00540 print "--Clean--" 00541 # 00542 # You can do a simple image-plane Hogbom clean with psfmode='hogbom' and imagermode='' 00543 # This will clean the IQUV planes consecutively 00544 # 00545 # We will do a Cotton-Schwab clean with psfmode='clark' and imagermode='csclean' 00546 # This will do a joint IQUV deconvolution 00547 # We will also set a threshold 00548 00549 if scriptmode: 00550 print "Use interactive clean to draw clean boxes or polygons on image" 00551 print "Increase npercycle as you clean deeper" 00552 clean(vis='3C129BC.ms',imagename='3C129BC.clean',psfmode='clark',imagermode='csclean',niter=50000,threshold='0.08mJy',imsize=[2048,2048],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=T,npercycle=1000) 00553 else: 00554 clean(vis='3C129BC.ms',imagename='3C129BC.clean',psfmode='clark',imagermode='csclean',niter=50000,threshold='0.16mJy',imsize=[2048,2048],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F); 00555 00556 # The non-interactive clean gives a S/N ratio around 850. 00557 # If you do custom clean boxes and clean more deeply you should 00558 # be able to reach a dynamic range over 1000. 00559 00560 print " Clean image is 3C129BC.clean.image" 00561 00562 # Variation: you could set up clean boxes around emission, e.g. 00563 # mask = [[992,988,1172,1139],[1029,1053,1212,1222],[1029,1081,1299,1222],[1065,1158,1397,1409],[1237,1153,1612,1583]] 00564 00565 #===================================================================== 00566 # Examine image in viewer... 00567 if scriptmode: 00568 print "--Viewer--" 00569 viewer('3C129BC.clean.image') 00570 00571 # Questions: How has the 3C129 image changed? Is it better? 00572 00573 # Variations: clean interactively (clean boxes) 00574 # try different weighting 00575 # split/concat/image dual-config calibrator data 00576 # self-cal? 00577 00578 #===================================================================== 00579 # Extract I, Q, U, V images 00580 print "--Immath--" 00581 immath(stokes='I', outfile='3C129BC.I', mode='evalexpr', expr='\'3C129BC.clean.image\''); 00582 immath(stokes='Q', outfile='3C129BC.Q', mode='evalexpr', expr='\'3C129BC.clean.image\''); 00583 immath(stokes='U', outfile='3C129BC.U', mode='evalexpr', expr='\'3C129BC.clean.image\''); 00584 immath(stokes='V', outfile='3C129BC.V', mode='evalexpr', expr='\'3C129BC.clean.image\''); 00585 00586 # Form poln intensity and pos ang 00587 immath(stokes='', outfile='3C129BC.P', mode='poli', imagename=['3C129BC.Q','3C129BC.U'], sigma='0.0mJy/beam'); 00588 immath(stokes='', outfile='3C129BC.X', mode='pola', imagename=['3C129BC.Q','3C129BC.U'], sigma='0.0mJy/beam'); 00589 00590 00591 #===================================================================== 00592 # Extract center of I image for testing 00593 00594 ia.open('3C129BC.I') 00595 ia.subimage(outfile='3C129BC.core.I',region=rg.box([1010,1040,0,0],[1025,1055,0,0])) 00596 ia.close() 00597 00598 #===================================================================== 00599 # Complex Linear Polarization 00600 print "--ComplexLinPol (toolkit)--" 00601 00602 # NOTE: The viewer can take complex images to make Vector plots, although 00603 # the image analysis tasks (and ia tool) cannot yet handle these. But we 00604 # can use the imagepol tool (which is not imported by default) to make 00605 # a complex image of the linear polarized intensity for display. 00606 # See CASA User Reference Manual: 00607 # http://casa.nrao.edu/docs/casaref/imagepol-Tool.html 00608 # 00609 # Make an imagepol tool and open the clean image 00610 po = casac.imagepol() 00611 po.open('3C129BC.clean.image') 00612 # Use complexlinpol to make a Q+iU image 00613 po.complexlinpol('3C129BC.cmplxlinpol') 00614 po.close() 00615 00616 # You can now display this in the viewer, in particular overlay this 00617 # over the intensity raster with the poli contours. The vector lengths 00618 # will be proportional to the polarized intensity. You can play with 00619 # the Data Display Options panel for vector spacing and length. 00620 # You will want to have this masked, like the pola image above, on 00621 # the polarized intensity. When you load the image, use the LEL: 00622 # '3C129BC.cmplxlinpol'['3C129BC.P'>0.0001] 00623 00624 #===================================================================== 00625 # View results 00626 # In viewer 'Load Data' window, use the following LEL expression to load the p.a. image: 00627 print "--Viewer--" 00628 if scriptmode: 00629 viewer('3C129BC.I') 00630 else: 00631 print " For viewer: " 00632 00633 print " In Load Data' window LEL box, use the following: " 00634 print " '3C129BC.X'['3C129BC.P'>0.00015] " 00635 print " or" 00636 print " '3C129BC.cmplxlinpol'['3C129BC.P'>0.0001] " 00637 00638 if scriptmode: 00639 user_check=raw_input('hit Return to continue script\n') 00640 00641 #===================================================================== 00642 # RESULTS 00643 #===================================================================== 00644 # Calculate image statistics 00645 # Outer source-free region 1100,100,1950,750 00646 outerbox = '1100,100,1950,750' 00647 # Inner source-free region 1315,940,1532,1143 (in clean bowl area) 00648 innerbox = '1315,940,1532,1143' 00649 print "--Imstat--" 00650 00651 # Stats on IPOL 00652 ipolstat = imstat('3C129BC.clean.image',stokes='I') 00653 ipolstat_offsrc = imstat('3C129BC.clean.image',stokes='I',box=outerbox) 00654 00655 print ' %40s : %12.7f ' % ('3C129 Combined I maximum (Jy)', ipolstat['max'][0]) 00656 print ' %40s : %12.7f ' % ('3C129 Combined I off-source rms (Jy)', ipolstat_offsrc['sigma'][0]) 00657 print ' %40s : %12.3f ' % ('3C129 Combined I dynamic range', 00658 ipolstat['max'][0]/ipolstat_offsrc['sigma'][0]) 00659 00660 # Stats on QPOL,UPOL 00661 qupolstat = imstat('3C129BC.clean.image',stokes='QU') 00662 qupolstat_offsrc = imstat('3C129BC.clean.image',stokes='QU',box=outerbox) 00663 quabsmax = max( qupolstat['max'][0], -qupolstat['min'][0] ) 00664 00665 print '' 00666 print ' %40s : %12.7f ' % ('3C129 Combined QU maximum (Jy)', quabsmax) 00667 print ' %40s : %12.7f ' % ('3C129 Combined QU off-source rms (Jy)', qupolstat_offsrc['sigma'][0]) 00668 print ' %40s : %12.3f ' % ('3C129 Combined QU dynamic range', 00669 quabsmax/qupolstat_offsrc['sigma'][0]) 00670 00671 # Stats on entire V image 00672 vpolstat = imstat('3C129BC.clean.image',stokes='V') 00673 00674 print '' 00675 print ' %40s : %12.7f ' % ('3C129 Combined VPOL image rms (Jy)', vpolstat['sigma'][0]) 00676 00677 #===================================================================== 00678 # DONE 00679 #===================================================================== 00680 if benchmarking: 00681 endProc=time.clock() 00682 endTime=time.time() 00683 00684 print "" 00685 print "Done with 3C129 Tutorial" 00686 00687 if benchmarking: 00688 walltime = (endTime - startTime) 00689 cputime = (endProc - startProc) 00690 00691 print '' 00692 print 'Regression PASSED' 00693 print '' 00694 print 'Total wall clock time was: %10.3f ' % walltime 00695 print 'Total CPU time was: %10.3f ' % cputime 00696 print '' 00697 00698 # 00699 ##########################################################################
1.8.0
