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