|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ###################################################################### 00002 # # 00003 # Use Case Script for Jupiter 6cm VLA # 00004 # Trimmed down from Use Case jupiter6cm_usecase.py # 00005 # # 00006 # Updated STM 2008-05-15 (Beta Patch 2.0) # 00007 # Updated STM 2008-06-11 (Beta Patch 2.0) # 00008 # Updated STM 2008-06-12 (Beta Patch 2.0) for summer school demo # 00009 # # 00010 # This is a VLA 6cm dataset that was observed in 1999 to set the # 00011 # flux scale for calibration of the VLA. Included in the program # 00012 # were observations of the planets, including Jupiter. # 00013 # # 00014 # This is D-configuration data, with resolution of around 14" # 00015 # # 00016 # Includes polarization imaging and analysis # 00017 # # 00018 ###################################################################### 00019 00020 import time 00021 import os 00022 00023 # 00024 #===================================================================== 00025 # 00026 # This script has some interactive commands: scriptmode = True 00027 # if you are running it and want it to stop during interactive parts. 00028 00029 scriptmode = True 00030 00031 #===================================================================== 00032 # 00033 # Set up some useful variables - these will be set during the script 00034 # also, but if you want to restart the script in the middle here 00035 # they are in one place: 00036 00037 # This will prefix all output file names 00038 prefix='jupiter6cm.demo' 00039 00040 # Clean up old files 00041 os.system('rm -rf '+prefix+'*') 00042 00043 # This is the output MS file name 00044 msfile = prefix + '.ms' 00045 00046 # 00047 #===================================================================== 00048 # Calibration variables 00049 # 00050 # Use same prefix as rest of script 00051 calprefix = prefix 00052 00053 # spectral windows to process 00054 usespw = '' 00055 usespwlist = ['0','1'] 00056 00057 # prior calibration to apply 00058 usegaincurve = True 00059 gainopacity = 0.0 00060 00061 # reference antenna 11 (11=VLA:N1) 00062 calrefant = '11' 00063 00064 gtable = calprefix + '.gcal' 00065 ftable = calprefix + '.fluxscale' 00066 atable = calprefix + '.accum' 00067 00068 # 00069 #===================================================================== 00070 # Polarization calibration setup 00071 # 00072 dopolcal = True 00073 00074 ptable = calprefix + '.pcal' 00075 xtable = calprefix + '.polx' 00076 00077 # Pol leakage calibrator 00078 poldfield = '0137+331' 00079 00080 # Pol angle calibrator 00081 polxfield = '1331+305' 00082 # At Cband the fractional polarization of this source is 0.112 and 00083 # the R-L PhaseDiff = 66deg (EVPA = 33deg) 00084 polxfpol = 0.112 00085 polxrlpd_deg = 66.0 00086 # Dictionary of IPOL in the spw 00087 polxipol = {'0' : 7.462, 00088 '1' : 7.510} 00089 00090 # Make Stokes lists for setjy 00091 polxiquv = {} 00092 for spw in ['0','1']: 00093 ipol = polxipol[spw] 00094 fpol = polxfpol 00095 ppol = ipol*fpol 00096 rlpd = polxrlpd_deg*pi/180.0 00097 qpol = ppol*cos(rlpd) 00098 upol = ppol*sin(rlpd) 00099 polxiquv[spw] = [ipol,qpol,upol,0.0] 00100 00101 # 00102 # Split output setup 00103 # 00104 srcname = 'JUPITER' 00105 srcsplitms = calprefix + '.' + srcname + '.split.ms' 00106 calname = '0137+331' 00107 calsplitms = calprefix + '.' + calname + '.split.ms' 00108 00109 # 00110 #===================================================================== 00111 # 00112 # Intensity imaging parameters 00113 # 00114 # Same prefix for this imaging demo output 00115 # 00116 imprefix = prefix 00117 00118 # This is D-config VLA 6cm (4.85GHz) obs 00119 # Check the observational status summary 00120 # Primary beam FWHM = 45'/f_GHz = 557" 00121 # Synthesized beam FWHM = 14" 00122 # RMS in 10min (600s) = 0.06 mJy (thats now, but close enough) 00123 00124 # Set the output image size and cell size (arcsec) 00125 # 4" will give 3.5x oversampling 00126 clncell = [4.,4.] 00127 00128 # 280 pix will cover to 2xPrimaryBeam 00129 # clean will say to use 288 (a composite integer) for efficiency 00130 clnalg = 'clark' 00131 clnmode = '' 00132 # For Cotton-Schwab use 00133 clnmode = 'csclean' 00134 clnimsize = [288,288] 00135 00136 # iterations 00137 clniter = 10000 00138 00139 # Also set flux residual threshold (0.04 mJy) 00140 # From our listobs: 00141 # Total integration time = 85133.2 seconds 00142 # With rms of 0.06 mJy in 600s ==> rms = 0.005 mJy 00143 # Set to 10x thermal rms 00144 clnthreshold=0.05 00145 00146 # 00147 # Filenames 00148 # 00149 imname1 = imprefix + '.clean1' 00150 clnimage1 = imname1+'.image' 00151 clnmodel1 = imname1+'.model' 00152 clnresid1 = imname1+'.residual' 00153 clnmask1 = imname1+'.clean_interactive.mask' 00154 00155 imname2 = imprefix + '.clean2' 00156 clnimage2 = imname2+'.image' 00157 clnmodel2 = imname2+'.model' 00158 clnresid2 = imname2+'.residual' 00159 clnmask2 = imname2+'.clean_interactive.mask' 00160 00161 imname3 = imprefix + '.clean3' 00162 clnimage3 = imname3+'.image' 00163 clnmodel3 = imname3+'.model' 00164 clnresid3 = imname3+'.residual' 00165 clnmask3 = imname3+'.clean_interactive.mask' 00166 00167 # 00168 # Selfcal parameters 00169 # 00170 # reference antenna 11 (11=VLA:N1) 00171 calrefant = '11' 00172 00173 # 00174 # Filenames 00175 # 00176 selfcaltab1 = imprefix + '.selfcal1.gtable' 00177 00178 selfcaltab2 = imprefix + '.selfcal2.gtable' 00179 smoothcaltab2 = imprefix + '.smoothcal2.gtable' 00180 00181 # 00182 #===================================================================== 00183 # 00184 # Polarization imaging parameters 00185 # 00186 # New prefix for polarization imaging output 00187 # 00188 polprefix = prefix + '.polimg' 00189 00190 # Set up clean slightly differently 00191 polclnalg = 'hogbom' 00192 polclnmode = 'csclean' 00193 00194 polimname = polprefix + '.clean' 00195 polimage = polimname+'.image' 00196 polmodel = polimname+'.model' 00197 polresid = polimname+'.residual' 00198 polmask = polimname+'.clean_interactive.mask' 00199 00200 # 00201 # Other files 00202 # 00203 ipolimage = polimage+'.I' 00204 qpolimage = polimage+'.Q' 00205 upolimage = polimage+'.U' 00206 00207 poliimage = polimage+'.poli' 00208 polaimage = polimage+'.pola' 00209 00210 # 00211 #===================================================================== 00212 # Start processing 00213 #===================================================================== 00214 # 00215 # Get to path to the CASA home and stip off the name 00216 pathname=os.environ.get('CASAPATH').split()[0] 00217 00218 # This is where the UVFITS data should be 00219 #fitsdata=pathname+'/data/demo/jupiter6cm.fits' 00220 # Or 00221 #fitsdata=pathname+'/data/nrao/VLA/planets_6cm.fits' 00222 #fitsdata='/home/ballista/casa/devel/data/nrao/VLA/planets_6cm.fits' 00223 # 00224 # Can also be found online at 00225 #http://casa.nrao.edu/Data/VLA/Planets6cm/planets_6cm.fits 00226 00227 # Use version in current directory 00228 fitsdata='planets_6cm.fits' 00229 00230 # 00231 #===================================================================== 00232 # Data Import and List 00233 #===================================================================== 00234 # 00235 # Import the data from FITS to MS 00236 # 00237 print '--Import--' 00238 00239 # Safest to start from task defaults 00240 default('importuvfits') 00241 00242 print "Use importuvfits to read UVFITS and make an MS" 00243 00244 # Set up the MS filename and save as new global variable 00245 msfile = prefix + '.ms' 00246 00247 print "MS will be called "+msfile 00248 00249 # Use task importuvfits 00250 fitsfile = fitsdata 00251 vis = msfile 00252 importuvfits() 00253 00254 #===================================================================== 00255 # 00256 # List a summary of the MS 00257 # 00258 print '--Listobs--' 00259 00260 # Don't default this one and make use of the previous setting of 00261 # vis. Remember, the variables are GLOBAL! 00262 00263 print "Use listobs to print verbose summary to logger" 00264 00265 # You may wish to see more detailed information, in this case 00266 # use the verbose = True option 00267 verbose = True 00268 00269 listobs() 00270 00271 # You should get in your logger window and in the casapy.log file 00272 # something like: 00273 # 00274 # Observer: FLUX99 Project: 00275 # Observation: VLA 00276 # 00277 # Data records: 2021424 Total integration time = 85133.2 seconds 00278 # Observed from 23:15:27 to 22:54:20 00279 # 00280 # ObservationID = 0 ArrayID = 0 00281 # Date Timerange Scan FldId FieldName SpwIds 00282 # 15-Apr-1999/23:15:26.7 - 23:16:10.0 1 0 0137+331 [0, 1] 00283 # 23:38:40.0 - 23:48:00.0 2 1 0813+482 [0, 1] 00284 # 23:53:40.0 - 23:55:20.0 3 2 0542+498 [0, 1] 00285 # 16-Apr-1999/00:22:10.1 - 00:23:49.9 4 3 0437+296 [0, 1] 00286 # 00:28:23.3 - 00:30:00.1 5 4 VENUS [0, 1] 00287 # 00:48:40.0 - 00:50:20.0 6 1 0813+482 [0, 1] 00288 # 00:56:13.4 - 00:57:49.9 7 2 0542+498 [0, 1] 00289 # 01:10:20.1 - 01:11:59.9 8 5 0521+166 [0, 1] 00290 # 01:23:29.9 - 01:25:00.1 9 3 0437+296 [0, 1] 00291 # 01:29:33.3 - 01:31:10.0 10 4 VENUS [0, 1] 00292 # 01:49:50.0 - 01:51:30.0 11 6 1411+522 [0, 1] 00293 # 02:03:00.0 - 02:04:30.0 12 7 1331+305 [0, 1] 00294 # 02:17:30.0 - 02:19:10.0 13 1 0813+482 [0, 1] 00295 # 02:24:20.0 - 02:26:00.0 14 2 0542+498 [0, 1] 00296 # 02:37:49.9 - 02:39:30.0 15 5 0521+166 [0, 1] 00297 # 02:50:50.1 - 02:52:20.1 16 3 0437+296 [0, 1] 00298 # 02:59:20.0 - 03:01:00.0 17 6 1411+522 [0, 1] 00299 # 03:12:30.0 - 03:14:10.0 18 7 1331+305 [0, 1] 00300 # 03:27:53.3 - 03:29:39.9 19 1 0813+482 [0, 1] 00301 # 03:35:00.0 - 03:36:40.0 20 2 0542+498 [0, 1] 00302 # 03:49:50.0 - 03:51:30.1 21 6 1411+522 [0, 1] 00303 # 04:03:10.0 - 04:04:50.0 22 7 1331+305 [0, 1] 00304 # 04:18:49.9 - 04:20:40.0 23 1 0813+482 [0, 1] 00305 # 04:25:56.6 - 04:27:39.9 24 2 0542+498 [0, 1] 00306 # 04:42:49.9 - 04:44:40.0 25 8 MARS [0, 1] 00307 # 04:56:50.0 - 04:58:30.1 26 6 1411+522 [0, 1] 00308 # 05:24:03.3 - 05:33:39.9 27 7 1331+305 [0, 1] 00309 # 05:48:00.0 - 05:49:49.9 28 1 0813+482 [0, 1] 00310 # 05:58:36.6 - 06:00:30.0 29 8 MARS [0, 1] 00311 # 06:13:20.1 - 06:14:59.9 30 6 1411+522 [0, 1] 00312 # 06:27:40.0 - 06:29:20.0 31 7 1331+305 [0, 1] 00313 # 06:44:13.4 - 06:46:00.0 32 1 0813+482 [0, 1] 00314 # 06:55:06.6 - 06:57:00.0 33 8 MARS [0, 1] 00315 # 07:10:40.0 - 07:12:20.0 34 6 1411+522 [0, 1] 00316 # 07:28:20.0 - 07:30:10.1 35 7 1331+305 [0, 1] 00317 # 07:42:49.9 - 07:44:30.0 36 8 MARS [0, 1] 00318 # 07:58:43.3 - 08:00:39.9 37 6 1411+522 [0, 1] 00319 # 08:13:30.0 - 08:15:19.9 38 7 1331+305 [0, 1] 00320 # 08:27:53.4 - 08:29:30.0 39 8 MARS [0, 1] 00321 # 08:42:59.9 - 08:44:50.0 40 6 1411+522 [0, 1] 00322 # 08:57:09.9 - 08:58:50.0 41 7 1331+305 [0, 1] 00323 # 09:13:03.3 - 09:14:50.1 42 9 NGC7027 [0, 1] 00324 # 09:26:59.9 - 09:28:40.0 43 6 1411+522 [0, 1] 00325 # 09:40:33.4 - 09:42:09.9 44 7 1331+305 [0, 1] 00326 # 09:56:19.9 - 09:58:10.0 45 9 NGC7027 [0, 1] 00327 # 10:12:59.9 - 10:14:50.0 46 8 MARS [0, 1] 00328 # 10:27:09.9 - 10:28:50.0 47 6 1411+522 [0, 1] 00329 # 10:40:30.0 - 10:42:00.0 48 7 1331+305 [0, 1] 00330 # 10:56:10.0 - 10:57:50.0 49 9 NGC7027 [0, 1] 00331 # 11:28:30.0 - 11:35:30.0 50 10 NEPTUNE [0, 1] 00332 # 11:48:20.0 - 11:50:10.0 51 6 1411+522 [0, 1] 00333 # 12:01:36.7 - 12:03:10.0 52 7 1331+305 [0, 1] 00334 # 12:35:33.3 - 12:37:40.0 53 11 URANUS [0, 1] 00335 # 12:46:30.0 - 12:48:10.0 54 10 NEPTUNE [0, 1] 00336 # 13:00:29.9 - 13:02:10.0 55 6 1411+522 [0, 1] 00337 # 13:15:23.3 - 13:17:10.1 56 9 NGC7027 [0, 1] 00338 # 13:33:43.3 - 13:35:40.0 57 11 URANUS [0, 1] 00339 # 13:44:30.0 - 13:46:10.0 58 10 NEPTUNE [0, 1] 00340 # 14:00:46.7 - 14:01:39.9 59 0 0137+331 [0, 1] 00341 # 14:10:40.0 - 14:12:09.9 60 12 JUPITER [0, 1] 00342 # 14:24:06.6 - 14:25:40.1 61 11 URANUS [0, 1] 00343 # 14:34:30.0 - 14:36:10.1 62 10 NEPTUNE [0, 1] 00344 # 14:59:13.4 - 15:00:00.0 63 0 0137+331 [0, 1] 00345 # 15:09:03.3 - 15:10:40.1 64 12 JUPITER [0, 1] 00346 # 15:24:30.0 - 15:26:20.1 65 9 NGC7027 [0, 1] 00347 # 15:40:10.0 - 15:45:00.0 66 11 URANUS [0, 1] 00348 # 15:53:50.0 - 15:55:20.0 67 10 NEPTUNE [0, 1] 00349 # 16:18:53.4 - 16:19:49.9 68 0 0137+331 [0, 1] 00350 # 16:29:10.1 - 16:30:49.9 69 12 JUPITER [0, 1] 00351 # 16:42:53.4 - 16:44:30.0 70 11 URANUS [0, 1] 00352 # 16:54:53.4 - 16:56:40.0 71 9 NGC7027 [0, 1] 00353 # 17:23:06.6 - 17:30:40.0 72 2 0542+498 [0, 1] 00354 # 17:41:50.0 - 17:43:20.0 73 3 0437+296 [0, 1] 00355 # 17:55:36.7 - 17:57:39.9 74 4 VENUS [0, 1] 00356 # 18:19:23.3 - 18:20:09.9 75 0 0137+331 [0, 1] 00357 # 18:30:23.3 - 18:32:00.0 76 12 JUPITER [0, 1] 00358 # 18:44:49.9 - 18:46:30.0 77 9 NGC7027 [0, 1] 00359 # 18:59:13.3 - 19:00:59.9 78 2 0542+498 [0, 1] 00360 # 19:19:10.0 - 19:21:20.1 79 5 0521+166 [0, 1] 00361 # 19:32:50.1 - 19:34:29.9 80 3 0437+296 [0, 1] 00362 # 19:39:03.3 - 19:40:40.1 81 4 VENUS [0, 1] 00363 # 20:08:06.7 - 20:08:59.9 82 0 0137+331 [0, 1] 00364 # 20:18:10.0 - 20:19:50.0 83 12 JUPITER [0, 1] 00365 # 20:33:53.3 - 20:35:40.1 84 1 0813+482 [0, 1] 00366 # 20:40:59.9 - 20:42:40.0 85 2 0542+498 [0, 1] 00367 # 21:00:16.6 - 21:02:20.1 86 5 0521+166 [0, 1] 00368 # 21:13:53.4 - 21:15:29.9 87 3 0437+296 [0, 1] 00369 # 21:20:43.4 - 21:22:30.0 88 4 VENUS [0, 1] 00370 # 21:47:26.7 - 21:48:20.1 89 0 0137+331 [0, 1] 00371 # 21:57:30.0 - 21:59:10.0 90 12 JUPITER [0, 1] 00372 # 22:12:13.3 - 22:14:00.1 91 2 0542+498 [0, 1] 00373 # 22:28:33.3 - 22:30:19.9 92 4 VENUS [0, 1] 00374 # 22:53:33.3 - 22:54:19.9 93 0 0137+331 [0, 1] 00375 # 00376 # Fields: 13 00377 # ID Name Right Ascension Declination Epoch 00378 # 0 0137+331 01:37:41.30 +33.09.35.13 J2000 00379 # 1 0813+482 08:13:36.05 +48.13.02.26 J2000 00380 # 2 0542+498 05:42:36.14 +49.51.07.23 J2000 00381 # 3 0437+296 04:37:04.17 +29.40.15.14 J2000 00382 # 4 VENUS 04:06:54.11 +22.30.35.91 J2000 00383 # 5 0521+166 05:21:09.89 +16.38.22.05 J2000 00384 # 6 1411+522 14:11:20.65 +52.12.09.14 J2000 00385 # 7 1331+305 13:31:08.29 +30.30.32.96 J2000 00386 # 8 MARS 14:21:41.37 -12.21.49.45 J2000 00387 # 9 NGC7027 21:07:01.59 +42.14.10.19 J2000 00388 # 10 NEPTUNE 20:26:01.14 -18.54.54.21 J2000 00389 # 11 URANUS 21:15:42.83 -16.35.05.59 J2000 00390 # 12 JUPITER 00:55:34.04 +04.45.44.71 J2000 00391 # 00392 # Spectral Windows: (2 unique spectral windows and 1 unique polarization setups) 00393 # SpwID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs 00394 # 0 1 TOPO 4885.1 50000 50000 4885.1 RR RL LR LL 00395 # 1 1 TOPO 4835.1 50000 50000 4835.1 RR RL LR LL 00396 # 00397 # Feeds: 28: printing first row only 00398 # Antenna Spectral Window # Receptors Polarizations 00399 # 1 -1 2 [ R, L] 00400 # 00401 # Antennas: 27: 00402 # ID Name Station Diam. Long. Lat. 00403 # 0 1 VLA:W9 25.0 m -107.37.25.1 +33.53.51.0 00404 # 1 2 VLA:N9 25.0 m -107.37.07.8 +33.54.19.0 00405 # 2 3 VLA:N3 25.0 m -107.37.06.3 +33.54.04.8 00406 # 3 4 VLA:N5 25.0 m -107.37.06.7 +33.54.08.0 00407 # 4 5 VLA:N2 25.0 m -107.37.06.2 +33.54.03.5 00408 # 5 6 VLA:E1 25.0 m -107.37.05.7 +33.53.59.2 00409 # 6 7 VLA:E2 25.0 m -107.37.04.4 +33.54.01.1 00410 # 7 8 VLA:N8 25.0 m -107.37.07.5 +33.54.15.8 00411 # 8 9 VLA:E8 25.0 m -107.36.48.9 +33.53.55.1 00412 # 9 10 VLA:W3 25.0 m -107.37.08.9 +33.54.00.1 00413 # 10 11 VLA:N1 25.0 m -107.37.06.0 +33.54.01.8 00414 # 11 12 VLA:E6 25.0 m -107.36.55.6 +33.53.57.7 00415 # 12 13 VLA:W7 25.0 m -107.37.18.4 +33.53.54.8 00416 # 13 14 VLA:E4 25.0 m -107.37.00.8 +33.53.59.7 00417 # 14 15 VLA:N7 25.0 m -107.37.07.2 +33.54.12.9 00418 # 15 16 VLA:W4 25.0 m -107.37.10.8 +33.53.59.1 00419 # 16 17 VLA:W5 25.0 m -107.37.13.0 +33.53.57.8 00420 # 17 18 VLA:N6 25.0 m -107.37.06.9 +33.54.10.3 00421 # 18 19 VLA:E7 25.0 m -107.36.52.4 +33.53.56.5 00422 # 19 20 VLA:E9 25.0 m -107.36.45.1 +33.53.53.6 00423 # 21 22 VLA:W8 25.0 m -107.37.21.6 +33.53.53.0 00424 # 22 23 VLA:W6 25.0 m -107.37.15.6 +33.53.56.4 00425 # 23 24 VLA:W1 25.0 m -107.37.05.9 +33.54.00.5 00426 # 24 25 VLA:W2 25.0 m -107.37.07.4 +33.54.00.9 00427 # 25 26 VLA:E5 25.0 m -107.36.58.4 +33.53.58.8 00428 # 26 27 VLA:N4 25.0 m -107.37.06.5 +33.54.06.1 00429 # 27 28 VLA:E3 25.0 m -107.37.02.8 +33.54.00.5 00430 # 00431 # Tables: 00432 # MAIN 2021424 rows 00433 # ANTENNA 28 rows 00434 # DATA_DESCRIPTION 2 rows 00435 # DOPPLER <absent> 00436 # FEED 28 rows 00437 # FIELD 13 rows 00438 # FLAG_CMD <empty> 00439 # FREQ_OFFSET <absent> 00440 # HISTORY 7058 rows 00441 # OBSERVATION 1 row 00442 # POINTING 2604 rows 00443 # POLARIZATION 1 row 00444 # PROCESSOR <empty> 00445 # SOURCE <empty> (see FIELD) 00446 # SPECTRAL_WINDOW 2 rows 00447 # STATE <empty> 00448 # SYSCAL <absent> 00449 # WEATHER <absent> 00450 00451 # 00452 #===================================================================== 00453 # Data Examination and Flagging 00454 #===================================================================== 00455 # 00456 # Use Plotxy to interactively flag the data 00457 # 00458 print '--Plotxy--' 00459 default('plotxy') 00460 00461 print "Now we use plotxy to examine and interactively flag data" 00462 00463 vis = msfile 00464 00465 # The fields we are interested in: 1331+305,JUPITER,0137+331 00466 selectdata = True 00467 00468 # First we do the primary calibrator 00469 field = '1331+305' 00470 00471 # Plot only the RR and LL for now 00472 correlation = 'RR LL' 00473 00474 # Plot amplitude vs. uvdist 00475 xaxis = 'uvdist' 00476 yaxis = 'amp' 00477 multicolor = 'both' 00478 00479 # Use the field name as the title 00480 selectplot = True 00481 title = field+" " 00482 00483 iteration = '' 00484 00485 plotxy() 00486 00487 print "" 00488 print "-----------------------------------------------------" 00489 print "Plotxy" 00490 print "Showing 1331+305 RR LL for all antennas" 00491 print "Use MarkRegion then draw boxes around points to flag" 00492 print "You can use ESC to drop last drawn box" 00493 print "When happy with boxes, hit Flag to flag" 00494 print "You can repeat as necessary" 00495 00496 # Pause script if you are running in scriptmode 00497 if scriptmode: 00498 user_check=raw_input('Return to continue script\n') 00499 00500 # You can also use flagdata to do this non-interactively 00501 # (see below) 00502 00503 # Now look at the cross-polar products 00504 correlation = 'RL LR' 00505 00506 plotxy() 00507 00508 print "" 00509 print "-----------------------------------------------------" 00510 print "Looking at RL LR" 00511 print "Now flag the bad data here" 00512 00513 # Pause script if you are running in scriptmode 00514 if scriptmode: 00515 user_check=raw_input('Return to continue script\n') 00516 00517 #--------------------------------------------------------------------- 00518 # Now do calibrater 0137+331 00519 field = '0137+331' 00520 correlation = 'RR LL' 00521 xaxis = 'uvdist' 00522 spw = '' 00523 iteration = '' 00524 antenna = '' 00525 00526 title = field+" " 00527 00528 plotxy() 00529 00530 # You'll see a bunch of bad data along the bottom near zero amp 00531 # Draw a box around some of it and use Locate 00532 # Looks like much of it is Antenna 9 (ID=8) in spw=1 00533 00534 print "" 00535 print "-----------------------------------------------------" 00536 print "Plotting 0137+331 RR LL all antennas" 00537 print "You see bad data along bottom" 00538 print "Mark a box around a bit of it and hit Locate" 00539 print "Look in logger to see what it is" 00540 print "You see much is Antenna 9 (ID=8) in spw 1" 00541 00542 # Pause script if you are running in scriptmode 00543 if scriptmode: 00544 user_check=raw_input('Return to continue script\n') 00545 00546 xaxis = 'time' 00547 spw = '1' 00548 correlation = '' 00549 00550 # Note that the strings like antenna='9' first try to match the 00551 # NAME which we see in listobs was the number '9' for ID=8. 00552 # So be careful here (why naming antennas as numbers is bad). 00553 antenna = '9' 00554 00555 plotxy() 00556 00557 # YES! the last 4 scans are bad. Box 'em and flag. 00558 00559 print "" 00560 print "-----------------------------------------------------" 00561 print "Plotting vs. time antenna='9' and spw='1' " 00562 print "Box up last 4 scans which are bad and Flag" 00563 00564 # Pause script if you are running in scriptmode 00565 if scriptmode: 00566 user_check=raw_input('Return to continue script\n') 00567 00568 # Go back and clean up 00569 xaxis = 'uvdist' 00570 spw = '' 00571 antenna = '' 00572 correlation = 'RR LL' 00573 00574 plotxy() 00575 00576 # Box up the bad low points (basically a clip below 0.52) and flag 00577 00578 # Note that RL,LR are too weak to clip on. 00579 00580 print "" 00581 print "-----------------------------------------------------" 00582 print "Back to all data" 00583 print "Clean up remaining bad points" 00584 00585 # Pause script if you are running in scriptmode 00586 if scriptmode: 00587 user_check=raw_input('Return to continue script\n') 00588 00589 #--------------------------------------------------------------------- 00590 # Finally, do JUPITER 00591 field = 'JUPITER' 00592 correlation = 'RR LL' 00593 iteration = '' 00594 xaxis = 'uvdist' 00595 00596 title = field+" " 00597 00598 plotxy() 00599 00600 # Here you will see that the final scan at 22:00:00 UT is bad 00601 # Draw a box around it and flag it! 00602 00603 print "" 00604 print "-----------------------------------------------------" 00605 print "Now plot JUPITER versus uvdist" 00606 print "Lots of bad stuff near bottom" 00607 print "Lets go and find it - try Locate" 00608 print "Looks like lots of different antennas but at same time" 00609 00610 # Pause script if you are running in scriptmode 00611 if scriptmode: 00612 user_check=raw_input('Return to continue script\n') 00613 00614 correlation = '' 00615 xaxis = 'time' 00616 00617 plotxy() 00618 00619 # Here you will see that the final scan at 22:00:00 UT is bad 00620 # Draw a box around it and flag it! 00621 00622 print "" 00623 print "-----------------------------------------------------" 00624 print "Now plotting vs. time" 00625 print "See bad scan at end - flag it!" 00626 00627 # Pause script if you are running in scriptmode 00628 if scriptmode: 00629 user_check=raw_input('Return to continue script\n') 00630 00631 # Now look at whats left 00632 correlation = 'RR LL' 00633 xaxis = 'uvdist' 00634 spw = '1' 00635 antenna = '' 00636 iteration = 'antenna' 00637 00638 plotxy() 00639 00640 # As you step through, you will see that Antenna 9 (ID=8) is often 00641 # bad in this spw. If you box and do Locate (or remember from 00642 # 0137+331) its probably a bad time. 00643 00644 print "" 00645 print "-----------------------------------------------------" 00646 print "Looking now at SPW 1" 00647 print "Now we set iteration to Antenna" 00648 print "Step through antennas with Next" 00649 print "See bad Antenna 9 (ID 8) as in 0137+331" 00650 00651 # Pause script if you are running in scriptmode 00652 if scriptmode: 00653 user_check=raw_input('Return to continue script\n') 00654 00655 # The easiset way to kill it: 00656 00657 antenna = '9' 00658 iteration = '' 00659 xaxis = 'time' 00660 correlation = '' 00661 00662 plotxy() 00663 00664 # Draw a box around all points in the last bad scans and flag 'em! 00665 00666 print "" 00667 print "-----------------------------------------------------" 00668 print "Now plotting vs. time antenna 9 spw 1" 00669 print "Box up the bad scans and Flag" 00670 00671 # Pause script if you are running in scriptmode 00672 if scriptmode: 00673 user_check=raw_input('Return to continue script\n') 00674 00675 # Now clean up the rest 00676 xaxis = 'uvdist' 00677 correlation = 'RR LL' 00678 antenna = '' 00679 spw = '' 00680 00681 # You will be drawing many tiny boxes, so remember you can 00682 # use the ESC key to get rid of the most recent box if you 00683 # make a mistake. 00684 00685 plotxy() 00686 00687 # Note that the end result is we've flagged lots of points 00688 # in RR and LL. We will rely upon imager to ignore the 00689 # RL LR for points with RR LL flagged! 00690 00691 print "" 00692 print "-----------------------------------------------------" 00693 print "Final cleanup of JUPITER data" 00694 print "Back to uvdist plot, see remaining bad data" 00695 print "You can draw little boxes around the outliers and Flag" 00696 print "Depends how patient you are in drawing boxes!" 00697 print "Could also use Locate to find where they come from" 00698 00699 # Pause script if you are running in scriptmode 00700 if scriptmode: 00701 user_check=raw_input('Return to continue script\n') 00702 00703 print "Done with plotxy!" 00704 00705 # 00706 #===================================================================== 00707 # 00708 # Use Flagmanager to save a copy of the flags so far 00709 # 00710 print '--Flagmanager--' 00711 default('flagmanager') 00712 00713 print "Now will use flagmanager to save a copy of the flags we just made" 00714 print "These are named xyflags" 00715 00716 vis = msfile 00717 mode = 'save' 00718 versionname = 'xyflags' 00719 comment = 'Plotxy flags' 00720 merge = 'replace' 00721 00722 flagmanager() 00723 00724 #===================================================================== 00725 # 00726 # Use Flagmanager to list all saved versions 00727 # 00728 print '--Flagmanager--' 00729 default('flagmanager') 00730 00731 print "Now will use flagmanager to list all the versions we saved" 00732 00733 vis = msfile 00734 mode = 'list' 00735 00736 flagmanager() 00737 00738 # 00739 # Done Flagging 00740 print '--Done with flagging--' 00741 00742 # 00743 #===================================================================== 00744 # Calibration 00745 #===================================================================== 00746 # 00747 # Set the fluxes of the primary calibrator(s) 00748 # 00749 print '--Setjy--' 00750 default('setjy') 00751 00752 print "Use setjy to set flux of 1331+305 (3C286)" 00753 00754 vis = msfile 00755 00756 # 00757 # 1331+305 = 3C286 is our primary calibrator 00758 field = '1331+305' 00759 00760 # Setjy knows about this source so we dont need anything more 00761 00762 setjy() 00763 00764 # 00765 # You should see something like this in the logger and casapy.log file: 00766 # 00767 # 1331+305 spwid= 0 [I=7.462, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00768 # 1331+305 spwid= 1 [I=7.51, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) 00769 # 00770 00771 print "Look in logger for the fluxes (should be 7.462 and 7.510 Jy)" 00772 00773 # 00774 #===================================================================== 00775 # 00776 # Initial gain calibration 00777 # 00778 print '--Gaincal--' 00779 default('gaincal') 00780 00781 print "Solve for antenna gains on 1331+305 and 0137+331" 00782 print "We have 2 single-channel continuum spw" 00783 print "Do not want bandpass calibration" 00784 00785 vis = msfile 00786 00787 # set the name for the output gain caltable 00788 caltable = gtable 00789 00790 print "Output gain cal table will be "+gtable 00791 00792 # Gain calibrators are 1331+305 and 0137+331 (FIELD_ID 7 and 0) 00793 # We have 2 IFs (SPW 0,1) with one channel each 00794 00795 # selection is via the field and spw strings 00796 field = '1331+305,0137+331' 00797 spw = '' 00798 00799 # a-priori calibration application 00800 gaincurve = usegaincurve 00801 opacity = gainopacity 00802 00803 # scan-based G solutions for both amplitude and phase 00804 gaintype = 'G' 00805 calmode = 'ap' 00806 00807 # one solution per scan 00808 solint = 'inf' 00809 combine = '' 00810 00811 # do not apply parallactic angle correction (yet) 00812 parang = False 00813 00814 # reference antenna 00815 refant = calrefant 00816 00817 # minimum SNR 3 00818 minsnr = 3 00819 00820 gaincal() 00821 00822 # 00823 #===================================================================== 00824 # 00825 # Bootstrap flux scale 00826 # 00827 print '--Fluxscale--' 00828 default('fluxscale') 00829 00830 print "Use fluxscale to rescale gain table to make new one" 00831 00832 vis = msfile 00833 00834 # set the name for the output rescaled caltable 00835 fluxtable = ftable 00836 00837 print "Output scaled gain cal table is "+ftable 00838 00839 # point to our first gain cal table 00840 caltable = gtable 00841 00842 # we will be using 1331+305 (the source we did setjy on) as 00843 # our flux standard reference 00844 reference = '1331+305' 00845 00846 # we want to transfer the flux to our other gain cal source 0137+331 00847 # to bring its gain amplitues in line with the absolute scale 00848 transfer = '0137+331' 00849 00850 fluxscale() 00851 00852 # You should see in the logger something like: 00853 #Flux density for 0137+331 in SpW=0 is: 00854 # 5.42575 +/- 0.00285011 (SNR = 1903.7, nAnt= 27) 00855 #Flux density for 0137+331 in SpW=1 is: 00856 # 5.46569 +/- 0.00301326 (SNR = 1813.88, nAnt= 27) 00857 00858 # 00859 #--------------------------------------------------------------------- 00860 # Plot calibration 00861 # 00862 print '--PlotCal--' 00863 default('plotcal') 00864 00865 showgui = True 00866 00867 caltable = ftable 00868 multiplot = True 00869 yaxis = 'amp' 00870 00871 showgui = True 00872 00873 plotcal() 00874 00875 print "" 00876 print "-------------------------------------------------" 00877 print "Plotcal" 00878 print "Looking at amplitude in cal-table "+caltable 00879 00880 # Pause script if you are running in scriptmode 00881 if scriptmode: 00882 user_check=raw_input('Return to continue script\n') 00883 00884 # 00885 # Now go back and plot to file 00886 # 00887 showgui = False 00888 00889 yaxis = 'amp' 00890 00891 figfile = caltable + '.plotcal.amp.png' 00892 print "Plotting calibration to file "+figfile 00893 #saveinputs('plotcal',caltable.plotcal.amp.saved') 00894 plotcal() 00895 00896 yaxis = 'phase' 00897 00898 figfile = caltable + '.plotcal.phase.png' 00899 print "Plotting calibration to file "+figfile 00900 #saveinputs('plotcal',caltable.plotcal.phase.saved') 00901 plotcal() 00902 00903 # 00904 #===================================================================== 00905 # Polarization Calibration 00906 #===================================================================== 00907 # 00908 if (dopolcal): 00909 print '--Polcal (D)--' 00910 default('polcal') 00911 00912 print "Solve for polarization leakage on 0137+331" 00913 print "Pretend it has unknown polarization" 00914 00915 vis = msfile 00916 00917 # Start with the un-fluxscaled gain table 00918 gaintable = gtable 00919 00920 # use settings from gaincal 00921 gaincurve = usegaincurve 00922 opacity = gainopacity 00923 00924 # Output table 00925 caltable = ptable 00926 00927 # Use a 3C48 tracked through a range of PA 00928 field = '0137+331' 00929 spw = '' 00930 00931 # No need for further selection 00932 selectdata=False 00933 00934 # Polcal mode (D+QU = unknown pol for D) 00935 poltype = 'D+QU' 00936 00937 # One solution for entire dataset 00938 solint = 'inf' 00939 combine = 'scan' 00940 00941 # reference antenna 00942 refant = calrefant 00943 00944 # minimum SNR 3 00945 minsnr = 3 00946 00947 #saveinputs('polcal',calprefix+'.polcal.saved') 00948 polcal() 00949 00950 #===================================================================== 00951 # 00952 # List polcal solutions 00953 # 00954 print '--Listcal (PolD)--' 00955 00956 listfile = caltable + '.list' 00957 00958 print "Listing calibration to file "+listfile 00959 00960 listcal() 00961 00962 #===================================================================== 00963 # 00964 # Plot polcal solutions 00965 # 00966 print '--Plotcal (PolD)--' 00967 00968 iteration = '' 00969 showgui = False 00970 00971 xaxis = 'real' 00972 yaxis = 'imag' 00973 figfile = caltable + '.plotcal.reim.png' 00974 print "Plotting calibration to file "+figfile 00975 #saveinputs('plotcal',caltable+'.plotcal.reim.saved') 00976 plotcal() 00977 00978 xaxis = 'antenna' 00979 yaxis = 'amp' 00980 figfile = caltable + '.plotcal.antamp.png' 00981 print "Plotting calibration to file "+figfile 00982 #saveinputs('plotcal',caltable+'.plotcal.antamp.saved') 00983 plotcal() 00984 00985 xaxis = 'antenna' 00986 yaxis = 'phase' 00987 figfile = caltable + '.plotcal.antphase.png' 00988 print "Plotting calibration to file "+figfile 00989 #saveinputs('plotcal',caltable+'.plotcal.antphase.saved') 00990 plotcal() 00991 00992 xaxis = 'antenna' 00993 yaxis = 'snr' 00994 figfile = caltable + '.plotcal.antsnr.png' 00995 print "Plotting calibration to file "+figfile 00996 #saveinputs('plotcal',caltable+'.plotcal.antsnr.saved') 00997 plotcal() 00998 00999 #===================================================================== 01000 # Do Chi (X) pol angle calibration 01001 #===================================================================== 01002 # First set the model 01003 print '--Setjy--' 01004 default('setjy') 01005 01006 vis = msfile 01007 01008 print "Use setjy to set IQU fluxes of "+polxfield 01009 field = polxfield 01010 01011 for spw in usespwlist: 01012 fluxdensity = polxiquv[spw] 01013 01014 #saveinputs('setjy',calprefix+'.setjy.polspw.'+spw+'.saved') 01015 setjy() 01016 01017 # 01018 # Polarization (X-term) calibration 01019 # 01020 print '--PolCal (X)--' 01021 default('polcal') 01022 01023 print "Polarization R-L Phase Calibration (linear approx)" 01024 01025 vis = msfile 01026 01027 # Start with the G and D tables 01028 gaintable = [gtable,ptable] 01029 01030 # use settings from gaincal 01031 gaincurve = usegaincurve 01032 opacity = gainopacity 01033 01034 # Output table 01035 caltable = xtable 01036 01037 # previously set with setjy 01038 field = polxfield 01039 spw = '' 01040 01041 selectdata=False 01042 01043 # Solve for Chi 01044 poltype = 'X' 01045 solint = 'inf' 01046 combine = 'scan' 01047 01048 # reference antenna 01049 refant = calrefant 01050 01051 # minimum SNR 3 01052 minsnr = 3 01053 01054 #saveinputs('polcal',calprefix+'.polcal.X.saved') 01055 polcal() 01056 01057 #===================================================================== 01058 # Apply the Calibration 01059 #===================================================================== 01060 # 01061 # Interpolate the gains onto Jupiter (and others) 01062 # 01063 # print '--Accum--' 01064 # default('accum') 01065 # 01066 # print "This will interpolate the gains onto Jupiter" 01067 # 01068 # vis = msfile 01069 # 01070 # tablein = '' 01071 # incrtable = ftable 01072 # calfield = '1331+305, 0137+331' 01073 # 01074 # # set the name for the output interpolated caltable 01075 # caltable = atable 01076 # 01077 # print "Output cumulative gain table will be "+atable 01078 # 01079 # # linear interpolation 01080 # interp = 'linear' 01081 # 01082 # # make 10s entries 01083 # accumtime = 10.0 01084 # 01085 # accum() 01086 # 01087 # NOTE: bypassing this during testing 01088 atable = ftable 01089 01090 # #===================================================================== 01091 # 01092 # Correct the data 01093 # (This will put calibrated data into the CORRECTED_DATA column) 01094 # 01095 print '--ApplyCal--' 01096 default('applycal') 01097 01098 print "This will apply the calibration to the DATA" 01099 print "Fills CORRECTED_DATA" 01100 01101 vis = msfile 01102 01103 # Start with the interpolated fluxscale/gain table 01104 gaintable = [atable,ptable,xtable] 01105 01106 # use settings from gaincal 01107 gaincurve = usegaincurve 01108 opacity = gainopacity 01109 01110 # select the fields 01111 field = '1331+305,0137+331,JUPITER' 01112 spw = '' 01113 selectdata = False 01114 01115 # IMPORTANT set parang=True for polarization 01116 parang = True 01117 01118 # do not need to select subset since we did accum 01119 # (note that correct only does 'nearest' interp) 01120 gainfield = '' 01121 01122 applycal() 01123 01124 # 01125 #===================================================================== 01126 # 01127 # Now split the Jupiter target data 01128 # 01129 print '--Split Jupiter--' 01130 default('split') 01131 01132 vis = msfile 01133 01134 # Now we write out the corrected data to a new MS 01135 01136 # Select the Jupiter field 01137 field = srcname 01138 spw = '' 01139 01140 # pick off the CORRECTED_DATA column 01141 datacolumn = 'corrected' 01142 01143 # Make an output vis file 01144 outputvis = srcsplitms 01145 01146 print "Split "+field+" data into new ms "+srcsplitms 01147 01148 split() 01149 01150 # Also split out 0137+331 as a check 01151 field = calname 01152 01153 outputvis = calsplitms 01154 01155 print "Split "+field+" data into new ms "+calsplitms 01156 01157 split() 01158 01159 #===================================================================== 01160 # Force scratch column creation so plotxy will work 01161 # 01162 vis = srcsplitms 01163 clearcal() 01164 01165 vis = calsplitms 01166 clearcal() 01167 01168 #===================================================================== 01169 # Use Plotxy to look at the split calibrated data 01170 # 01171 print '--Plotxy--' 01172 default('plotxy') 01173 01174 vis = srcsplitms 01175 selectdata = True 01176 01177 # Plot only the RR and LL for now 01178 correlation = 'RR LL' 01179 01180 # Plot amplitude vs. uvdist 01181 xaxis = 'uvdist' 01182 datacolumn = 'data' 01183 multicolor = 'both' 01184 01185 iteration = '' 01186 selectplot = True 01187 interactive = True 01188 01189 field = 'JUPITER' 01190 yaxis = 'amp' 01191 # Use the field name as the title 01192 title = field+" " 01193 01194 plotxy() 01195 01196 print "" 01197 print "-----------------------------------------------------" 01198 print "Plotting JUPITER corrected visibilities" 01199 print "Look for outliers" 01200 01201 # Pause script if you are running in scriptmode 01202 if scriptmode: 01203 user_check=raw_input('Return to continue script\n') 01204 01205 # Now go back and plot to files 01206 interactive = False 01207 01208 # 01209 # First the target 01210 # 01211 vis = srcsplitms 01212 field = srcname 01213 yaxis = 'amp' 01214 # Use the field name as the title 01215 title = field+" " 01216 figfile = vis + '.plotxy.amp.png' 01217 print "Plotting to file "+figfile 01218 #saveinputs('plotxy',vis+'.plotxy.amp.saved') 01219 01220 plotxy() 01221 01222 yaxis = 'phase' 01223 # Use the field name as the title 01224 figfile = vis + '.plotxy.phase.png' 01225 print "Plotting to file "+figfile 01226 #saveinputs('plotxy',vis+'.plotxy.phase.saved') 01227 01228 plotxy() 01229 01230 # 01231 # Now the calibrator 01232 # 01233 vis = calsplitms 01234 field = calname 01235 yaxis = 'amp' 01236 # Use the field name as the title 01237 title = field+" " 01238 figfile = vis + '.plotxy.amp.png' 01239 print "Plotting to file "+figfile 01240 #saveinputs('plotxy',vis+'.plotxy.amp.saved') 01241 01242 plotxy() 01243 01244 yaxis = 'phase' 01245 # Use the field name as the title 01246 figfile = vis + '.plotxy.phase.png' 01247 print "Plotting to file "+figfile 01248 #saveinputs('plotxy',vis+'.plotxy.phase.saved') 01249 01250 plotxy() 01251 01252 print 'Calibration completed' 01253 # 01254 #===================================================================== 01255 # 01256 # Intensity Imaging/Selfcal 01257 # 01258 #===================================================================== 01259 # 01260 # Make the scratch columns in the split ms 01261 # 01262 print '--Clearcal--' 01263 default('clearcal') 01264 01265 vis = srcsplitms 01266 01267 clearcal() 01268 01269 print "Created scratch columns for MS "+vis 01270 print "" 01271 # 01272 #===================================================================== 01273 # FIRST CLEAN / SELFCAL CYCLE 01274 #===================================================================== 01275 # 01276 # Now clean an image of Jupiter 01277 # NOTE: this uses the new combined invert/clean/mosaic task Patch 2 01278 # 01279 print '--Clean 1--' 01280 default('clean') 01281 01282 # Pick up our split source data 01283 vis = srcsplitms 01284 01285 # Make an image root file name 01286 imagename = imname1 01287 01288 print "Output images will be prefixed with "+imname1 01289 01290 # Set up the output continuum image (single plane mfs) 01291 mode = 'mfs' 01292 stokes = 'I' 01293 01294 print "Will be a single MFS continuum image" 01295 01296 # NOTE: current version field='' doesnt work 01297 field = '*' 01298 01299 # Combine all spw 01300 spw = '' 01301 01302 # Imaging mode params 01303 psfmode = clnalg 01304 imagermode = clnmode 01305 01306 # Imsize and cell 01307 imsize = clnimsize 01308 cell = clncell 01309 01310 # NOTE: will eventually have an imadvise task to give you this 01311 # information 01312 01313 # Standard gain factor 0.1 01314 gain = 0.1 01315 01316 # Fix maximum number of iterations and threshold 01317 niter = clniter 01318 threshold = clnthreshold 01319 01320 # Note - we can change niter and threshold interactively 01321 # during clean 01322 01323 # Set up the weighting 01324 # Use Briggs weighting (a moderate value, on the uniform side) 01325 weighting = 'briggs' 01326 robust = 0.5 01327 01328 # No clean mask or box 01329 mask = '' 01330 01331 # Use interactive clean mode 01332 interactive = True 01333 01334 # Moderate number of iter per interactive cycle 01335 npercycle = 100 01336 01337 saveinputs('clean',imagename+'.clean.saved') 01338 clean() 01339 01340 # When the interactive clean window comes up, use the right-mouse 01341 # to draw rectangles around obvious emission double-right-clicking 01342 # inside them to add to the flag region. You can also assign the 01343 # right-mouse to polygon region drawing by right-clicking on the 01344 # polygon drawing icon in the toolbar. When you are happy with 01345 # the region, click 'Done Flagging' and it will go and clean another 01346 # 100 iterations. When done, click 'Stop'. 01347 01348 print "" 01349 print "----------------------------------------------------" 01350 print "Clean" 01351 print "Final clean model is "+clnmodel1 01352 print "Final restored clean image is "+clnimage1 01353 print "The clean residual image is "+clnresid1 01354 print "Your final clean mask is "+clnmask1 01355 01356 print "" 01357 print "This is the final restored clean image in the viewer" 01358 print "Zoom in and set levels to see faint emission" 01359 print "Use rectangle drawing tool to box off source" 01360 print "Double-click inside to print statistics" 01361 print "Move box on-source and get the max" 01362 print "Calcualte DynRange = MAXon/RMSoff" 01363 print "I got 1.060/0.004 = 270" 01364 print "Still not as good as it can be - lets selfcal" 01365 print "Close viewer panel when done" 01366 01367 # 01368 #--------------------------------------------------------------------- 01369 # 01370 # If you did not do interactive clean, bring up viewer manually 01371 viewer(clnimage1,'image') 01372 01373 # Pause script if you are running in scriptmode 01374 if scriptmode: 01375 user_check=raw_input('Return to continue script\n') 01376 01377 # You can use the right-mouse to draw a box in the lower right 01378 # corner of the image away from emission, the double-click inside 01379 # to bring up statistics. Use the right-mouse to grab this box 01380 # and move it up over Jupiter and double-click again. You should 01381 # see stuff like this in the terminal: 01382 # 01383 # jupiter6cm.usecase.clean1.image (Jy/beam) 01384 # 01385 # n Std Dev RMS Mean Variance Sum 01386 # 4712 0.003914 0.003927 0.0003205 1.532e-05 1.510 01387 # 01388 # Flux Med |Dev| IntQtlRng Median Min Max 01389 # 0.09417 0.002646 0.005294 0.0001885 -0.01125 0.01503 01390 # 01391 # 01392 # On Jupiter: 01393 # 01394 # n Std Dev RMS Mean Variance Sum 01395 # 3640 0.1007 0.1027 0.02023 0.01015 73.63 01396 # 01397 # Flux Med |Dev| IntQtlRng Median Min Max 01398 # 4.592 0.003239 0.007120 0.0001329 -0.01396 1.060 01399 # 01400 # Estimated dynamic range = 1.060 / 0.003927 = 270 (poor) 01401 # 01402 # Note that the exact numbers you get will depend on how deep you 01403 # take the interactive clean and how you draw the box for the stats. 01404 01405 #===================================================================== 01406 # 01407 # Do some non-interactive image statistics 01408 print '--Imstat--' 01409 default('imstat') 01410 01411 imagename = clnimage1 01412 on_statistics1 = imstat() 01413 01414 # Now do stats in the lower right corner of the image 01415 # remember clnimsize = [288,288] 01416 box = '216,1,287,72' 01417 off_statistics1 = imstat() 01418 01419 # Pull the max and rms from the clean image 01420 thistest_immax=on_statistics1['max'][0] 01421 print ' Found : Max in image = ',thistest_immax 01422 thistest_imrms=off_statistics1['rms'][0] 01423 print ' Found : rms in image = ',thistest_imrms 01424 print ' Clean image Dynamic Range = ',thistest_immax/thistest_imrms 01425 print '' 01426 # 01427 #--------------------------------------------------------------------- 01428 # 01429 # Self-cal using clean model 01430 # 01431 # Note: clean will have left FT of model in the MODEL_DATA column 01432 # If you've done something in between, can use the ft task to 01433 # do this manually. 01434 # 01435 print '--SelfCal 1--' 01436 default('gaincal') 01437 01438 vis = srcsplitms 01439 01440 print "Will self-cal using MODEL_DATA left in MS by clean" 01441 01442 # New gain table 01443 caltable = selfcaltab1 01444 01445 print "Will write gain table "+selfcaltab1 01446 01447 # Don't need a-priori cals 01448 selectdata = False 01449 gaincurve = False 01450 opacity = 0.0 01451 01452 # This choice seemed to work 01453 refant = calrefant 01454 01455 # Do amp and phase 01456 gaintype = 'G' 01457 calmode = 'ap' 01458 01459 # Do 30s solutions with SNR>1 01460 solint = 30.0 01461 minsnr = 1.0 01462 print "Calibrating amplitudes and phases on 30s timescale" 01463 01464 # Do not need to normalize (let gains float) 01465 solnorm = False 01466 01467 gaincal() 01468 01469 # 01470 #--------------------------------------------------------------------- 01471 # It is useful to put this up in plotcal 01472 # 01473 # 01474 print '--PlotCal--' 01475 default('plotcal') 01476 01477 caltable = selfcaltab1 01478 multiplot = True 01479 yaxis = 'amp' 01480 01481 plotcal() 01482 01483 print "" 01484 print "-------------------------------------------------" 01485 print "Plotcal" 01486 print "Looking at amplitude in self-cal table "+caltable 01487 01488 # Pause script if you are running in scriptmode 01489 if scriptmode: 01490 user_check=raw_input('Return to continue script\n') 01491 01492 yaxis = 'phase' 01493 01494 plotcal() 01495 01496 print "" 01497 print "-------------------------------------------------" 01498 print "Plotcal" 01499 print "Looking at phases in self-cal table "+caltable 01500 01501 # 01502 # Pause script if you are running in scriptmode 01503 if scriptmode: 01504 user_check=raw_input('Return to continue script\n') 01505 01506 # 01507 #--------------------------------------------------------------------- 01508 # 01509 # Correct the data (no need for interpolation this stage) 01510 # 01511 print '--ApplyCal--' 01512 default('applycal') 01513 01514 vis = srcsplitms 01515 01516 print "Will apply self-cal table to over-write CORRECTED_DATA in MS" 01517 01518 gaintable = selfcaltab1 01519 01520 gaincurve = False 01521 opacity = 0.0 01522 field = '' 01523 spw = '' 01524 selectdata = False 01525 01526 calwt = True 01527 01528 applycal() 01529 01530 # Self-cal is now in CORRECTED_DATA column of split ms 01531 #===================================================================== 01532 # Use Plotxy to look at the self-calibrated data 01533 # 01534 print '--Plotxy--' 01535 default('plotxy') 01536 01537 vis = srcsplitms 01538 selectdata = True 01539 field = 'JUPITER' 01540 01541 correlation = 'RR LL' 01542 xaxis = 'uvdist' 01543 yaxis = 'amp' 01544 datacolumn = 'corrected' 01545 multicolor = 'both' 01546 01547 # Use the field name as the title 01548 selectplot = True 01549 title = field+" " 01550 01551 iteration = '' 01552 01553 plotxy() 01554 01555 print "" 01556 print "-----------------------------------------------------" 01557 print "Plotting JUPITER self-corrected visibilities" 01558 print "Look for outliers, and you can flag them" 01559 01560 # Pause script if you are running in scriptmode 01561 if scriptmode: 01562 user_check=raw_input('Return to continue script\n') 01563 01564 # 01565 #===================================================================== 01566 # SECOND CLEAN / SELFCAL CYCLE 01567 #===================================================================== 01568 # 01569 print '--Clean 2--' 01570 default('clean') 01571 01572 print "Now clean on self-calibrated data" 01573 01574 vis = srcsplitms 01575 01576 imagename = imname2 01577 01578 field = '*' 01579 spw = '' 01580 mode = 'mfs' 01581 gain = 0.1 01582 01583 # Imaging mode params 01584 psfmode = clnalg 01585 imagermode = clnmode 01586 imsize = clnimsize 01587 cell = clncell 01588 niter = clniter 01589 threshold = clnthreshold 01590 01591 weighting = 'briggs' 01592 robust = 0.5 01593 01594 mask = '' 01595 interactive = True 01596 npercycle = 100 01597 01598 saveinputs('clean',imagename+'.clean.saved') 01599 clean() 01600 01601 print "" 01602 print "----------------------------------------------------" 01603 print "Clean" 01604 print "Final clean model is "+clnmodel2 01605 print "Final restored clean image is "+clnimage2 01606 print "The clean residual image is "+clnresid2 01607 print "Your final clean mask is "+clnmask2 01608 01609 print "" 01610 print "This is the final restored clean image in the viewer" 01611 print "Zoom in and set levels to see faint emission" 01612 print "Use rectangle drawing tool to box off source" 01613 print "Double-click inside to print statistics" 01614 print "Move box on-source and get the max" 01615 print "Calcualte DynRange = MAXon/RMSoff" 01616 print "This time I got 1.076 / 0.001389 = 775 (better)" 01617 print "Still not as good as it can be - lets selfcal again" 01618 print "Close viewer panel when done" 01619 01620 # 01621 #--------------------------------------------------------------------- 01622 # 01623 # If you did not do interactive clean, bring up viewer manually 01624 viewer(clnimage2,'image') 01625 01626 # Pause script if you are running in scriptmode 01627 if scriptmode: 01628 user_check=raw_input('Return to continue script\n') 01629 01630 # jupiter6cm.usecase.clean2.image (Jy/beam) 01631 # 01632 # n Std Dev RMS Mean Variance Sum 01633 # 5236 0.001389 0.001390 3.244e-05 1.930e-06 0.1699 01634 # 01635 # Flux Med |Dev| IntQtlRng Median Min Max 01636 # 0.01060 0.0009064 0.001823 -1.884e-05 -0.004015 0.004892 01637 # 01638 # 01639 # On Jupiter: 01640 # 01641 # n Std Dev RMS Mean Variance Sum 01642 # 5304 0.08512 0.08629 0.01418 0.007245 75.21 01643 # 01644 # Flux Med |Dev| IntQtlRng Median Min Max 01645 # 4.695 0.0008142 0.001657 0.0001557 -0.004526 1.076 01646 # 01647 # Estimated dynamic range = 1.076 / 0.001389 = 775 (better) 01648 # 01649 # Note that the exact numbers you get will depend on how deep you 01650 # take the interactive clean and how you draw the box for the stats. 01651 # 01652 print "" 01653 print "--------------------------------------------------" 01654 print "After this script is done you can continue on with" 01655 print "more self-cal, or try different cleaning options" 01656 01657 # 01658 #===================================================================== 01659 # Image Analysis 01660 #===================================================================== 01661 # 01662 # Can do some image statistics if you wish 01663 print '--Imstat (Cycle 2)--' 01664 default('imstat') 01665 01666 imagename = clnimage2 01667 on_statistics2 = imstat() 01668 01669 # Now do stats in the lower right corner of the image 01670 # remember clnimsize = [288,288] 01671 box = '216,1,287,72' 01672 off_statistics2 = imstat() 01673 01674 # Pull the max and rms from the clean image 01675 thistest_immax=on_statistics2['max'][0] 01676 print ' Found : Max in image = ',thistest_immax 01677 thistest_imrms=off_statistics2['rms'][0] 01678 print ' Found : rms in image = ',thistest_imrms 01679 print ' Clean image Dynamic Range = ',thistest_immax/thistest_imrms 01680 print '' 01681 01682 #===================================================================== 01683 # 01684 # Print results and regression versus previous runs 01685 # 01686 print "" 01687 print ' Final Jupiter results ' 01688 print ' ===================== ' 01689 print '' 01690 # Pull the max and rms from the clean image 01691 thistest_immax=on_statistics2['max'][0] 01692 oldtest_immax = 1.07732224464 01693 print ' Clean image ON-SRC max = ',thistest_immax 01694 print ' Previously found to be = ',oldtest_immax 01695 diff_immax = abs((oldtest_immax-thistest_immax)/oldtest_immax) 01696 print ' Difference (fractional) = ',diff_immax 01697 01698 print '' 01699 thistest_imrms=off_statistics2['rms'][0] 01700 oldtest_imrms = 0.0010449 01701 print ' Clean image OFF-SRC rms = ',thistest_imrms 01702 print ' Previously found to be = ',oldtest_imrms 01703 diff_imrms = abs((oldtest_imrms-thistest_imrms)/oldtest_imrms) 01704 print ' Difference (fractional) = ',diff_imrms 01705 01706 print '' 01707 print ' Final Clean image Dynamic Range = ',thistest_immax/thistest_imrms 01708 print '' 01709 print '--- Done with I Imaging and Selfcal---' 01710 01711 # 01712 #===================================================================== 01713 # Polarization Imaging 01714 #===================================================================== 01715 # 01716 print '--Clean (Polarization)--' 01717 default('clean') 01718 01719 print "Now clean polarized data" 01720 01721 vis = srcsplitms 01722 01723 imagename = polimname 01724 01725 field = '*' 01726 spw = '' 01727 mode = 'mfs' 01728 gain = 0.1 01729 01730 # Polarization 01731 stokes = 'IQUV' 01732 01733 psfmode = polclnalg 01734 imagermode = polclnmode 01735 01736 niter = clniter 01737 threshold = clnthreshold 01738 01739 imsize = clnimsize 01740 cell = clncell 01741 01742 weighting = 'briggs' 01743 robust = 0.5 01744 01745 interactive = True 01746 npercycle = 100 01747 01748 saveinputs('clean',imagename+'.clean.saved') 01749 clean() 01750 01751 print "" 01752 print "----------------------------------------------------" 01753 print "Clean" 01754 print "Final restored clean image is "+polimage 01755 print "Final clean model is "+polmodel 01756 print "The clean residual image is "+polresid 01757 print "Your final clean mask is "+polmask 01758 01759 # 01760 #===================================================================== 01761 # Image Analysis 01762 #===================================================================== 01763 # 01764 # Polarization statistics 01765 print '--Final Pol Imstat--' 01766 default('imstat') 01767 01768 imagename = polimage 01769 01770 on_statistics = {} 01771 off_statistics = {} 01772 01773 # lower right corner of the image (clnimsize = [288,288]) 01774 onbox = '' 01775 # lower right corner of the image (clnimsize = [288,288]) 01776 offbox = '216,1,287,72' 01777 01778 for stokes in ['I','Q','U','V']: 01779 box = onbox 01780 on_statistics[stokes] = imstat() 01781 box = offbox 01782 off_statistics[stokes] = imstat() 01783 01784 # 01785 # Peel off some Q and U planes 01786 # 01787 print '--Immath--' 01788 default('immath') 01789 01790 mode = 'evalexpr' 01791 01792 stokes = 'I' 01793 outfile = ipolimage 01794 expr = '\"'+polimage+'\"' 01795 01796 immath() 01797 print "Created I image "+outfile 01798 01799 stokes = 'Q' 01800 outfile = qpolimage 01801 expr = '\"'+polimage+'\"' 01802 01803 immath() 01804 print "Created Q image "+outfile 01805 01806 stokes = 'U' 01807 outfile = upolimage 01808 expr = '\"'+polimage+'\"' 01809 01810 immath() 01811 print "Created U image "+outfile 01812 01813 # 01814 #--------------------------------------------------------------------- 01815 # Now make POLI and POLA images 01816 # 01817 stokes = '' 01818 outfile = poliimage 01819 mode = 'poli' 01820 imagename = [qpolimage,upolimage] 01821 # Use our rms above for debiasing 01822 mysigma = 0.5*( off_statistics['Q']['rms'][0] + off_statistics['U']['rms'][0] ) 01823 #sigma = str(mysigma)+'Jy/beam' 01824 # This does not work well yet 01825 sigma = '0.0Jy/beam' 01826 01827 immath() 01828 print "Created POLI image "+outfile 01829 01830 outfile = polaimage 01831 mode = 'pola' 01832 01833 immath() 01834 print "Created POLA image "+outfile 01835 01836 # 01837 #--------------------------------------------------------------------- 01838 # Save statistics of these images 01839 default('imstat') 01840 01841 imagename = poliimage 01842 stokes = '' 01843 box = onbox 01844 on_statistics['POLI'] = imstat() 01845 box = offbox 01846 off_statistics['POLI'] = imstat() 01847 01848 # 01849 # 01850 #--------------------------------------------------------------------- 01851 # Display clean I image in viewer but with polarization vectors 01852 # 01853 # If you did not do interactive clean, bring up viewer manually 01854 viewer(polimage,'image') 01855 01856 print "Displaying pol I now. You should overlay pola vectors" 01857 print "Bring up the Load Data panel:" 01858 print "" 01859 print "Use LEL for POLA VECTOR with cut above 6*mysigma in POLI = "+str(6*mysigma) 01860 print "For example:" 01861 print "\'"+polaimage+"\'[\'"+poliimage+"\'>0.0048]" 01862 print "" 01863 print "In the Data Display Options for the vector plot:" 01864 print " Set the x,y increments to 2 (default is 3)" 01865 print " Use an extra rotation this 90deg to get B field" 01866 print "Note the lengths are all equal. You can fiddle these." 01867 print "" 01868 print "You can also load the poli image as contours" 01869 01870 # Pause script if you are running in scriptmode 01871 if scriptmode: 01872 user_check=raw_input('Return to continue script\n') 01873 01874 # NOTE: the LEL will be something like 01875 # 'jupiter6cm.usecase.polimg.clean.image.pola'['jupiter6cm.usecase.polimg.clean.image.poli'>0.005] 01876 01877 # 01878 # NOTE: The viewer can take complex images to make Vector plots, although 01879 # the image analysis tasks (and ia tool) cannot yet handle these. But we 01880 # can use the imagepol tool (which is not imported by default) to make 01881 # a complex image of the linear polarized intensity for display. 01882 # See CASA User Reference Manual: 01883 # http://casa.nrao.edu/docs/casaref/imagepol-Tool.html 01884 # 01885 # Make an imagepol tool and open the clean image 01886 po = casac.imagepol() 01887 po.open(polimage) 01888 # Use complexlinpol to make a Q+iU image 01889 complexlinpolimage = polimname + '.cmplxlinpol' 01890 po.complexlinpol(complexlinpolimage) 01891 po.close() 01892 01893 # You can now display this in the viewer, in particular overlay this 01894 # over the intensity raster with the poli contours. The vector lengths 01895 # will be proportional to the polarized intensity. You can play with 01896 # the Data Display Options panel for vector spacing and length. 01897 # You will want to have this masked, like the pola image above, on 01898 # the polarized intensity. When you load the image, use the LEL: 01899 # 'jupiter6cm.usecase.polimg.clean.cmplxlinpol'['jupiter6cm.usecase.polimg.clean.image.poli'>0.005] 01900 01901 #===================================================================== 01902 # 01903 # Print results 01904 # 01905 print "" 01906 print ' Jupiter polarization results ' 01907 print ' ============================ ' 01908 print '' 01909 for stokes in ['I','Q','U','V','POLI']: 01910 print '' 01911 print ' =============== ' 01912 print '' 01913 print ' Polarization (Stokes '+stokes+'):' 01914 mymax = on_statistics[stokes]['max'][0] 01915 mymin = on_statistics[stokes]['min'][0] 01916 myrms = off_statistics[stokes]['rms'][0] 01917 absmax = max(mymax,mymin) 01918 mydra = absmax/myrms 01919 print ' Clean image ON-SRC max = ',mymax 01920 print ' Clean image ON-SRC min = ',mymin 01921 print ' Clean image OFF-SRC rms = ',myrms 01922 print ' Clean image dynamic rng = ',mydra 01923 01924 01925 print '--- Done ---' 01926 01927 # 01928 #=====================================================================
1.8.0
