|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ########################################################################## 00002 # # 00003 # Use Case Script for POLCAL 6cm Data # 00004 # Using POLCA data 20080224 BnC-config C-band # 00005 # # 00006 # Last Updated STM 2008-05-23 (Beta Patch 2) # 00007 # Updated STM 2008-06-11 (Beta Patch 2.0) # 00008 # Uses new clean task # 00009 # # 00010 ########################################################################## 00011 00012 import time 00013 import os 00014 import pickle 00015 00016 # 00017 #===================================================================== 00018 # 00019 # This script may have some interactive commands: scriptmode = True 00020 # if you are running it and want it to stop during interactive parts. 00021 00022 scriptmode = True 00023 00024 #===================================================================== 00025 # SET UP THE SCRIPT CONTROL PARAMETERS HERE 00026 #===================================================================== 00027 # 00028 # Set up some useful variables to control subsequent actions: 00029 00030 pathname=os.environ.get('CASAPATH').split()[0] 00031 00032 # This name will prefix all output files 00033 prefix = 'polcal_20080224.cband.all' 00034 00035 # Clean up old files 00036 os.system('rm -rf '+prefix+'*') 00037 00038 #===================================================================== 00039 00040 # Import data from export or use already existing MS? Or UVFITS? 00041 importmode = 'vla' # 'vla','fits','ms' 00042 # This is the name of the datafile used in import 00043 # or the name of a previously made ms that will be copied 00044 # NOTE: if an ms name must be different than prefix + '.ms' 00045 #datafile = 'polcal_20080224.cband.edited.ms' 00046 #datafile = '20080224C.UVF' 00047 # 00048 # NOTE: This file may be obtained from the CASA repository: 00049 # http://casa.nrao.edu/Data/VLA/Polcal/POLCA_20080224_1 00050 datafile = ['POLCA_20080224_1'] 00051 # 00052 # If from export set these: 00053 exportproject = 'POLCA' 00054 exportband = 'C' 00055 # 00056 # Spectral windows to use in ms (usually 0,1) 00057 usespw = '' 00058 usespwlist = ['0','1'] 00059 00060 # The ms will have this name 00061 msfile = prefix + '.ms' 00062 00063 # These are names of calibration tables 00064 gtable = prefix + '.gcal' 00065 ftable = prefix + '.fluxscale' 00066 ptable = prefix + '.pcal' 00067 xtable = prefix + '.polx' 00068 00069 # Flagging: 00070 myquackinterval = 14.0 # if >0 then quack scan beginnings 00071 # Flagging these antennas (if blank then no flagging) 00072 # NOTE: This script uses NEW names, so VLA ants are VAxx 00073 flagants = '' 00074 #flagants = 'EA*' # keep only VLA antennas 00075 #flagants = 'VA*' # keep only EVLA antennas 00076 00077 # 00078 # List of sources in ms 00079 # 00080 # 0 A 1924-292 19:24:51.06 -29.14.30.12 J2000 00081 # 1 A 1743-038 17:43:58.86 -03.50.04.62 J2000 00082 # 2 A 2202+422 22:02:43.29 +42.16.39.98 J2000 00083 # 3 A 2253+161 22:53:57.75 +16.08.53.56 J2000 00084 # 4 B 2136+006 21:36:38.59 +00.41.54.21 J2000 00085 # 5 B 0137+331 01:37:41.30 +33.09.35.13 J2000 00086 # 6 A 2355+498 23:55:09.46 +49.50.08.34 J2000 00087 # 7 B 0319+415 03:19:48.16 +41.30.42.10 J2000 00088 # 8 B 0359+509 03:59:29.75 +50.57.50.16 J2000 00089 # 00090 # These sources are the gain calibrators 00091 gaincalfield = ['0137+331','2202+422','1743-038','1924-292','2136+006','2253+161','2355+498','0319+415','0359+509'] 00092 # 00093 # These sources will have calibration transferred from srclist 00094 targets = [] 00095 00096 # Assemble field strings from lists 00097 fieldgain = '' 00098 if ( len(gaincalfield) > 0 ): 00099 for fn in range(len(gaincalfield)): 00100 if ( fn > 0 ): 00101 fieldgain += ',' 00102 fieldgain += gaincalfield[fn] 00103 00104 fieldtargets = '' 00105 if ( len(targets) > 0 ): 00106 for fn in range(len(targets)): 00107 if ( fn > 0 ): 00108 fieldtargets += ',' 00109 fieldtargets += targets[fn] 00110 00111 # 00112 # This list is used for final clean and stats 00113 srclist = gaincalfield + targets 00114 00115 # Location of Cal Models 00116 # e.g. for MacOSX 00117 #fluxcaldir = '/opt/casa/data/nrao/VLA/CalModels/' 00118 # or standard distro 00119 fluxcaldir = pathname + '/data/nrao/VLA/CalModels/' 00120 # or in place 00121 #fluxcaldir = './' 00122 00123 # Calibration parameters: 00124 fluxcalfield = '0137+331' # primary calibrator for setjy 00125 fluxcalmodel = '3C48_C.im' # if non-blank use this model image 00126 gaincalfield = '' # names of gain calibrators (''=all fields) 00127 usegaincurve = False # use a-priori antenna gain-elevation curve? 00128 gainopacity = 0.0 # a-priori atmospheric optical depth (Tau) 00129 calrefant = 'VA15' # reference antenna name for calibration (VA15,EA19) 00130 gainsolint = 20.0 # 20s for gaincal solutions 00131 polcalfield = '2202+422' # polarization (D-term) calibrator 00132 polcalmode = 'D+QU' # polarization (D-term) calibration mode 00133 polduvrange = '' # uvrange for polcal D 00134 setpolmodel = True # if true then use setjy to set pol model 00135 polxfield = '0137+331' # polarization angle (X) calibrator 00136 polxuvrange = '' # uvrange for polcal X 00137 # 00138 setjymode = 'set' # mode for fluxcal setyjy: 'set', 'flux', 'ft' 00139 00140 # This is the name of the split file for corrected data 00141 srcsplitms = prefix + '.split.ms' 00142 00143 # 00144 # Set up general clean parameters 00145 00146 # This is BnC-config VLA 6cm (4.85GHz) obs 00147 # Check the observational status summary 00148 # Primary beam FWHM = 45'/f_GHz = 557" 00149 # Synthesized beam for VLA/EVLA at C-Band: 00150 # A-config FWHM = 0.4" 00151 # B-config FWHM = 1.2" 00152 # C-config FWHM = 3.9" 00153 # D-config FWHM = 14.0" 00154 # RMS in 10min (600s) = 0.06 mJy (thats now, but close enough) 00155 # 00156 # Set the output image size and cell size (arcsec) 00157 # 0.4" will give 3x oversampling at least 00158 # clean will say to use a composite integer (e.g.288) for efficiency 00159 #clnalg = 'clark' 00160 clnalg = 'hogbom' 00161 usecsclean = False 00162 clnimsize = 288 00163 clncell = 0.4 00164 # Fix maximum number of iterations 00165 clniter = 200 00166 # Also set flux residual threshold (0.04 mJy) 00167 # Our scans are around 120s 00168 # With rms of 0.06 mJy in 600s ==> rms = 0.13 mJy 00169 # Set to 10x thermal rms 00170 clthreshold = 1.3 00171 00172 # Set up a clean box in the center (1/8 of image) 00173 clncenter = clnimsize/2 00174 clnblc = clncenter - clnimsize/8 00175 clntrc = clncenter + clnimsize/8 00176 # For poor uv coverage, use tigher box (6 x SynthBeam = 18xcell) 00177 clnblc = clncenter - 10 00178 clntrc = clncenter + 10 00179 centerbox = [clnblc,clnblc,clntrc,clntrc] 00180 00181 myclnbox = centerbox 00182 # Can also force interactive cleaning 00183 #myclnbox = 'interactive' 00184 00185 aipsmodel = {} 00186 # 00187 #===================================================================== 00188 # 00189 # Polarization of X angle calibrator 0137+331 00190 # If setpolmodel = True 00191 # 00192 # Set up fluxcalmodel 00193 # 00194 fcalmodel = {} 00195 # 00196 # The flux model for 0137+331 (C-band) 00197 fcalfield = {} 00198 # NOTE: you must have entries for all spw in usespwlist 00199 # I,Q,U,V 00200 fcalfield['0'] = [5.405,0,0,0] 00201 fcalfield['1'] = [5.458,0,0,0] 00202 fcalmodel['0137+331'] = fcalfield 00203 # Put in 2202+422 00204 # These values from AIPS (http://www.vla.nrao.edu/astro/calib/polar/2004/) 00205 fcalfield = {} 00206 fcalfield['0'] = [2.465,0,0,0] 00207 fcalfield['1'] = [2.461,0,0,0] 00208 fcalmodel['2202+422'] = fcalfield 00209 # 00210 # Set up pcalmodel 00211 # 00212 pcalmodel = {} 00213 # 00214 # The polarization model for 0137+331 00215 pcalfield = {} 00216 # NOTE: you must have entries for all spw in usespwlist 00217 # From calibrator manual: C-band RLPD=-148deg P/I=0.041 00218 # IPOL,FPOL,RLPHASE 00219 pcalfield['0'] = [5.405,0.041,-148.0] 00220 pcalfield['1'] = [5.458,0.041,-148.0] 00221 pcalmodel['0137+331'] = pcalfield 00222 # Put in 2202+422 (with effective flux of 1.0 before fluxscale) 00223 # These values from AIPS (http://www.vla.nrao.edu/astro/calib/polar/2004/) 00224 pcalfield = {} 00225 pcalfield['0'] = [1.0,0.072,-55.00] 00226 pcalfield['1'] = [1.0,0.072,-55.00] 00227 pcalmodel['2202+422'] = pcalfield 00228 # 00229 # Set the polmodel from pcalmodel 00230 # 00231 print '--Setting up Polarization models--' 00232 00233 polmodel = {} 00234 for field in pcalmodel.keys() : 00235 spwmodel = {} 00236 # the RLPD is atan2(U,Q) so Q=I*P/I*cos(RLPD) U=I*P/I*sin(RLPD) 00237 for spw in usespwlist: 00238 ipol = pcalmodel[field][spw][0] 00239 fpol = pcalmodel[field][spw][1] 00240 rlpd_deg = pcalmodel[field][spw][2] 00241 rlpd = rlpd_deg*pl.pi/180.0 00242 ppol = ipol*fpol 00243 qpol = ppol*cos(rlpd) 00244 upol = ppol*sin(rlpd) 00245 fluxdensity=[ipol,qpol,upol,0.0] 00246 00247 pmodel = {} 00248 pmodel['rlpd_deg'] = rlpd_deg 00249 pmodel['rlpd'] = rlpd 00250 pmodel['fpol'] = fpol 00251 00252 fmodel = {} 00253 fmodel['flux'] = fluxdensity 00254 fmodel['poln'] = pmodel 00255 spwmodel[spw] = fmodel 00256 00257 polmodel[field] = spwmodel 00258 00259 print "Created polmodel dictionary" 00260 print polmodel 00261 # 00262 #===================================================================== 00263 # Data Import and List 00264 #===================================================================== 00265 # 00266 if ( importmode == 'vla' ): 00267 # 00268 # Import the data from VLA Export to MS 00269 # 00270 print '--ImportVLA--' 00271 default('importvla') 00272 00273 print "Use importvla to read VLA Export and make an MS" 00274 00275 archivefiles = datafile 00276 vis = msfile 00277 bandname = exportband 00278 autocorr = False 00279 antnamescheme = 'new' 00280 project = exportproject 00281 00282 saveinputs('importvla',prefix+'.importvla.saved') 00283 importvla() 00284 elif ( importmode == 'fits' ): 00285 # 00286 # Import the data from VLA Export to MS 00287 # 00288 print '--ImportUVFITS--' 00289 default('importuvfits') 00290 00291 print "Use importuvfits to read UVFITS and make an MS" 00292 00293 fitsfile = datafile 00294 vis = msfile 00295 async = False 00296 00297 saveinputs('importuvfits',prefix+'.importuvfits.saved') 00298 importuvfits() 00299 else: 00300 # 00301 # Copy from msfile 00302 # 00303 print '--MS Copy--' 00304 print "Copying "+datafile+" to "+msfile 00305 00306 os.system('cp -r '+datafile+' '+msfile) 00307 vis = msfile 00308 00309 # 00310 #===================================================================== 00311 # 00312 print '--Listobs--' 00313 00314 print "List summary of MS" 00315 00316 listobs() 00317 00318 ############################################### 00319 ### Begin Task: listobs ### 00320 # 00321 # MeasurementSet Name: 00322 # /home/sandrock/smyers/Testing/2008-03/polcal_20080224/polcal_20080224.cband.raw.ms 00323 # MS Version 2 00324 # 00325 # Observer: unavailable Project: POLCA 00326 # Observation: VLA 00327 # Data records: 318708 Total integration time = 9836.67 seconds 00328 # Observed from 17:10:52 to 19:54:48 00329 # 00330 # ObservationID = 0 ArrayID = 0 00331 # Date Timerange Scan FldId FieldName SpwIds 00332 # 24-Feb-2008/17:10:51.7 - 17:12:08.3 1 0 1924-292 [0, 1] 00333 # 17:21:01.7 - 17:22:18.3 2 1 1743-038 [0, 1] 00334 # 17:34:31.7 - 17:35:48.3 3 2 2202+422 [0, 1] 00335 # 17:45:01.7 - 17:46:18.3 4 3 2253+161 [0, 1] 00336 # 17:55:11.7 - 17:56:28.3 5 4 2136+006 [0, 1] 00337 # 18:08:01.7 - 18:09:18.3 6 5 0137+331 [0, 1] 00338 # 18:22:11.7 - 18:23:58.3 7 6 2355+498 [0, 1] 00339 # 18:32:51.7 - 19:07:58.3 8 2 2202+422 [0, 1] 00340 # 19:20:51.7 - 19:22:18.3 9 5 0137+331 [0, 1] 00341 # 19:32:11.7 - 19:33:48.3 10 7 0319+415 [0, 1] 00342 # 19:42:01.7 - 19:43:18.3 11 8 0359+509 [0, 1] 00343 # 19:53:31.7 - 19:54:48.3 12 2 2202+422 [0, 1] 00344 # Fields: 9 00345 # ID Code Name Right Ascension Declination Epoch 00346 # 0 A 1924-292 19:24:51.06 -29.14.30.12 J2000 00347 # 1 A 1743-038 17:43:58.86 -03.50.04.62 J2000 00348 # 2 A 2202+422 22:02:43.29 +42.16.39.98 J2000 00349 # 3 A 2253+161 22:53:57.75 +16.08.53.56 J2000 00350 # 4 B 2136+006 21:36:38.59 +00.41.54.21 J2000 00351 # 5 B 0137+331 01:37:41.30 +33.09.35.13 J2000 00352 # 6 A 2355+498 23:55:09.46 +49.50.08.34 J2000 00353 # 7 B 0319+415 03:19:48.16 +41.30.42.10 J2000 00354 # 8 B 0359+509 03:59:29.75 +50.57.50.16 J2000 00355 # Spectral Windows: (2 unique spectral windows and 1 unique polarization setups) 00356 # SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 00357 # 0 1 TOPO 4885.1 50000 50000 4885.1 RR RL LR LL 00358 # 1 1 TOPO 4835.1 50000 50000 4835.1 RR RL LR LL 00359 # Feeds: 27: printing first row only 00360 # Antenna Spectral Window # Receptors Polarizations 00361 # 1 -1 2 [ R, L] 00362 # Antennas: 27: 00363 # ID Name Station Diam. Long. Lat. 00364 # 0 EA24 VLA:W12 25.0 m -107.37.37.4 +33.53.44.2 00365 # 1 EA16 VLA:W6 25.0 m -107.37.15.6 +33.53.56.4 00366 # 2 EA01 VLA:W10 25.0 m -107.37.28.9 +33.53.48.9 00367 # 3 EA19 VLA:W4 25.0 m -107.37.10.8 +33.53.59.1 00368 # 4 VA08 VLA:W16 25.0 m -107.37.57.4 +33.53.33.0 00369 # 5 EA17 VLA:W14 25.0 m -107.37.46.9 +33.53.38.9 00370 # 6 VA06 VLA:W8 25.0 m -107.37.21.6 +33.53.53.0 00371 # 7 VA22 VLA:W2 25.0 m -107.37.07.4 +33.54.00.9 00372 # 8 EA04 UNKNOWN 25.0 m -107.37.41.3 +33.53.42.0 00373 # 9 VA20 VLA:E12 25.0 m -107.36.31.7 +33.53.48.5 00374 # 10 VA15 VLA:E4 25.0 m -107.37.00.8 +33.53.59.7 00375 # 11 VA28 VLA:E6 25.0 m -107.36.55.6 +33.53.57.7 00376 # 12 VA10 VLA:E8 25.0 m -107.36.48.9 +33.53.55.1 00377 # 13 EA14 VLA:E16 25.0 m -107.36.09.8 +33.53.40.0 00378 # 14 EA11 VLA:E10 25.0 m -107.36.40.9 +33.53.52.0 00379 # 15 VA03 VLA:E14 25.0 m -107.36.21.3 +33.53.44.5 00380 # 16 EA23 VLA:E18 25.0 m -107.35.57.2 +33.53.35.1 00381 # 17 EA21 VLA:E2 25.0 m -107.37.04.4 +33.54.01.1 00382 # 18 VA12 VLA:N4 25.0 m -107.37.06.5 +33.54.06.1 00383 # 19 VA02 VLA:N20 25.0 m -107.37.13.2 +33.55.09.5 00384 # 20 EA13 VLA:N16 25.0 m -107.37.10.9 +33.54.48.0 00385 # 21 EA26 VLA:N32 25.0 m -107.37.22.0 +33.56.33.6 00386 # 22 EA25 VLA:N24 25.0 m -107.37.16.1 +33.55.37.7 00387 # 23 VA09 VLA:N8 25.0 m -107.37.07.5 +33.54.15.8 00388 # 24 EA18 VLA:N12 25.0 m -107.37.09.0 +33.54.30.0 00389 # 25 VA07 VLA:N36 25.0 m -107.37.25.6 +33.57.07.6 00390 # 26 VA27 VLA:N28 25.0 m -107.37.18.7 +33.56.02.5 00391 # 00392 # 00393 # Tables: 00394 # MAIN 318708 rows 00395 # ANTENNA 27 rows 00396 # DATA_DESCRIPTION 2 rows 00397 # DOPPLER 2 rows 00398 # FEED 27 rows 00399 # FIELD 9 rows 00400 # FLAG_CMD <empty> 00401 # FREQ_OFFSET <absent> 00402 # HISTORY 6 rows 00403 # OBSERVATION 1 row 00404 # POINTING <empty> 00405 # POLARIZATION 1 row 00406 # PROCESSOR <empty> 00407 # SOURCE 9 rows 00408 # SPECTRAL_WINDOW 2 rows 00409 # STATE <empty> 00410 # SYSCAL <absent> 00411 # WEATHER <absent> 00412 # 00413 ### End Task: listobs ### 00414 ############################################### 00415 00416 # Note that the antennas are out of order as loaded by importvla 00417 00418 # 00419 #===================================================================== 00420 # Data Flagging if needed 00421 #===================================================================== 00422 # 00423 if ( myquackinterval > 0.0 ): 00424 # 00425 # First quack the data 00426 # 00427 print '--Flagdata--' 00428 default('tflagdata') 00429 00430 print "Quacking scan beginnings using interval "+str(myquackinterval) 00431 00432 vis = msfile 00433 correlation = '' 00434 field = '' 00435 antenna = '' 00436 spw = usespw 00437 mode = 'quack' 00438 quackinterval = myquackinterval 00439 00440 saveinputs('flagdata',prefix+'.flagdata.quack.saved') 00441 tflagdata() 00442 00443 # 00444 # Use Flagmanager to save a copy of the flags so far 00445 # 00446 default('flagmanager') 00447 00448 print "Now will use flagmanager to save the flags" 00449 00450 vis = msfile 00451 mode = 'save' 00452 versionname = 'quack' 00453 comment = 'Quack '+str(myquackinterval) 00454 merge = 'replace' 00455 00456 saveinputs('flagmanager',prefix+'.flagmanager.quack.saved') 00457 flagmanager() 00458 00459 00460 # 00461 if (flagants != '' and not flagants.isspace() ): 00462 print '--Flagdata--' 00463 default('tflagdata') 00464 00465 print "Flag all data to AN "+flagants 00466 00467 vis = msfile 00468 correlation = '' 00469 field = '' 00470 spw = usespw 00471 mode = 'manual' 00472 antenna = flagants 00473 00474 saveinputs('flagdata',prefix+'.flagdata.ants.saved') 00475 tflagdata() 00476 00477 # 00478 # Use Flagmanager to save a copy of the flags so far 00479 # 00480 default('flagmanager') 00481 00482 print "Now will use flagmanager to save the flags" 00483 00484 vis = msfile 00485 mode = 'save' 00486 versionname = 'antflags' 00487 comment = 'flag AN '+flagants 00488 merge = 'replace' 00489 00490 saveinputs('flagmanager',prefix+'.flagmanager.ants.saved') 00491 flagmanager() 00492 00493 00494 # 00495 #===================================================================== 00496 # Calibration 00497 #===================================================================== 00498 # 00499 # Set the fluxes of the primary calibrator(s) 00500 # 00501 if ( setjymode == 'flux' ): 00502 print '--Setjy--' 00503 default('setjy') 00504 00505 vis = msfile 00506 00507 print "Use setjy to set flux of "+fluxcalfield+" to point model" 00508 field = fluxcalfield 00509 spw = usespw 00510 00511 # If we need a model for flux calibrator then put this here 00512 modimage = fluxcaldir + fluxcalmodel 00513 00514 # Loop over spw 00515 for spw in usespwlist: 00516 fluxdensity = fcalmodel[fluxcalfield][spw] 00517 print "Setting SPW "+spw+" to "+str(fluxdensity) 00518 saveinputs('setjy',prefix+'.setjy.'+spw+'.saved') 00519 setjy() 00520 00521 elif ( setjymode == 'ft' ): 00522 print '--FT--' 00523 00524 default('ft') 00525 vis = msfile 00526 field = fluxcalfield 00527 00528 for spw in usespwlist: 00529 model = fluxcaldir + fluxcalmodel+'_'+spw+'_IQUV.model' 00530 print "Use FT to set model"+model 00531 saveinputs('ft',prefix+'.ft.0.saved') 00532 ft() 00533 00534 else: 00535 print '--Setjy--' 00536 default('setjy') 00537 00538 vis = msfile 00539 00540 print "Use setjy to set flux of "+fluxcalfield 00541 field = fluxcalfield 00542 spw = usespw 00543 00544 # If we need a model or fluxdensities then put those here 00545 modimage = fluxcaldir + fluxcalmodel 00546 00547 saveinputs('setjy',prefix+'.setjy.saved') 00548 setjy() 00549 # 00550 # You should see something like this in the logger and casapy.log file: 00551 # 00552 # 0137+331 spwid= 0 [I=5.405, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00553 # 0137+331 spwid= 1 [I=5.458, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00554 00555 # cf. AIPS 00556 # SETJY '0137+331 ' IF = 1 FLUX = 5.4054 (Jy calcd) 00557 # SETJY '0137+331 ' IF = 2 FLUX = 5.4585 (Jy calcd) 00558 00559 print "Look in logger for the fluxes (should be 5.405 and 5.458 Jy)" 00560 00561 #===================================================================== 00562 # 00563 # Initial gain calibration 00564 # 00565 print '--Gaincal--' 00566 default('gaincal') 00567 00568 print "Solve for antenna gains on sources "+gaincalfield 00569 print "We have 2 single-channel continuum spw" 00570 00571 vis = msfile 00572 00573 # set the name for the output gain caltable 00574 print "Output gain table name is "+gtable 00575 caltable = gtable 00576 00577 # All fields are calibrators 00578 # We have 2 IFs (SPW 0,1) with one channel each 00579 00580 # Assemble field string from gaincalfield list 00581 field = fieldgain 00582 print "Calibrating using fields "+field 00583 00584 # Calibrate these spw 00585 spw = usespw 00586 00587 # a-priori calibration application 00588 gaincurve = usegaincurve 00589 opacity = gainopacity 00590 00591 # do not apply parallactic angle correction 00592 parang = False 00593 00594 # G solutions for both amplitude and phase using gainsolint 00595 gaintype = 'G' 00596 solint = gainsolint 00597 calmode = 'ap' 00598 00599 # reference antenna 00600 refant = calrefant 00601 00602 # minimum SNR 3 00603 minsnr = 3 00604 00605 saveinputs('gaincal',prefix+'.gaincal.saved') 00606 gaincal() 00607 00608 # use plotcal to view or listcal to list 00609 00610 #===================================================================== 00611 # 00612 # List gain calibration 00613 # 00614 print '--Listcal--' 00615 00616 listfile = caltable + '.list' 00617 00618 print "Listing calibration to file "+listfile 00619 00620 listcal() 00621 00622 # 00623 #===================================================================== 00624 # 00625 # Bootstrap flux scale 00626 # 00627 print '--Fluxscale--' 00628 default('fluxscale') 00629 00630 print "Use fluxscale to rescale gain table to make new one" 00631 00632 vis = msfile 00633 00634 # set the name for the output rescaled caltable 00635 ftable = prefix + '.fluxscale' 00636 fluxtable = ftable 00637 00638 print "Output scaled gain cal table is "+ftable 00639 00640 # point to our first gain cal table 00641 caltable = gtable 00642 00643 # use the source we did setjy on as our flux standard reference 00644 reference = fluxcalfield 00645 00646 # transfer the flux to all our other sources 00647 # to bring amplitues in line with the absolute scale 00648 transfer = fieldgain 00649 00650 saveinputs('fluxscale',prefix+'.fluxscale.saved') 00651 fluxscale() 00652 00653 # You should see in the logger something like: 00654 # Found reference field(s): 0137+331 00655 # Found transfer field(s): 1924-292 1743-038 2202+422 2253+161 2136+006 2355+498 0319+415 0359+509 00656 # Flux density for 1924-292 in SpW=0 is: 8.25145 +/- 0.00988121 (SNR = 835.065, nAnt= 13) 00657 # Flux density for 1924-292 in SpW=1 is: 8.22457 +/- 0.0140951 (SNR = 583.505, nAnt= 13) 00658 # Flux density for 1743-038 in SpW=0 is: 5.31336 +/- 0.00603626 (SNR = 880.239, nAnt= 13) 00659 # Flux density for 1743-038 in SpW=1 is: 5.3184 +/- 0.00480634 (SNR = 1106.54, nAnt= 13) 00660 # Flux density for 2202+422 in SpW=0 is: 2.46545 +/- 0.00335055 (SNR = 735.833, nAnt= 13) 00661 # Flux density for 2202+422 in SpW=1 is: 2.46072 +/- 0.00353799 (SNR = 695.512, nAnt= 13) 00662 # Flux density for 2253+161 in SpW=0 is: 8.74607 +/- 0.0142334 (SNR = 614.474, nAnt= 13) 00663 # Flux density for 2253+161 in SpW=1 is: 8.77219 +/- 0.0102289 (SNR = 857.587, nAnt= 13) 00664 # Flux density for 2136+006 in SpW=0 is: 9.97863 +/- 0.013815 (SNR = 722.303, nAnt= 13) 00665 # Flux density for 2136+006 in SpW=1 is: 9.99001 +/- 0.0170089 (SNR = 587.339, nAnt= 13) 00666 # Flux density for 2355+498 in SpW=0 is: 1.29395 +/- 0.00181169 (SNR = 714.221, nAnt= 13) 00667 # Flux density for 2355+498 in SpW=1 is: 1.29893 +/- 0.00217214 (SNR = 597.995, nAnt= 13) 00668 # Flux density for 0319+415 in SpW=0 is: 13.5742 +/- 0.0221722 (SNR = 612.218, nAnt= 13) 00669 # Flux density for 0319+415 in SpW=1 is: 13.5481 +/- 0.0230828 (SNR = 586.932, nAnt= 13) 00670 # Flux density for 0359+509 in SpW=0 is: 5.13982 +/- 0.00906505 (SNR = 566.993, nAnt= 13) 00671 # Flux density for 0359+509 in SpW=1 is: 5.10322 +/- 0.00990264 (SNR = 515.339, nAnt= 13) 00672 # Storing result in polcal_20080224.cband.vla_3c84.fluxscale 00673 # Writing solutions to table: polcal_20080224.cband.vla_3c84.fluxscale 00674 00675 #===================================================================== 00676 # 00677 # List fluxscale table 00678 # 00679 print '--Listcal--' 00680 00681 caltable = ftable 00682 listfile = caltable + '.list' 00683 00684 print "Listing calibration to file "+listfile 00685 00686 listcal() 00687 00688 #===================================================================== 00689 # 00690 # Plot final gain calibration 00691 # 00692 print '--Plotcal--' 00693 00694 iteration = '' 00695 showgui = False 00696 00697 xaxis = 'time' 00698 yaxis = 'amp' 00699 figfile = caltable + '.plot.amp.png' 00700 print "Plotting calibration to file "+figfile 00701 saveinputs('plotcal',prefix+'.plotcal.fluxscale.amp.saved') 00702 plotcal() 00703 00704 xaxis = 'time' 00705 yaxis = 'phase' 00706 figfile = caltable + '.plot.phase.png' 00707 print "Plotting calibration to file "+figfile 00708 saveinputs('plotcal',prefix+'.plotcal.fluxscale.phase.saved') 00709 plotcal() 00710 00711 xaxis = 'antenna' 00712 yaxis = 'amp' 00713 figfile = caltable + '.plot.antamp.png' 00714 print "Plotting calibration to file "+figfile 00715 saveinputs('plotcal',prefix+'.plotcal.fluxscale.antamp.saved') 00716 plotcal() 00717 00718 if ( setpolmodel and polcalmode.count('X') > 0 ): 00719 # 00720 # ===================================================================== 00721 # 00722 # Now run setjy to (re)set model for polxfield 00723 # 00724 print '--Setjy--' 00725 default('setjy') 00726 00727 vis = msfile 00728 00729 print "Use setjy to set IQU fluxes of "+polxfield 00730 field = polxfield 00731 00732 for spw in usespwlist: 00733 fluxdensity = polmodel[field][spw]['flux'] 00734 00735 saveinputs('setjy',prefix+'.setjy.polspw.'+spw+'.saved') 00736 setjy() 00737 00738 00739 #===================================================================== 00740 # 00741 # Polarization (D-term) calibration 00742 # 00743 print '--PolCal--' 00744 default('polcal') 00745 00746 print "Polarization D-term Calibration (linear approx) on "+polcalfield 00747 00748 vis = msfile 00749 00750 # Start with the un-fluxscaled gain table 00751 gaintable = gtable 00752 00753 # use settings from gaincal 00754 gaincurve = usegaincurve 00755 opacity = gainopacity 00756 00757 # Output table 00758 ptable = prefix + '.pcal' 00759 caltable = ptable 00760 00761 # Use an unpolarized source or a source tracked through a range of PA 00762 field = polcalfield 00763 spw = usespw 00764 00765 selectdata=True 00766 uvrange = polduvrange 00767 00768 # Polcal mode 00769 poltype = polcalmode 00770 00771 # Currently 1-day timescale is hardwired 00772 solint = 86400. 00773 00774 # reference antenna 00775 refant = calrefant 00776 00777 # minimum SNR 3 00778 minsnr = 3 00779 00780 saveinputs('polcal',prefix+'.polcal.saved') 00781 polcal() 00782 00783 # You should see something like: 00784 # Fractional polarization solution for 2202+422 (spw = 0): 00785 # : Q = 0.00356182, U = 0.0717148 (P = 0.0718032, X = 43.5783 deg) 00786 # Fractional polarization solution for 2202+422 (spw = 1): 00787 # : Q = -0.00561314, U = -0.0720833 (P = 0.0723015, X = -47.2263 deg) 00788 00789 #===================================================================== 00790 # 00791 # List polcal solutions 00792 # 00793 print '--Listcal--' 00794 00795 listfile = caltable + '.list' 00796 00797 print "Listing calibration to file "+listfile 00798 00799 listcal() 00800 00801 #===================================================================== 00802 # 00803 # Plot polcal solutions 00804 # 00805 print '--Plotcal--' 00806 00807 iteration = '' 00808 showgui = False 00809 00810 xaxis = 'real' 00811 yaxis = 'imag' 00812 figfile = caltable + '.plot.reim.png' 00813 print "Plotting calibration to file "+figfile 00814 saveinputs('plotcal',prefix+'.plotcal.polcal.d.reim.saved') 00815 plotcal() 00816 00817 xaxis = 'antenna' 00818 yaxis = 'amp' 00819 figfile = caltable + '.plot.antamp.png' 00820 print "Plotting calibration to file "+figfile 00821 saveinputs('plotcal',prefix+'.plotcal.polcal.d.antamp.saved') 00822 plotcal() 00823 00824 xaxis = 'antenna' 00825 yaxis = 'phase' 00826 figfile = caltable + '.plot.antphase.png' 00827 print "Plotting calibration to file "+figfile 00828 saveinputs('plotcal',prefix+'.plotcal.polcal.d.antphase.saved') 00829 plotcal() 00830 00831 xaxis = 'antenna' 00832 yaxis = 'snr' 00833 figfile = caltable + '.plot.antsnr.png' 00834 print "Plotting calibration to file "+figfile 00835 saveinputs('plotcal',prefix+'.plotcal.polcal.d.antsnr.saved') 00836 plotcal() 00837 00838 #===================================================================== 00839 # Do Chi (X) pol angle calibration if possible 00840 #===================================================================== 00841 # 00842 dopolx = False 00843 if ( pcalmodel.has_key(polxfield) ): 00844 dopolx = True 00845 00846 if ( setpolmodel and not polcalmode.count('X') > 0 ): 00847 # 00848 # ============================================================= 00849 # 00850 # Now run setjy if we havent already 00851 # 00852 00853 print '--Setjy--' 00854 default('setjy') 00855 00856 vis = msfile 00857 00858 print "Use setjy to set IQU fluxes of "+polxfield 00859 field = polxfield 00860 00861 for spw in usespwlist: 00862 fluxdensity = polmodel[field][spw]['flux'] 00863 00864 saveinputs('setjy',prefix+'.setjy.polspw.'+spw+'.saved') 00865 setjy() 00866 00867 00868 # 00869 # ===================================================================== 00870 # 00871 # Polarization (X-term) calibration 00872 # 00873 print '--PolCal--' 00874 default('polcal') 00875 00876 print "Polarization R-L Phase Calibration (linear approx)" 00877 00878 vis = msfile 00879 00880 # Start with the G and D tables 00881 gaintable = [gtable,ptable] 00882 00883 # use settings from gaincal 00884 gaincurve = usegaincurve 00885 opacity = gainopacity 00886 00887 # Output table 00888 xtable = prefix + '.polx' 00889 caltable = xtable 00890 00891 # previously set with setjy 00892 field = polxfield 00893 spw = usespw 00894 00895 selectdata=True 00896 uvrange = polxuvrange 00897 00898 # Solve for Chi 00899 poltype = 'X' 00900 solint = 86400. 00901 00902 # reference antenna 00903 refant = calrefant 00904 00905 # minimum SNR 3 00906 minsnr = 3 00907 00908 saveinputs('polcal',prefix+'.polcal.X.saved') 00909 polcal() 00910 00911 # You should get something like: 00912 # Position angle offset solution for 0137+331 (spw = 0) = 72.437 deg. 00913 # Position angle offset solution for 0137+331 (spw = 1) = -21.0703 deg. 00914 00915 # 00916 # ===================================================================== 00917 # 00918 # List polcal solutions 00919 # 00920 #print '--Listcal--' 00921 00922 #listfile = caltable + '.list' 00923 00924 #print "Listing calibration to file "+listfile 00925 00926 #listcal() 00927 00928 # 00929 # ===================================================================== 00930 # 00931 # Plot polcal solutions 00932 # 00933 print '--Plotcal--' 00934 00935 xaxis = 'antenna' 00936 yaxis = 'phase' 00937 iteration = '' 00938 00939 showgui = False 00940 figfile = caltable + '.plot.png' 00941 00942 print "Plotting calibration to file "+figfile 00943 saveinputs('plotcal',prefix+'.plotcal.polcal.x.antphase.saved') 00944 plotcal() 00945 00946 else: 00947 if (polxfield != '' and not polxfield.isspace() ): 00948 print "DO NOT HAVE PCALMODEL FOR "+polxfield 00949 print "PCALMODEL = ",pcalmodel 00950 00951 00952 #===================================================================== 00953 # 00954 # Correct the data 00955 # (This will put calibrated data into the CORRECTED_DATA column) 00956 # 00957 # First using gaincalfield 00958 # 00959 print '--ApplyCal--' 00960 default('applycal') 00961 00962 print "This will apply the calibration to the DATA" 00963 print "Fills CORRECTED_DATA" 00964 00965 vis = msfile 00966 00967 # Start with the fluxscaled G table, the D table, and the X table 00968 if (dopolx): 00969 gaintable = [ftable,ptable,xtable] 00970 else: 00971 gaintable = [ftable,ptable] 00972 00973 # use settings from gaincal 00974 gaincurve = usegaincurve 00975 opacity = gainopacity 00976 00977 # select all the data 00978 spw = usespw 00979 selectdata = False 00980 00981 # IMPORTANT set parang=True for polarization 00982 parang = True 00983 00984 # use the list of gain calibrators, apply to themselves 00985 field = fieldgain 00986 gainselect = field 00987 print "Applying calibration to gain calibrators "+field 00988 00989 saveinputs('applycal',prefix+'.applycal.saved') 00990 applycal() 00991 00992 if ( len(targets) > 0 ): 00993 # 00994 # Now with targets if any (transfer from gaincalfield) 00995 # 00996 # Assemble field string from target list 00997 field = fieldtargets 00998 print "Applying calibration to targets "+field 00999 01000 saveinputs('applycal',prefix+'.applycal.targets.saved') 01001 applycal() 01002 01003 # 01004 #===================================================================== 01005 # 01006 # Now write out the corrected data 01007 # 01008 print '--Split--' 01009 default('split') 01010 01011 vis = msfile 01012 01013 # Now we write out the corrected data to a new MS 01014 01015 # Make an output vis file 01016 srcsplitms = prefix + '.split.ms' 01017 outputvis = srcsplitms 01018 01019 # Select all data 01020 field = '' 01021 01022 # Have to split all spw to preserve numbering 01023 spw = '' 01024 01025 # pick off the CORRECTED_DATA column 01026 datacolumn = 'corrected' 01027 01028 print "Split CORRECTED_DATA into DATA in new ms "+srcsplitms 01029 01030 saveinputs('split',prefix+'.split.saved') 01031 split() 01032 01033 # 01034 #===================================================================== 01035 # 01036 # Plot up the visibilities for the main calibrators 01037 # 01038 print '--Plotxy--' 01039 default('plotxy') 01040 01041 vis = srcsplitms 01042 01043 field = fluxcalfield 01044 spw = '' 01045 01046 selectdata=True 01047 01048 xaxis = 'uvdist' 01049 01050 interactive=False 01051 01052 correlation='RR LL' 01053 yaxis = 'amp' 01054 figfile = prefix+'.split.'+field+'.uvplot.amp.png' 01055 saveinputs('plotxy',prefix+'.plotxy.'+field+'.amp.saved') 01056 plotxy() 01057 01058 correlation='RL LR' 01059 yaxis = 'phase' 01060 figfile = prefix+'.split.'+field+'.uvplot.rlphase.png' 01061 saveinputs('plotxy',prefix+'.plotxy.'+field+'.rlphase.saved') 01062 plotxy() 01063 01064 if ( polcalfield != fluxcalfield ): 01065 # Now the poln calibrator 01066 field = polcalfield 01067 01068 correlation='RR LL' 01069 yaxis = 'amp' 01070 figfile = prefix+'.split.'+field+'.uvplot.amp.png' 01071 saveinputs('plotxy',prefix+'.plotxy.'+field+'.amp.saved') 01072 plotxy() 01073 01074 correlation='RL LR' 01075 yaxis = 'phase' 01076 figfile = prefix+'.split.'+field+'.uvplot.rlphase.png' 01077 saveinputs('plotxy',prefix+'.plotxy.'+field+'.rlphase.saved') 01078 plotxy() 01079 01080 # 01081 #===================================================================== 01082 # CLEAN the sources 01083 #===================================================================== 01084 01085 clnmodel = {} 01086 # 01087 #===================================================================== 01088 # Loop over sources and spw 01089 # Set up for new clean in patch 2 01090 # 01091 for src in srclist: 01092 01093 srcmodel = {} 01094 01095 for spwid in usespwlist: 01096 01097 print '-- Clean '+src+' spw '+spwid+' --' 01098 default('clean') 01099 01100 field = src 01101 spw = spwid 01102 01103 # Pick up our split source data 01104 vis = srcsplitms 01105 01106 # Make an image root file name 01107 imname1 = prefix + '.' + src + '.' + spwid + '.clean' 01108 imagename = imname1 01109 01110 print " Output images will be prefixed with "+imname1 01111 01112 # Set up the output continuum image (single plane mfs) 01113 mode = 'mfs' 01114 01115 # All polarizations 01116 stokes = 'IQUV' 01117 01118 # Use chose clean style 01119 psfmode = clnalg 01120 csclean = usecsclean 01121 01122 imsize = [clnimsize,clnimsize] 01123 cell = [clncell,clncell] 01124 01125 # Standard gain factor 0.1 01126 gain = 0.1 01127 01128 niter = clniter 01129 01130 threshold = clthreshold 01131 01132 # Set up the weighting 01133 # Use Briggs weighting (a moderate value, on the uniform side) 01134 weighting = 'briggs' 01135 robust = 0.5 01136 # Use natural weighting 01137 weighting = 'natural' 01138 01139 # Use the cleanbox 01140 mask = myclnbox 01141 01142 saveinputs('clean',prefix+'.clean.'+src+'.'+spwid+'.saved') 01143 clean() 01144 01145 # Set up variables 01146 clnimage1 = imname1+'.image' 01147 clnmodel1 = imname1+'.model' 01148 clnresid1 = imname1+'.residual' 01149 clnmask1 = imname1+'.mask' 01150 clnpsf1 = imname1+'.psf' 01151 clnflux1 = imname1+'.flux' 01152 01153 # 01154 # ===================================================================== 01155 # 01156 # Get some statistics of the clean image 01157 # 01158 default('imstat') 01159 01160 field = src 01161 spw = spwid 01162 01163 # Use the clean box 01164 mybox = str(clnblc)+','+str(clnblc)+','+str(clntrc)+','+str(clntrc) 01165 01166 spwmodel = {} 01167 01168 spwstats = {} 01169 spwfluxes = {} 01170 spwsum = {} 01171 spwmod = {} 01172 01173 for stokes in ['I','Q','U','V']: 01174 01175 # Use the clean image 01176 imagename = clnimage1 01177 box = mybox 01178 01179 saveinputs('imstat',prefix+'.imstat.'+src+'.'+spwid+'.'+stokes+'.saved') 01180 xstat = imstat() 01181 01182 spwstats[stokes] = xstat 01183 01184 # Peak (max or min) in box 01185 xmax = xstat['max'][0] 01186 xmin = xstat['min'][0] 01187 if( abs(xmin) > abs(xmax) ): 01188 xpol = xmin 01189 else: 01190 xpol = xmax 01191 01192 spwfluxes[stokes]= xpol 01193 01194 # Integrated flux in box 01195 xsum = xstat['flux'][0] 01196 spwsum[stokes]= xsum 01197 01198 # Use the clean model and no box 01199 imagename = clnmodel1 01200 box = '' 01201 01202 saveinputs('imstat',prefix+'.imstat.'+src+'.'+spwid+'.'+stokes+'.model.saved') 01203 xstat = imstat() 01204 # Integrated flux in image 01205 xmod = xstat['sum'][0] 01206 spwmod[stokes]= xmod 01207 01208 # Done with stokes 01209 01210 spwmodel['stat'] = spwstats 01211 spwmodel['flux'] = spwfluxes 01212 spwmodel['integ'] = spwsum 01213 spwmodel['model'] = spwmod 01214 01215 # Use ia tool for pixel values in the restored image 01216 imagename = clnimage1 01217 # Get image values at the reference pixel 01218 spwref = {} 01219 ia.open(imagename) 01220 # 01221 # Stokes I 01222 ipix = ia.pixelvalue() 01223 # Get reference pixel 01224 xref = ipix['pixel'][0] 01225 yref = ipix['pixel'][1] 01226 iflx = ipix['value']['value'] 01227 spwref['I'] = iflx 01228 # 01229 # Stokes Q 01230 qpix = ia.pixelvalue([xref,yref,1,0]) 01231 qflx = qpix['value']['value'] 01232 spwref['Q'] = qflx 01233 # 01234 # Stokes U 01235 upix = ia.pixelvalue([xref,yref,2,0]) 01236 uflx = upix['value']['value'] 01237 spwref['U'] = uflx 01238 # 01239 # Stokes V 01240 vpix = ia.pixelvalue([xref,yref,3,0]) 01241 vflx = vpix['value']['value'] 01242 spwref['V'] = vflx 01243 # 01244 # Polarization quantities 01245 pflx = sqrt( qflx**2 + uflx**2 ) 01246 fflx = pflx/iflx 01247 xflx = atan2(uflx,qflx)*180.0/pi 01248 spwref['P'] = pflx 01249 spwref['F'] = fflx 01250 spwref['X'] = xflx 01251 spwref['xref'] = xref 01252 spwref['yref'] = yref 01253 # 01254 01255 # Now the values at the maximum of I 01256 spwmax = {} 01257 # 01258 # Pull the maxpos of I 01259 xref = spwstats['I']['maxpos'][0] 01260 yref = spwstats['I']['maxpos'][1] 01261 # 01262 # Stokes I 01263 iflx = spwstats['I']['max'][0] 01264 spwmax['I'] = iflx 01265 # 01266 # Stokes Q 01267 qpix = ia.pixelvalue([xref,yref,1,0]) 01268 qflx = qpix['value']['value'] 01269 spwmax['Q'] = qflx 01270 # 01271 # Stokes U 01272 upix = ia.pixelvalue([xref,yref,2,0]) 01273 uflx = upix['value']['value'] 01274 spwmax['U'] = uflx 01275 # 01276 # Stokes V 01277 vpix = ia.pixelvalue([xref,yref,3,0]) 01278 vflx = vpix['value']['value'] 01279 spwmax['V'] = vflx 01280 01281 spwmax['xref'] = xref 01282 spwmax['yref'] = yref 01283 # Done with ia tool 01284 ia.close() 01285 01286 spwmodel['refval'] = spwref 01287 spwmodel['maxval'] = spwmax 01288 01289 srcmodel[spwid] = spwmodel 01290 01291 # Done with spw 01292 01293 clnmodel[src] = srcmodel 01294 01295 # Done with srcs 01296 # 01297 #===================================================================== 01298 # Report Final Stats 01299 #===================================================================== 01300 # 01301 print 'Results for '+prefix+' :' 01302 print "" 01303 01304 import datetime 01305 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 01306 01307 outfile = 'out.'+prefix+'.'+datestring+'.log' 01308 logfile=open(outfile,'w') 01309 print >>logfile,'Results for '+prefix+' :' 01310 print >>logfile,"" 01311 01312 if ( polmodel.has_key(polxfield) ): 01313 # Check RL phase offset on X calibrator 01314 print "R-L phase residual from image of "+polxfield 01315 print "" 01316 print >>logfile,"R-L phase residual from image of "+polxfield+" :" 01317 print >>logfile,"" 01318 01319 src = polxfield 01320 rlcor = {} 01321 01322 for spwid in usespwlist: 01323 ipol = clnmodel[src][spwid]['flux']['I'] 01324 qpol = clnmodel[src][spwid]['flux']['Q'] 01325 upol = clnmodel[src][spwid]['flux']['U'] 01326 vpol = clnmodel[src][spwid]['flux']['V'] 01327 rlpd = atan2(upol,qpol) 01328 rlpdcal = polmodel[src][spwid]['poln']['rlpd'] 01329 rlpcor = rlpdcal - rlpd 01330 scor = sin(rlpcor); ccor = cos(rlpcor); rlpcor = atan2(scor,ccor) 01331 rlcor[spwid] = rlpcor 01332 rlpcor_deg = rlpcor*180.0/pl.pi 01333 01334 print "R-L Phase Correction SPW "+spwid+" = %7.2f deg" % rlpcor_deg 01335 print >>logfile,"R-L Phase Correction SPW "+spwid+" = %7.2f deg" % rlpcor_deg 01336 01337 01338 # 01339 #===================================================================== 01340 # 01341 # Loop over sources and spw 01342 # 01343 print "" 01344 print "Final Stats:" 01345 print "" 01346 01347 print >>logfile,"" 01348 print >>logfile,"Final Stats:" 01349 print >>logfile,"" 01350 01351 for src in srclist: 01352 01353 print "Source "+src+" :" 01354 print >>logfile,"Source "+src+" :" 01355 01356 for spwid in usespwlist: 01357 01358 field = src 01359 spw = spwid 01360 01361 # Get fluxes from images 01362 01363 ipol = clnmodel[src][spwid]['flux']['I'] 01364 qpol = clnmodel[src][spwid]['flux']['Q'] 01365 upol = clnmodel[src][spwid]['flux']['U'] 01366 vpol = clnmodel[src][spwid]['flux']['V'] 01367 01368 # Now get polarization results 01369 01370 ppol = sqrt(qpol**2 + upol**2) 01371 fpol = ppol/ipol 01372 rlpd = atan2(upol,qpol) 01373 rlpd_deg = rlpd*180.0/pl.pi 01374 01375 #print ' spw %s CASA I = %7.3f Q = %7.3f U = %7.3f V = %7.3f ' % (spwid,ipol,qpol,upol,vpol) 01376 print ' spw %s CASA I = %7.3f P = %7.3f F = %7.4f X = %7.2f deg' % (spwid,ipol,ppol,fpol,rlpd_deg) 01377 print >>logfile,' spw %s CASA I = %7.3f P = %7.3f F = %7.4f X = %7.2f deg' % (spwid,ipol,ppol,fpol,rlpd_deg) 01378 01379 if (aipsmodel.has_key(src)): 01380 iflx = aipsmodel[src][spwid][0]/1000.0 01381 fflx = aipsmodel[src][spwid][1] 01382 rlaips_deg = aipsmodel[src][spwid][2] 01383 rlaips = rlaips_deg*pl.pi/180.0 01384 01385 pflx = iflx*fflx 01386 qflx = pflx*cos(rlaips) 01387 uflx = pflx*sin(rlaips) 01388 vflx = 0.0 01389 01390 print ' spw %s AIPS I = %7.3f P = %7.3f F = %7.4f X = %7.2f deg' % (spwid,iflx,pflx,fflx,rlaips_deg) 01391 print >>logfile,' spw %s AIPS I = %7.3f P = %7.3f F = %7.4f X = %7.2f deg' % (spwid,iflx,pflx,fflx,rlaips_deg) 01392 01393 01394 # Done with spw 01395 print "" 01396 print >>logfile,"" 01397 01398 # Should see something like: 01399 # 01400 # R-L phase residual from image of 0137+331 : 01401 # 01402 # R-L Phase Correction SPW 0 = 0.50 deg 01403 # R-L Phase Correction SPW 1 = 0.17 deg 01404 # 01405 # Final Stats: 01406 # 01407 # Source 0137+331 : 01408 # spw 0 CASA I = 5.320 P = 0.229 F = 0.0430 X = -148.50 deg 01409 # spw 1 CASA I = 5.362 P = 0.236 F = 0.0439 X = -148.17 deg 01410 # 01411 # Source 2202+422 : 01412 # spw 0 CASA I = 2.570 P = 0.185 F = 0.0719 X = -57.71 deg 01413 # spw 1 CASA I = 2.571 P = 0.188 F = 0.0733 X = -53.43 deg 01414 # 01415 # Source 1743-038 : 01416 # spw 0 CASA I = 5.545 P = 0.073 F = 0.0131 X = -1.69 deg 01417 # spw 1 CASA I = 5.543 P = 0.078 F = 0.0141 X = 4.38 deg 01418 # 01419 # Source 1924-292 : 01420 # spw 0 CASA I = 8.243 P = 0.058 F = 0.0070 X = 12.19 deg 01421 # spw 1 CASA I = 8.184 P = 0.073 F = 0.0089 X = 7.27 deg 01422 # 01423 # Source 2136+006 : 01424 # spw 0 CASA I = 10.491 P = 0.153 F = 0.0146 X = -171.61 deg 01425 # spw 1 CASA I = 10.521 P = 0.145 F = 0.0138 X = -160.12 deg 01426 # 01427 # Source 2253+161 : 01428 # spw 0 CASA I = 9.113 P = 0.527 F = 0.0578 X = 4.42 deg 01429 # spw 1 CASA I = 9.096 P = 0.524 F = 0.0576 X = 4.72 deg 01430 # 01431 # Source 2355+498 : 01432 # spw 0 CASA I = 1.340 P = 0.003 F = 0.0022 X = -162.58 deg 01433 # spw 1 CASA I = 1.355 P = 0.003 F = 0.0023 X = 151.08 deg 01434 # 01435 # Source 0319+415 : 01436 # spw 0 CASA I = 14.147 P = 0.041 F = 0.0029 X = -79.67 deg 01437 # spw 1 CASA I = 14.229 P = 0.043 F = 0.0030 X = -130.02 deg 01438 # 01439 # Source 0359+509 : 01440 # spw 0 CASA I = 5.351 P = 0.090 F = 0.0168 X = -124.92 deg 01441 # spw 1 CASA I = 5.333 P = 0.094 F = 0.0177 X = -128.78 deg 01442 # 01443 #===================================================================== 01444 # 01445 # Done 01446 # 01447 logfile.close() 01448 print "Results are in "+outfile 01449 01450 # 01451 #===================================================================== 01452 # 01453 # Now save stat dictionaries using Pickle 01454 pickfile = prefix + '.pickle' 01455 f = open(pickfile,'w') 01456 p = pickle.Pickler(f) 01457 p.dump(clnmodel) 01458 p.dump(polmodel) 01459 f.close() 01460 01461 print "" 01462 print "Dictionaries clnmodel,polmodel saved in "+pickfile 01463 print "Use Pickle to retrieve" 01464 print "" 01465 01466 # e.g. 01467 # f = open(pickfile) 01468 # u = Unpickler(f) 01469 # clnmodel = u.load() 01470 # polmodel = u.load() 01471 # f.close() 01472 01473 print "" 01474 print "Completed Processing"
1.8.0
