|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ###################################################################### 00002 # # 00003 # Use Case Script for Jupiter 6cm VLA # 00004 # # 00005 # Last Updated STM 2007-10-10 (Beta) # 00006 # # 00007 ###################################################################### 00008 00009 import time 00010 import os 00011 00012 # 00013 #===================================================================== 00014 # 00015 # This script has some interactive commands: scriptmode = True 00016 # if you are running it and want it to stop during interactive parts. 00017 00018 scriptmode = True 00019 00020 #===================================================================== 00021 # 00022 # Set up some useful variables - these will be set during the script 00023 # also, but if you want to restart the script in the middle here 00024 # they are in one place: 00025 00026 pathname=os.environ.get('CASAPATH').split()[0] 00027 prefix='jupiter6cm.usecase' 00028 00029 msfile = prefix + '.ms' 00030 00031 gtable = prefix + '.gcal' 00032 ftable = prefix + '.fluxscale' 00033 atable = prefix + '.accum' 00034 00035 srcsplitms = prefix + '.split.ms' 00036 00037 clnimsize = [288,288] 00038 clncell = [4.,4.] 00039 00040 imname1 = prefix + '.clean1' 00041 clnimage1 = imname1+'.image' 00042 clnmodel1 = imname1+'.model' 00043 clnresid1 = imname1+'.residual' 00044 clnmask1 = imname1+'.clean_interactive.mask' 00045 00046 selfcaltab1 = srcsplitms + '.selfcal1' 00047 00048 imname2 = prefix + '.clean2' 00049 clnimage2 = imname2+'.image' 00050 clnmodel2 = imname2+'.model' 00051 clnresid2 = imname2+'.residual' 00052 clnmask2 = imname2+'.clean_interactive.mask' 00053 00054 selfcaltab2 = srcsplitms + '.selfcal2' 00055 smoothcaltab2 = srcsplitms + '.smoothcal2' 00056 00057 imname3 = prefix + '.clean3' 00058 clnimage3 = imname3+'.image' 00059 clnmodel3 = imname3+'.model' 00060 clnresid3 = imname3+'.residual' 00061 clnmask3 = imname3+'.clean_interactive.mask' 00062 00063 # 00064 #===================================================================== 00065 # 00066 # Get to path to the CASA home and stip off the name 00067 pathname=os.environ.get('CASAPATH').split()[0] 00068 00069 # This is where the UVFITS data will be 00070 #fitsdata=pathname+'/data/demo/jupiter6cm.fits' 00071 fitsdata='/home/sandrock2/smyers/NAUG2/Data/VLA_CONT/FLUX99-6CM.CBAND' 00072 00073 # The prefix to use for all output files 00074 prefix='jupiter6cm.usecase' 00075 00076 # Clean up old files 00077 os.system('rm -rf '+prefix+'*') 00078 00079 # 00080 #===================================================================== 00081 # Data Import and List 00082 #===================================================================== 00083 # 00084 # Import the data from FITS to MS 00085 # 00086 print '--Import--' 00087 00088 # Safest to start from task defaults 00089 default('importuvfits') 00090 00091 # Set up the MS filename and save as new global variable 00092 msfile = prefix + '.ms' 00093 00094 # Use task importuvfits 00095 fitsfile = fitsdata 00096 vis = msfile 00097 importuvfits() 00098 00099 #===================================================================== 00100 # 00101 # List a summary of the MS 00102 # 00103 print '--Listobs--' 00104 00105 # Don't default this one and make use of the previous setting of 00106 # vis. Remember, the variables are GLOBAL! 00107 00108 # You may wish to see more detailed information, in this case 00109 # use the verbose = True option 00110 verbose = True 00111 00112 listobs() 00113 00114 # You should get in your logger window and in the casapy.log file 00115 # something like: 00116 # 00117 # Observer: FLUX99 Project: 00118 # Observation: VLA 00119 # 00120 # Data records: 2021424 Total integration time = 85133.2 seconds 00121 # Observed from 23:15:27 to 22:54:20 00122 # 00123 # ObservationID = 0 ArrayID = 0 00124 # Date Timerange Scan FldId FieldName SpwIds 00125 # 15-Apr-1999/23:15:26.7 - 23:16:10.0 1 0 0137+331 [0, 1] 00126 # 23:38:40.0 - 23:48:00.0 2 1 0813+482 [0, 1] 00127 # 23:53:40.0 - 23:55:20.0 3 2 0542+498 [0, 1] 00128 # 16-Apr-1999/00:22:10.1 - 00:23:49.9 4 3 0437+296 [0, 1] 00129 # 00:28:23.3 - 00:30:00.1 5 4 VENUS [0, 1] 00130 # 00:48:40.0 - 00:50:20.0 6 1 0813+482 [0, 1] 00131 # 00:56:13.4 - 00:57:49.9 7 2 0542+498 [0, 1] 00132 # 01:10:20.1 - 01:11:59.9 8 5 0521+166 [0, 1] 00133 # 01:23:29.9 - 01:25:00.1 9 3 0437+296 [0, 1] 00134 # 01:29:33.3 - 01:31:10.0 10 4 VENUS [0, 1] 00135 # 01:49:50.0 - 01:51:30.0 11 6 1411+522 [0, 1] 00136 # 02:03:00.0 - 02:04:30.0 12 7 1331+305 [0, 1] 00137 # 02:17:30.0 - 02:19:10.0 13 1 0813+482 [0, 1] 00138 # 02:24:20.0 - 02:26:00.0 14 2 0542+498 [0, 1] 00139 # 02:37:49.9 - 02:39:30.0 15 5 0521+166 [0, 1] 00140 # 02:50:50.1 - 02:52:20.1 16 3 0437+296 [0, 1] 00141 # 02:59:20.0 - 03:01:00.0 17 6 1411+522 [0, 1] 00142 # 03:12:30.0 - 03:14:10.0 18 7 1331+305 [0, 1] 00143 # 03:27:53.3 - 03:29:39.9 19 1 0813+482 [0, 1] 00144 # 03:35:00.0 - 03:36:40.0 20 2 0542+498 [0, 1] 00145 # 03:49:50.0 - 03:51:30.1 21 6 1411+522 [0, 1] 00146 # 04:03:10.0 - 04:04:50.0 22 7 1331+305 [0, 1] 00147 # 04:18:49.9 - 04:20:40.0 23 1 0813+482 [0, 1] 00148 # 04:25:56.6 - 04:27:39.9 24 2 0542+498 [0, 1] 00149 # 04:42:49.9 - 04:44:40.0 25 8 MARS [0, 1] 00150 # 04:56:50.0 - 04:58:30.1 26 6 1411+522 [0, 1] 00151 # 05:24:03.3 - 05:33:39.9 27 7 1331+305 [0, 1] 00152 # 05:48:00.0 - 05:49:49.9 28 1 0813+482 [0, 1] 00153 # 05:58:36.6 - 06:00:30.0 29 8 MARS [0, 1] 00154 # 06:13:20.1 - 06:14:59.9 30 6 1411+522 [0, 1] 00155 # 06:27:40.0 - 06:29:20.0 31 7 1331+305 [0, 1] 00156 # 06:44:13.4 - 06:46:00.0 32 1 0813+482 [0, 1] 00157 # 06:55:06.6 - 06:57:00.0 33 8 MARS [0, 1] 00158 # 07:10:40.0 - 07:12:20.0 34 6 1411+522 [0, 1] 00159 # 07:28:20.0 - 07:30:10.1 35 7 1331+305 [0, 1] 00160 # 07:42:49.9 - 07:44:30.0 36 8 MARS [0, 1] 00161 # 07:58:43.3 - 08:00:39.9 37 6 1411+522 [0, 1] 00162 # 08:13:30.0 - 08:15:19.9 38 7 1331+305 [0, 1] 00163 # 08:27:53.4 - 08:29:30.0 39 8 MARS [0, 1] 00164 # 08:42:59.9 - 08:44:50.0 40 6 1411+522 [0, 1] 00165 # 08:57:09.9 - 08:58:50.0 41 7 1331+305 [0, 1] 00166 # 09:13:03.3 - 09:14:50.1 42 9 NGC7027 [0, 1] 00167 # 09:26:59.9 - 09:28:40.0 43 6 1411+522 [0, 1] 00168 # 09:40:33.4 - 09:42:09.9 44 7 1331+305 [0, 1] 00169 # 09:56:19.9 - 09:58:10.0 45 9 NGC7027 [0, 1] 00170 # 10:12:59.9 - 10:14:50.0 46 8 MARS [0, 1] 00171 # 10:27:09.9 - 10:28:50.0 47 6 1411+522 [0, 1] 00172 # 10:40:30.0 - 10:42:00.0 48 7 1331+305 [0, 1] 00173 # 10:56:10.0 - 10:57:50.0 49 9 NGC7027 [0, 1] 00174 # 11:28:30.0 - 11:35:30.0 50 10 NEPTUNE [0, 1] 00175 # 11:48:20.0 - 11:50:10.0 51 6 1411+522 [0, 1] 00176 # 12:01:36.7 - 12:03:10.0 52 7 1331+305 [0, 1] 00177 # 12:35:33.3 - 12:37:40.0 53 11 URANUS [0, 1] 00178 # 12:46:30.0 - 12:48:10.0 54 10 NEPTUNE [0, 1] 00179 # 13:00:29.9 - 13:02:10.0 55 6 1411+522 [0, 1] 00180 # 13:15:23.3 - 13:17:10.1 56 9 NGC7027 [0, 1] 00181 # 13:33:43.3 - 13:35:40.0 57 11 URANUS [0, 1] 00182 # 13:44:30.0 - 13:46:10.0 58 10 NEPTUNE [0, 1] 00183 # 14:00:46.7 - 14:01:39.9 59 0 0137+331 [0, 1] 00184 # 14:10:40.0 - 14:12:09.9 60 12 JUPITER [0, 1] 00185 # 14:24:06.6 - 14:25:40.1 61 11 URANUS [0, 1] 00186 # 14:34:30.0 - 14:36:10.1 62 10 NEPTUNE [0, 1] 00187 # 14:59:13.4 - 15:00:00.0 63 0 0137+331 [0, 1] 00188 # 15:09:03.3 - 15:10:40.1 64 12 JUPITER [0, 1] 00189 # 15:24:30.0 - 15:26:20.1 65 9 NGC7027 [0, 1] 00190 # 15:40:10.0 - 15:45:00.0 66 11 URANUS [0, 1] 00191 # 15:53:50.0 - 15:55:20.0 67 10 NEPTUNE [0, 1] 00192 # 16:18:53.4 - 16:19:49.9 68 0 0137+331 [0, 1] 00193 # 16:29:10.1 - 16:30:49.9 69 12 JUPITER [0, 1] 00194 # 16:42:53.4 - 16:44:30.0 70 11 URANUS [0, 1] 00195 # 16:54:53.4 - 16:56:40.0 71 9 NGC7027 [0, 1] 00196 # 17:23:06.6 - 17:30:40.0 72 2 0542+498 [0, 1] 00197 # 17:41:50.0 - 17:43:20.0 73 3 0437+296 [0, 1] 00198 # 17:55:36.7 - 17:57:39.9 74 4 VENUS [0, 1] 00199 # 18:19:23.3 - 18:20:09.9 75 0 0137+331 [0, 1] 00200 # 18:30:23.3 - 18:32:00.0 76 12 JUPITER [0, 1] 00201 # 18:44:49.9 - 18:46:30.0 77 9 NGC7027 [0, 1] 00202 # 18:59:13.3 - 19:00:59.9 78 2 0542+498 [0, 1] 00203 # 19:19:10.0 - 19:21:20.1 79 5 0521+166 [0, 1] 00204 # 19:32:50.1 - 19:34:29.9 80 3 0437+296 [0, 1] 00205 # 19:39:03.3 - 19:40:40.1 81 4 VENUS [0, 1] 00206 # 20:08:06.7 - 20:08:59.9 82 0 0137+331 [0, 1] 00207 # 20:18:10.0 - 20:19:50.0 83 12 JUPITER [0, 1] 00208 # 20:33:53.3 - 20:35:40.1 84 1 0813+482 [0, 1] 00209 # 20:40:59.9 - 20:42:40.0 85 2 0542+498 [0, 1] 00210 # 21:00:16.6 - 21:02:20.1 86 5 0521+166 [0, 1] 00211 # 21:13:53.4 - 21:15:29.9 87 3 0437+296 [0, 1] 00212 # 21:20:43.4 - 21:22:30.0 88 4 VENUS [0, 1] 00213 # 21:47:26.7 - 21:48:20.1 89 0 0137+331 [0, 1] 00214 # 21:57:30.0 - 21:59:10.0 90 12 JUPITER [0, 1] 00215 # 22:12:13.3 - 22:14:00.1 91 2 0542+498 [0, 1] 00216 # 22:28:33.3 - 22:30:19.9 92 4 VENUS [0, 1] 00217 # 22:53:33.3 - 22:54:19.9 93 0 0137+331 [0, 1] 00218 # 00219 # Fields: 13 00220 # ID Name Right Ascension Declination Epoch 00221 # 0 0137+331 01:37:41.30 +33.09.35.13 J2000 00222 # 1 0813+482 08:13:36.05 +48.13.02.26 J2000 00223 # 2 0542+498 05:42:36.14 +49.51.07.23 J2000 00224 # 3 0437+296 04:37:04.17 +29.40.15.14 J2000 00225 # 4 VENUS 04:06:54.11 +22.30.35.91 J2000 00226 # 5 0521+166 05:21:09.89 +16.38.22.05 J2000 00227 # 6 1411+522 14:11:20.65 +52.12.09.14 J2000 00228 # 7 1331+305 13:31:08.29 +30.30.32.96 J2000 00229 # 8 MARS 14:21:41.37 -12.21.49.45 J2000 00230 # 9 NGC7027 21:07:01.59 +42.14.10.19 J2000 00231 # 10 NEPTUNE 20:26:01.14 -18.54.54.21 J2000 00232 # 11 URANUS 21:15:42.83 -16.35.05.59 J2000 00233 # 12 JUPITER 00:55:34.04 +04.45.44.71 J2000 00234 # 00235 # Spectral Windows: (2 unique spectral windows and 1 unique polarization setups) 00236 # SpwID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs 00237 # 0 1 TOPO 4885.1 50000 50000 4885.1 RR RL LR LL 00238 # 1 1 TOPO 4835.1 50000 50000 4835.1 RR RL LR LL 00239 # 00240 # Feeds: 28: printing first row only 00241 # Antenna Spectral Window # Receptors Polarizations 00242 # 1 -1 2 [ R, L] 00243 # 00244 # Antennas: 27: 00245 # ID Name Station Diam. Long. Lat. 00246 # 0 1 VLA:W9 25.0 m -107.37.25.1 +33.53.51.0 00247 # 1 2 VLA:N9 25.0 m -107.37.07.8 +33.54.19.0 00248 # 2 3 VLA:N3 25.0 m -107.37.06.3 +33.54.04.8 00249 # 3 4 VLA:N5 25.0 m -107.37.06.7 +33.54.08.0 00250 # 4 5 VLA:N2 25.0 m -107.37.06.2 +33.54.03.5 00251 # 5 6 VLA:E1 25.0 m -107.37.05.7 +33.53.59.2 00252 # 6 7 VLA:E2 25.0 m -107.37.04.4 +33.54.01.1 00253 # 7 8 VLA:N8 25.0 m -107.37.07.5 +33.54.15.8 00254 # 8 9 VLA:E8 25.0 m -107.36.48.9 +33.53.55.1 00255 # 9 10 VLA:W3 25.0 m -107.37.08.9 +33.54.00.1 00256 # 10 11 VLA:N1 25.0 m -107.37.06.0 +33.54.01.8 00257 # 11 12 VLA:E6 25.0 m -107.36.55.6 +33.53.57.7 00258 # 12 13 VLA:W7 25.0 m -107.37.18.4 +33.53.54.8 00259 # 13 14 VLA:E4 25.0 m -107.37.00.8 +33.53.59.7 00260 # 14 15 VLA:N7 25.0 m -107.37.07.2 +33.54.12.9 00261 # 15 16 VLA:W4 25.0 m -107.37.10.8 +33.53.59.1 00262 # 16 17 VLA:W5 25.0 m -107.37.13.0 +33.53.57.8 00263 # 17 18 VLA:N6 25.0 m -107.37.06.9 +33.54.10.3 00264 # 18 19 VLA:E7 25.0 m -107.36.52.4 +33.53.56.5 00265 # 19 20 VLA:E9 25.0 m -107.36.45.1 +33.53.53.6 00266 # 21 22 VLA:W8 25.0 m -107.37.21.6 +33.53.53.0 00267 # 22 23 VLA:W6 25.0 m -107.37.15.6 +33.53.56.4 00268 # 23 24 VLA:W1 25.0 m -107.37.05.9 +33.54.00.5 00269 # 24 25 VLA:W2 25.0 m -107.37.07.4 +33.54.00.9 00270 # 25 26 VLA:E5 25.0 m -107.36.58.4 +33.53.58.8 00271 # 26 27 VLA:N4 25.0 m -107.37.06.5 +33.54.06.1 00272 # 27 28 VLA:E3 25.0 m -107.37.02.8 +33.54.00.5 00273 # 00274 # Tables: 00275 # MAIN 2021424 rows 00276 # ANTENNA 28 rows 00277 # DATA_DESCRIPTION 2 rows 00278 # DOPPLER <absent> 00279 # FEED 28 rows 00280 # FIELD 13 rows 00281 # FLAG_CMD <empty> 00282 # FREQ_OFFSET <absent> 00283 # HISTORY 7058 rows 00284 # OBSERVATION 1 row 00285 # POINTING 2604 rows 00286 # POLARIZATION 1 row 00287 # PROCESSOR <empty> 00288 # SOURCE <empty> (see FIELD) 00289 # SPECTRAL_WINDOW 2 rows 00290 # STATE <empty> 00291 # SYSCAL <absent> 00292 # WEATHER <absent> 00293 00294 # 00295 #===================================================================== 00296 # Data Examination and Flagging 00297 #===================================================================== 00298 # 00299 # Get rid of the autocorrelations from the MS 00300 # 00301 print '--Flag auto-correlations--' 00302 00303 default(tflagdata) 00304 vis = msfile 00305 mode = 'manual' 00306 autocorr = True 00307 tflagdata() 00308 00309 # 00310 #===================================================================== 00311 # 00312 # Use Flagmanager to save a copy of the flags 00313 # 00314 print '--Flagmanager--' 00315 default('flagmanager') 00316 00317 vis = msfile 00318 00319 # Save a copy of the MAIN table flags 00320 00321 mode = 'save' 00322 versionname = 'flagautocorr' 00323 comment = 'flagged autocorr' 00324 merge = 'replace' 00325 00326 flagmanager() 00327 00328 # If you look in the 'jupiter6cm.usecase.ms.flagversions/ 00329 # you'll see flags.flagautocorr there along with the 00330 # flags.Original that importuvfits made for you 00331 # Or use 00332 00333 mode = 'list' 00334 00335 flagmanager() 00336 00337 # In the logger you will see something like: 00338 # 00339 # MS : /home/sandrock2/smyers/Testing2/Aug07/jupiter6cm.usecase.ms 00340 # 00341 # main : working copy in main table 00342 # Original : Original flags at import into CASA 00343 # flagautocorr : flagged autocorr 00344 # See logger for flag versions for this file 00345 00346 # 00347 #===================================================================== 00348 # 00349 # Use Plotxy to interactively flag the data 00350 # 00351 print '--Plotxy--' 00352 default('plotxy') 00353 00354 vis = msfile 00355 00356 # The fields we are interested in: 1331+305,JUPITER,0137+331 00357 selectdata = True 00358 00359 # First we do the primary calibrator 00360 field = '1331+305' 00361 00362 # Plot only the RR and LL for now 00363 correlation = 'RR LL' 00364 00365 # Plot amplitude vs. uvdist 00366 xaxis = 'uvdist' 00367 yaxis = 'amp' 00368 multicolor = 'both' 00369 00370 # The easiest thing is to iterate over antennas 00371 iteration = 'antenna' 00372 00373 plotxy() 00374 00375 # Pause script if you are running in scriptmode 00376 if scriptmode: 00377 user_check=raw_input('Return to continue script\n') 00378 00379 # You'll see lots of low points as you step through RR LL RL LR 00380 # A basic clip at 0.75 for RR LL and 0.055 for RL LR will work 00381 # If you want to do this interactively, set 00382 iteration = '' 00383 00384 plotxy() 00385 00386 # You can also use flagdata to do this non-interactively 00387 # (see below) 00388 00389 # Now look at the cross-polar products 00390 correlation = 'RL LR' 00391 00392 plotxy() 00393 00394 # Pause script if you are running in scriptmode 00395 if scriptmode: 00396 user_check=raw_input('Return to continue script\n') 00397 00398 #--------------------------------------------------------------------- 00399 # Now do calibrater 0137+331 00400 field = '0137+331' 00401 correlation = 'RR LL' 00402 xaxis = 'uvdist' 00403 spw = '' 00404 iteration = '' 00405 antenna = '' 00406 00407 plotxy() 00408 00409 # You'll see a bunch of bad data along the bottom near zero amp 00410 # Draw a box around some of it and use Locate 00411 # Looks like much of it is Antenna 9 (ID=8) in spw=1 00412 00413 # Pause script if you are running in scriptmode 00414 if scriptmode: 00415 user_check=raw_input('Return to continue script\n') 00416 00417 xaxis = 'time' 00418 spw = '1' 00419 correlation = '' 00420 00421 # Note that the strings like antenna='9' first try to match the 00422 # NAME which we see in listobs was the number '9' for ID=8. 00423 # So be careful here (why naming antennas as numbers is bad). 00424 antenna = '9' 00425 00426 plotxy() 00427 00428 # YES! the last 4 scans are bad. Box 'em and flag. 00429 00430 # Pause script if you are running in scriptmode 00431 if scriptmode: 00432 user_check=raw_input('Return to continue script\n') 00433 00434 # Go back and clean up 00435 xaxis = 'uvdist' 00436 spw = '' 00437 antenna = '' 00438 correlation = 'RR LL' 00439 00440 plotxy() 00441 00442 # Box up the bad low points (basically a clip below 0.52) and flag 00443 00444 # Note that RL,LR are too weak to clip on. 00445 00446 # Pause script if you are running in scriptmode 00447 if scriptmode: 00448 user_check=raw_input('Return to continue script\n') 00449 00450 #--------------------------------------------------------------------- 00451 # Finally, do JUPITER 00452 field = 'JUPITER' 00453 correlation = '' 00454 iteration = '' 00455 xaxis = 'time' 00456 00457 plotxy() 00458 00459 # Here you will see that the final scan at 22:00:00 UT is bad 00460 # Draw a box around it and flag it! 00461 00462 # Pause script if you are running in scriptmode 00463 if scriptmode: 00464 user_check=raw_input('Return to continue script\n') 00465 00466 # Now look at whats left 00467 correlation = 'RR LL' 00468 xaxis = 'uvdist' 00469 spw = '1' 00470 antenna = '' 00471 iteration = 'antenna' 00472 00473 plotxy() 00474 00475 # As you step through, you will see that Antenna 9 (ID=8) is often 00476 # bad in this spw. If you box and do Locate (or remember from 00477 # 0137+331) its probably a bad time. 00478 00479 # Pause script if you are running in scriptmode 00480 if scriptmode: 00481 user_check=raw_input('Return to continue script\n') 00482 00483 # The easiset way to kill it: 00484 00485 antenna = '9' 00486 iteration = '' 00487 xaxis = 'time' 00488 correlation = '' 00489 00490 plotxy() 00491 00492 # Draw a box around all points in the last bad scans and flag 'em! 00493 00494 # Pause script if you are running in scriptmode 00495 if scriptmode: 00496 user_check=raw_input('Return to continue script\n') 00497 00498 # Now clean up the rest 00499 xaxis = 'uvdist' 00500 correlation = 'RR LL' 00501 antenna = '' 00502 spw = '' 00503 00504 # You will be drawing many tiny boxes, so remember you can 00505 # use the ESC key to get rid of the most recent box if you 00506 # make a mistake. 00507 00508 plotxy() 00509 00510 # Note that the end result is we've flagged lots of points 00511 # in RR and LL. We will rely upon imager to ignore the 00512 # RL LR for points with RR LL flagged! 00513 00514 # Pause script if you are running in scriptmode 00515 if scriptmode: 00516 user_check=raw_input('Return to continue script\n') 00517 00518 # 00519 #===================================================================== 00520 # 00521 # Use Flagmanager to save a copy of the flags so far 00522 # 00523 print '--Flagmanager--' 00524 default('flagmanager') 00525 00526 vis = msfile 00527 mode = 'save' 00528 versionname = 'xyflags' 00529 comment = 'Plotxy flags' 00530 merge = 'replace' 00531 00532 flagmanager() 00533 00534 # 00535 #===================================================================== 00536 # 00537 # You can use Flagdata to explicitly clip the data also 00538 # 00539 print '--Flagdata--' 00540 default('tflagdata') 00541 00542 vis = msfile 00543 00544 00545 # Set some clipping regions 00546 mode = 'clip' 00547 clipcolumn = 'DATA' 00548 clipoutside = False 00549 00550 # Clip calibraters 00551 field = '1331+305' 00552 correlation = 'ABS_RR,LL' 00553 clipminmax = [0.0,0.75] 00554 tflagdata() 00555 00556 correlation = 'ABS_RL,LR' 00557 clipminmax = [0.0,0.055] 00558 tflagdata() 00559 00560 field = '0137+331' 00561 correlation = 'ABS_RR,LL' 00562 clipminmax = [0.0,0.55] 00563 tflagdata() 00564 00565 # You can also do the antenna edits on 0137+331 and JUPITER 00566 # with flagdata 00567 00568 # 00569 #===================================================================== 00570 # Calibration 00571 #===================================================================== 00572 # 00573 # Set the fluxes of the primary calibrator(s) 00574 # 00575 print '--Setjy--' 00576 default('setjy') 00577 00578 vis = msfile 00579 00580 # 00581 # 1331+305 = 3C286 is our primary calibrator 00582 field = '1331+305' 00583 00584 # Setjy knows about this source so we dont need anything more 00585 00586 setjy() 00587 00588 # 00589 # You should see something like this in the logger and casapy.log file: 00590 # 00591 # 1331+305 spwid= 0 [I=7.462, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00592 # 1331+305 spwid= 1 [I=7.51, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00593 # 00594 00595 # 00596 #===================================================================== 00597 # 00598 # Initial gain calibration 00599 # 00600 print '--Gaincal--' 00601 default('gaincal') 00602 00603 vis = msfile 00604 00605 # set the name for the output gain caltable 00606 gtable = prefix + '.gcal' 00607 caltable = gtable 00608 00609 # Gain calibrators are 1331+305 and 0137+331 (FIELD_ID 7 and 0) 00610 # We have 2 IFs (SPW 0,1) with one channel each 00611 00612 # selection is via the field and spw strings 00613 field = '1331+305,0137+331' 00614 spw = '' 00615 00616 # a-priori calibration application 00617 # atmospheric optical depth (turn off) 00618 gaincurve = True 00619 opacity = 0.0 00620 00621 # scan-based G solutions for both amplitude and phase 00622 gaintype = 'G' 00623 solint = 0. 00624 calmode = 'ap' 00625 00626 # reference antenna 11 (11=VLA:N1) 00627 refant = '11' 00628 00629 # minimum SNR 3 00630 minsnr = 3 00631 00632 gaincal() 00633 00634 # 00635 #===================================================================== 00636 # 00637 # Bootstrap flux scale 00638 # 00639 print '--Fluxscale--' 00640 default('fluxscale') 00641 00642 vis = msfile 00643 00644 # set the name for the output rescaled caltable 00645 ftable = prefix + '.fluxscale' 00646 fluxtable = ftable 00647 00648 # point to our first gain cal table 00649 caltable = gtable 00650 00651 # we will be using 1331+305 (the source we did setjy on) as 00652 # our flux standard reference 00653 reference = '1331+305' 00654 00655 # we want to transfer the flux to our other gain cal source 0137+331 00656 # to bring its gain amplitues in line with the absolute scale 00657 transfer = '0137+331' 00658 00659 fluxscale() 00660 00661 # You should see in the logger something like: 00662 #Flux density for 0137+331 in SpW=0 is: 00663 # 5.42575 +/- 0.00285011 (SNR = 1903.7, nAnt= 27) 00664 #Flux density for 0137+331 in SpW=1 is: 00665 # 5.46569 +/- 0.00301326 (SNR = 1813.88, nAnt= 27) 00666 00667 #===================================================================== 00668 # 00669 # Interpolate the gains onto Jupiter (and others) 00670 # 00671 print '--Accum--' 00672 default('accum') 00673 00674 vis = msfile 00675 00676 tablein = '' 00677 incrtable = ftable 00678 calfield = '1331+305, 0137+331' 00679 00680 # set the name for the output interpolated caltable 00681 atable = prefix + '.accum' 00682 caltable = atable 00683 00684 # linear interpolation 00685 interp = 'linear' 00686 00687 # make 10s entries 00688 accumtime = 10.0 00689 00690 accum() 00691 00692 #===================================================================== 00693 # 00694 # Correct the data 00695 # (This will put calibrated data into the CORRECTED_DATA column) 00696 # 00697 print '--ApplyCal--' 00698 default('applycal') 00699 00700 vis = msfile 00701 00702 # Start with the interpolated fluxscale/gain table 00703 bptable = '' 00704 gaintable = atable 00705 00706 # Since we did gaincurve=True in gaincal, we need it here also 00707 gaincurve = True 00708 opacity=0.0 00709 00710 # select the fields 00711 field = '1331+305,0137+331,JUPITER' 00712 spw = '' 00713 selectdata = False 00714 00715 # do not need to select subset since we did accum 00716 # (note that correct only does 'nearest' interp) 00717 gainselect = '' 00718 00719 applycal() 00720 00721 # 00722 #===================================================================== 00723 # 00724 # Now split the Jupiter target data 00725 # 00726 print '--Split Jupiter--' 00727 default('split') 00728 00729 vis = msfile 00730 00731 # Now we write out the corrected data for the calibrator 00732 00733 # Make an output vis file 00734 srcsplitms = prefix + '.split.ms' 00735 outputvis = srcsplitms 00736 00737 # Select the Jupiter field 00738 field = 'JUPITER' 00739 spw = '' 00740 00741 # pick off the CORRECTED_DATA column 00742 datacolumn = 'corrected' 00743 00744 split() 00745 00746 #===================================================================== 00747 # 00748 # Export the Jupiter data as UVFITS 00749 # Start with the split file. 00750 # 00751 print '--Export UVFITS--' 00752 default('exportuvfits') 00753 00754 srcuvfits = prefix + '.split.uvfits' 00755 00756 vis = srcsplitms 00757 fitsfile = srcuvfits 00758 00759 # Since this is a split dataset, the calibrated data is 00760 # in the DATA column already. 00761 datacolumn = 'data' 00762 00763 # Write as a multisource UVFITS (with SU table) 00764 # even though it will have only one field in it 00765 multisource = True 00766 00767 # Run asynchronously so as not to interfere with other tasks 00768 # (BETA: also avoids crash on next importuvfits) 00769 async = True 00770 00771 exportuvfits() 00772 00773 # 00774 #===================================================================== 00775 # FIRST CLEAN / SELFCAL CYCLE 00776 #===================================================================== 00777 # 00778 # Now clean an image of Jupiter 00779 # 00780 print '--Clean 1--' 00781 default('clean') 00782 00783 # Pick up our split source data 00784 vis = srcsplitms 00785 00786 # Make an image root file name 00787 imname1 = prefix + '.clean1' 00788 imagename = imname1 00789 00790 # Set up the output continuum image (single plane mfs) 00791 mode = 'mfs' 00792 stokes = 'I' 00793 00794 # NOTE: current version field='' doesnt work 00795 field = '*' 00796 00797 # Combine all spw 00798 spw = '' 00799 00800 # This is D-config VLA 6cm (4.85GHz) obs 00801 # Check the observational status summary 00802 # Primary beam FWHM = 45'/f_GHz = 557" 00803 # Synthesized beam FWHM = 14" 00804 # RMS in 10min (600s) = 0.06 mJy (thats now, but close enough) 00805 00806 # Set the output image size and cell size (arcsec) 00807 # 4" will give 3.5x oversampling 00808 # 280 pix will cover to 2xPrimaryBeam 00809 # clean will say to use 288 (a composite integer) for efficiency 00810 clnalg = 'clark' 00811 clnimsize = [288,288] 00812 00813 # double for CS Clean 00814 #clnalg = 'csclean' 00815 #clnimsize = [576,576] 00816 00817 clncell = [4.,4.] 00818 00819 alg = clnalg 00820 imsize = clnimsize 00821 cell = clncell 00822 00823 # NOTE: will eventually have an imadvise task to give you this 00824 # information 00825 00826 # Standard gain factor 0.1 00827 gain = 0.1 00828 00829 # Fix maximum number of iterations 00830 niter = 10000 00831 00832 # Also set flux residual threshold (0.04 mJy) 00833 # From our listobs: 00834 # Total integration time = 85133.2 seconds 00835 # With rms of 0.06 mJy in 600s ==> rms = 0.005 mJy 00836 # Set to 10x thermal rms 00837 threshold=0.05 00838 00839 # Note - we can change niter and threshold interactively 00840 # during clean 00841 00842 # Set up the weighting 00843 # Use Briggs weighting (a moderate value, on the uniform side) 00844 weighting = 'briggs' 00845 rmode = 'norm' 00846 robust = 0.5 00847 00848 # No clean mask 00849 mask = '' 00850 00851 # Use interactive clean mode 00852 cleanbox = 'interactive' 00853 00854 # Moderate number of iter per interactive cycle 00855 npercycle = 100 00856 00857 clean() 00858 00859 # When the interactive clean window comes up, use the right-mouse 00860 # to draw rectangles around obvious emission double-right-clicking 00861 # inside them to add to the flag region. You can also assign the 00862 # right-mouse to polygon region drawing by right-clicking on the 00863 # polygon drawing icon in the toolbar. When you are happy with 00864 # the region, click 'Done Flagging' and it will go and clean another 00865 # 100 iterations. When done, click 'Stop'. 00866 00867 # Set up variables 00868 clnimage1 = imname1+'.image' 00869 clnmodel1 = imname1+'.model' 00870 clnresid1 = imname1+'.residual' 00871 clnmask1 = imname1+'.clean_interactive.mask' 00872 00873 # 00874 #--------------------------------------------------------------------- 00875 # 00876 # Look at this in viewer 00877 viewer(clnimage1,'image') 00878 00879 # You can use the right-mouse to draw a box in the lower right 00880 # corner of the image away from emission, the double-click inside 00881 # to bring up statistics. Use the right-mouse to grab this box 00882 # and move it up over Jupiter and double-click again. You should 00883 # see stuff like this in the terminal: 00884 # 00885 # jupiter6cm.usecase.clean1.image (Jy/beam) 00886 # 00887 # n Std Dev RMS Mean Variance Sum 00888 # 4712 0.003914 0.003927 0.0003205 1.532e-05 1.510 00889 # 00890 # Flux Med |Dev| IntQtlRng Median Min Max 00891 # 0.09417 0.002646 0.005294 0.0001885 -0.01125 0.01503 00892 # 00893 # 00894 # On Jupiter: 00895 # 00896 # n Std Dev RMS Mean Variance Sum 00897 # 3640 0.1007 0.1027 0.02023 0.01015 73.63 00898 # 00899 # Flux Med |Dev| IntQtlRng Median Min Max 00900 # 4.592 0.003239 0.007120 0.0001329 -0.01396 1.060 00901 # 00902 # Estimated dynamic range = 1.060 / 0.003927 = 270 (poor) 00903 # 00904 # Note that the exact numbers you get will depend on how deep you 00905 # take the interactive clean and how you draw the box for the stats. 00906 # 00907 #--------------------------------------------------------------------- 00908 # 00909 # Self-cal using clean model 00910 # 00911 # Note: clean will have left FT of model in the MODEL_DATA column 00912 # If you've done something in between, can use the ft task to 00913 # do this manually. 00914 # 00915 print '--SelfCal 1--' 00916 default('gaincal') 00917 00918 vis = srcsplitms 00919 00920 # New gain table 00921 selfcaltab1 = srcsplitms + '.selfcal1' 00922 caltable = selfcaltab1 00923 00924 # Don't need a-priori cals 00925 selectdata = False 00926 gaincurve = False 00927 opacity = 0.0 00928 00929 # This choice seemed to work 00930 refant = '11' 00931 00932 # Lets do phase-only first time around 00933 gaintype = 'G' 00934 calmode = 'p' 00935 00936 # Do scan-based solutions with SNR>3 00937 solint = 0.0 00938 minsnr = 3.0 00939 00940 # Do not need to normalize (let gains float) 00941 solnorm = False 00942 00943 gaincal() 00944 00945 # 00946 #--------------------------------------------------------------------- 00947 # 00948 # Correct the data (no need for interpolation this stage) 00949 # 00950 print '--ApplyCal--' 00951 default('applycal') 00952 00953 vis = srcsplitms 00954 00955 gaintable = selfcaltab1 00956 00957 gaincurve = False 00958 opacity = 0.0 00959 field = '' 00960 spw = '' 00961 selectdata = False 00962 00963 calwt = True 00964 00965 applycal() 00966 00967 # Self-cal is now in CORRECTED_DATA column of split ms 00968 # 00969 #===================================================================== 00970 # SECOND CLEAN / SELFCAL CYCLE 00971 #===================================================================== 00972 # 00973 print '--Clean 2--' 00974 default('clean') 00975 00976 vis = srcsplitms 00977 00978 imname2 = prefix + '.clean2' 00979 imagename = imname2 00980 00981 field = '*' 00982 spw = '' 00983 mode = 'mfs' 00984 gain = 0.1 00985 niter = 10000 00986 threshold=0.04 00987 00988 alg = clnalg 00989 imsize = clnimsize 00990 cell = clncell 00991 00992 weighting = 'briggs' 00993 rmode = 'norm' 00994 robust = 0.5 00995 00996 cleanbox = 'interactive' 00997 npercycle = 100 00998 00999 clean() 01000 01001 # Set up variables 01002 clnimage2 = imname2+'.image' 01003 clnmodel2 = imname2+'.model' 01004 clnresid2 = imname2+'.residual' 01005 clnmask2 = imname2+'.clean_interactive.mask' 01006 01007 # 01008 #--------------------------------------------------------------------- 01009 # 01010 # Look at this in viewer 01011 viewer(clnimage2,'image') 01012 01013 # jupiter6cm.usecase.clean2.image (Jy/beam) 01014 # 01015 # n Std Dev RMS Mean Variance Sum 01016 # 5236 0.001389 0.001390 3.244e-05 1.930e-06 0.1699 01017 # 01018 # Flux Med |Dev| IntQtlRng Median Min Max 01019 # 0.01060 0.0009064 0.001823 -1.884e-05 -0.004015 0.004892 01020 # 01021 # 01022 # On Jupiter: 01023 # 01024 # n Std Dev RMS Mean Variance Sum 01025 # 5304 0.08512 0.08629 0.01418 0.007245 75.21 01026 # 01027 # Flux Med |Dev| IntQtlRng Median Min Max 01028 # 4.695 0.0008142 0.001657 0.0001557 -0.004526 1.076 01029 # 01030 # Estimated dynamic range = 1.076 / 0.001389 = 775 (better) 01031 # 01032 # Note that the exact numbers you get will depend on how deep you 01033 # take the interactive clean and how you draw the box for the stats. 01034 # 01035 #--------------------------------------------------------------------- 01036 # 01037 # Next self-cal cycle 01038 # 01039 print '--SelfCal 2--' 01040 default('gaincal') 01041 01042 vis = srcsplitms 01043 01044 selfcaltab2 = srcsplitms + '.selfcal2' 01045 caltable = selfcaltab2 01046 01047 selectdata = False 01048 gaincurve = False 01049 opacity = 0.0 01050 refant = '11' 01051 01052 # This time amp+phase on 10s timescales SNR>1 01053 gaintype = 'G' 01054 calmode = 'ap' 01055 solint = 10.0 01056 minsnr = 1.0 01057 solnorm = False 01058 01059 gaincal() 01060 01061 # 01062 # It is useful to put this up in plotcal 01063 # 01064 #--------------------------------------------------------------------- 01065 # 01066 print '--PlotCal--' 01067 default('plotcal') 01068 01069 tablein = selfcaltab2 01070 multiplot = True 01071 yaxis = 'amp' 01072 01073 plotcal() 01074 01075 # Use the Next button to iterate over antennas 01076 01077 # Pause script if you are running in scriptmode 01078 if scriptmode: 01079 user_check=raw_input('Return to continue script\n') 01080 01081 yaxis = 'phase' 01082 01083 plotcal() 01084 01085 # 01086 # You can see it is not too noisy. 01087 # 01088 # Pause script if you are running in scriptmode 01089 if scriptmode: 01090 user_check=raw_input('Return to continue script\n') 01091 01092 # Lets do some smoothing anyway. 01093 # 01094 #--------------------------------------------------------------------- 01095 # 01096 # Smooth calibration solutions 01097 # 01098 print '--Smooth--' 01099 default('smoothcal') 01100 01101 vis = srcsplitms 01102 01103 tablein = selfcaltab2 01104 01105 smoothcaltab2 = srcsplitms + '.smoothcal2' 01106 caltable = smoothcaltab2 01107 01108 # Do a 30s boxcar average 01109 smoothtype = 'mean' 01110 smoothtime = 30.0 01111 01112 smoothcal() 01113 01114 # If you put into plotcal you'll see the results 01115 # For example, you can grap the inputs from the last 01116 # time you ran plotcal, set the new tablename, and plot! 01117 #run plotcal.last 01118 #tablein = smoothcaltab2 01119 #plotcal() 01120 01121 # 01122 #--------------------------------------------------------------------- 01123 # 01124 # Correct the data 01125 # 01126 print '--ApplyCal--' 01127 default('applycal') 01128 01129 vis = srcsplitms 01130 01131 gaintable = smoothcaltab2 01132 01133 gaincurve = False 01134 opacity = 0.0 01135 field = '' 01136 spw = '' 01137 selectdata = False 01138 calwt = True 01139 01140 applycal() 01141 01142 # 01143 #===================================================================== 01144 # THIRD CLEAN / SELFCAL CYCLE 01145 #===================================================================== 01146 # 01147 print '--Clean 3--' 01148 default('clean') 01149 01150 vis = srcsplitms 01151 01152 imname3 = prefix + '.clean3' 01153 imagename = imname3 01154 01155 field = '*' 01156 spw = '' 01157 mode = 'mfs' 01158 gain = 0.1 01159 niter = 10000 01160 threshold=0.04 01161 01162 alg = clnalg 01163 imsize = clnimsize 01164 cell = clncell 01165 01166 weighting = 'briggs' 01167 rmode = 'norm' 01168 robust = 0.5 01169 01170 cleanbox = 'interactive' 01171 npercycle = 100 01172 01173 clean() 01174 01175 # Cleans alot deeper 01176 # You can change the npercycle to larger numbers 01177 # (like 250 or so) as you get deeper also. 01178 01179 # Set up variables 01180 clnimage3 = imname3+'.image' 01181 clnmodel3 = imname3+'.model' 01182 clnresid3 = imname3+'.residual' 01183 clnmask3 = imname3+'.clean_interactive.mask' 01184 01185 # 01186 #--------------------------------------------------------------------- 01187 # 01188 # Look at this in viewer 01189 viewer(clnimage3,'image') 01190 01191 # jupiter6cm.usecase.clean3.image (Jy/beam) 01192 # 01193 # n Std Dev RMS Mean Variance Sum 01194 # 5848 0.001015 0.001015 -4.036e-06 1.029e-06 -0.02360 01195 # 01196 # Flux Med |Dev| IntQtlRng Median Min Max 01197 # -0.001470 0.0006728 0.001347 8.245e-06 -0.003260 0.003542 01198 # 01199 # 01200 # On Jupiter: 01201 # 01202 # n Std Dev RMS Mean Variance Sum 01203 # 6003 0.08012 0.08107 0.01245 0.006419 74.72 01204 # 01205 # Flux Med |Dev| IntQtlRng Median Min Max 01206 # 4.653 0.0006676 0.001383 -1.892e-06 -0.002842 1.076 01207 # 01208 # Estimated dynamic range = 1.076 / 0.001015 = 1060 (even better!) 01209 # 01210 # Note that the exact numbers you get will depend on how deep you 01211 # take the interactive clean and how you draw the box for the stats. 01212 # 01213 # Greg Taylor got 1600:1 so we still have some ways to go 01214 # This will probably take several more careful self-cal cycles. 01215 01216 # Set up final variables 01217 clnimage = clnimage3 01218 clnmodel = clnmodel3 01219 clnresid = clnresid3 01220 clnmask = clnmask3 01221 01222 #===================================================================== 01223 # 01224 # Export the Final CLEAN Image as FITS 01225 # 01226 print '--Final Export CLEAN FITS--' 01227 default('exportfits') 01228 01229 clnfits = prefix + '.clean.fits' 01230 01231 imagename = clnimage 01232 fitsimage = clnfits 01233 01234 # Run asynchronously so as not to interfere with other tasks 01235 # (BETA: also avoids crash on next importfits) 01236 async = True 01237 01238 exportfits() 01239 01240 #===================================================================== 01241 # 01242 # Export the Final Self-Calibrated Jupiter data as UVFITS 01243 # 01244 print '--Final Export UVFITS--' 01245 default('exportuvfits') 01246 01247 caluvfits = prefix + '.selfcal.uvfits' 01248 01249 vis = srcsplitms 01250 fitsfile = caluvfits 01251 01252 # The self-calibrated data is in the CORRECTED_DATA column 01253 datacolumn = 'corrected' 01254 01255 # Write as a multisource UVFITS (with SU table) 01256 # even though it will have only one field in it 01257 multisource = True 01258 01259 # Run asynchronously so as not to interfere with other tasks 01260 # (BETA: also avoids crash on next importuvfits) 01261 async = True 01262 01263 exportuvfits() 01264 01265 # 01266 #===================================================================== 01267 # Image Analysis 01268 #===================================================================== 01269 # 01270 # Can do some image statistics if you wish 01271 # Treat this like a regression script 01272 # WARNING: currently requires toolkit 01273 # 01274 print ' Jupiter results ' 01275 print ' =============== ' 01276 01277 print '' 01278 # Pull the max src amp value out of the MS 01279 ms.open(srcsplitms) 01280 thistest_src = max(ms.range(["amplitude"]).get('amplitude')) 01281 oldtest_src = 4.92000198364 01282 print ' MS max amplitude should be ',oldtest_src 01283 print ' Found : Max in MS = ',thistest_src 01284 diff_src = abs((oldtest_src-thistest_src)/oldtest_src) 01285 print ' Difference (fractional) = ',diff_src 01286 01287 ms.close() 01288 01289 print '' 01290 # Pull the max and rms from the clean image 01291 ia.open(clnimage) 01292 on_statistics=ia.statistics() 01293 thistest_immax=on_statistics['max'][0] 01294 oldtest_immax = 1.07732224464 01295 print ' Clean image ON-SRC max should be ',oldtest_immax 01296 print ' Found : Max in image = ',thistest_immax 01297 diff_immax = abs((oldtest_immax-thistest_immax)/oldtest_immax) 01298 print ' Difference (fractional) = ',diff_immax 01299 01300 print '' 01301 # Now do stats in the lower right corner of the image 01302 box = ia.setboxregion([0.75,0.00],[1.00,0.25],frac=true) 01303 off_statistics=ia.statistics(region=box) 01304 thistest_imrms=off_statistics['rms'][0] 01305 oldtest_imrms = 0.0010449 01306 print ' Clean image OFF-SRC rms should be ',oldtest_imrms 01307 print ' Found : rms in image = ',thistest_imrms 01308 diff_imrms = abs((oldtest_imrms-thistest_imrms)/oldtest_imrms) 01309 print ' Difference (fractional) = ',diff_imrms 01310 01311 print '' 01312 print ' Final Clean image Dynamic Range = ',thistest_immax/thistest_imrms 01313 print '' 01314 print ' =============== ' 01315 01316 ia.close() 01317 01318 print '' 01319 print '--- Done ---' 01320 01321 # 01322 #=====================================================================
1.8.0
