Testing Date - 2002.05.23 through 2002.05.30 Tester - S.T. Myers (AOC) Platform - Linux (sandrock) Last Revised - 2002.05.31 Note: bugs are denoted by lines with prefix: >>>BUG: and queries are denoted by lines with prefix: >>>QUERY: while comments are prefixed by: >>>COMMENT: Some enhancement requests are given after >>>REQUEST: Notes added later for clarification are indicated by >>>NOTE: ------------------------------------------------------------------------------- GOAL: start going through tutorial jupiter dataset =============================================================================== Date - 2002.05.23-24 Version - Weekly (1.7 build 373) Directory - /home/sandrock/smyers/Testing/2002-05-23/ Data - /aips++/data/nrao/VLA/planets_6cm.fits =============================================================================== George's script is jupiter.g for reference. I will attempt to use the GUI to follow through the VLA-specific cookbook procedure on this data. This is what the students will have to work from. The discussion is meant to be representative of what I would say to the students. There are numerous comments/queries denoted by >>>BUG: here, you might want to search through this looking at them, not all were submitted as bugs but are for NAUG discussion etc. ------------------------------------------------------------------------------- First - give quick orientation on how to start up aips++ they may well kill it inadvertantly during the session) and to navigate the Tool Manager, both Manager and Tools in use. Also bring up a Netscape and point to the documentation tree, show them how to get to the VLA cookbook section of in Volume 3 of Getting Results, Volume 1 of Getting Results, and the User Reference Guide. Describe some glish functionality (show manual?) and explain GUI vs scripting. And say what a measurement set is and maybe look at one in the catalog tool. We will use the GUIs from the Tool manager window to carry out the commands. The cookbook describes the procedures in terms of the aips++ glish commands as if you were running from a script. One useful thing is to turn on the scripter window, to see the commands that punching GUI buttons is actually executing. To do this, at the top of the Tool manager window, choose 'Show scripter' from Options menu. A scripter window should pop up. Note that the 'Copy commands to scripter' Option should be selected by default, if not choose it now. Now you should see the commands appear in the scripter window as you exercise the GUI. Note that you can Save this script for later use (it is executable from glish)! Now we follow through the delineated steps in the cookbook http://aips2.nrao.edu/stable/docs/gettingresults/grvol3/node2.html The numbered steps are taken directly from the cookbook, with the sections noted, e.g. >1.5.2 Filling from UVFITS format ------------------------------------------------------------------------------- 1. Import the data. Data may be imported from either a VLA archive tape using the vlafiller tool, or from a UVFITS disk file (e.g., as written by AIPS) using the ms tool's fitstoms constructor. The result is a Measurement Set table (the AIPS++ internal data format). In either case, the ms tool is useful for obtaining a basic summary of the dataset. >1.5.2 Filling from UVFITS format Use FITS dataset '/aips++/data/nrao/VLA/planets_6cm.fits' Create: general->ms->ms Construct: fitstoms->myms( fitsfile='/aips++/data/nrao/VLA/planets_6cm.fits', msfile='planets.ms') Watch the progress meter window and the logger, and choose the summary function, then Go myms.summary(verbose=T) Watch the logger and the red 'Running: abort?' button to see that it is in progress. Output does not appear in the glish window but in the logger. You should keep an online log, you can cut and paste from the logger window, e.g. ObservationID = 1 ArrayID = 1 Date Timerange Field DataDescIds 15-Apr-1999/23:15:26.7 - 23:16:10.0 0137+331 [1, 2] 23:38:40.0 - 23:48:00.0 0813+482 [1, 2] 23:53:40.0 - 23:55:20.0 0542+498 [1, 2] 16-Apr-1999/00:22:10.1 - 00:23:49.9 0437+296 [1, 2] 00:28:23.3 - 00:30:00.1 VENUS [1, 2] 00:48:40.0 - 00:50:20.0 0813+482 [1, 2] 00:56:13.4 - 00:57:49.9 0542+498 [1, 2] 01:10:20.1 - 01:11:59.9 0521+166 [1, 2] 01:23:29.9 - 01:25:00.1 0437+296 [1, 2] 01:29:33.3 - 01:31:10.0 VENUS [1, 2] 01:49:50.0 - 01:51:30.0 1411+522 [1, 2] 02:03:00.0 - 02:04:30.0 1331+305 [1, 2] 02:17:30.0 - 02:19:10.0 0813+482 [1, 2] 02:24:20.0 - 02:26:00.0 0542+498 [1, 2] 02:37:49.9 - 02:39:30.0 0521+166 [1, 2] 02:50:50.1 - 02:52:20.1 0437+296 [1, 2] 02:59:20.0 - 03:01:00.0 1411+522 [1, 2] 03:12:30.0 - 03:14:10.0 1331+305 [1, 2] 03:27:53.3 - 03:29:39.9 0813+482 [1, 2] 03:35:00.0 - 03:36:40.0 0542+498 [1, 2] 03:49:50.0 - 03:51:30.1 1411+522 [1, 2] 04:03:10.0 - 04:04:50.0 1331+305 [1, 2] 04:18:49.9 - 04:20:40.0 0813+482 [1, 2] 04:25:56.6 - 04:27:39.9 0542+498 [1, 2] 04:42:49.9 - 04:44:40.0 MARS [1, 2] 04:56:50.0 - 04:58:30.1 1411+522 [1, 2] 05:24:03.3 - 05:33:39.9 1331+305 [1, 2] 05:48:00.0 - 05:49:49.9 0813+482 [1, 2] 05:58:36.6 - 06:00:30.0 MARS [1, 2] 06:13:20.1 - 06:14:59.9 1411+522 [1, 2] 06:27:40.0 - 06:29:20.0 1331+305 [1, 2] 06:44:13.4 - 06:46:00.0 0813+482 [1, 2] 06:55:06.6 - 06:57:00.0 MARS [1, 2] 07:10:40.0 - 07:12:20.0 1411+522 [1, 2] 07:28:20.0 - 07:30:10.1 1331+305 [1, 2] 07:42:49.9 - 07:44:30.0 MARS [1, 2] 07:58:43.3 - 08:00:39.9 1411+522 [1, 2] 08:13:30.0 - 08:15:19.9 1331+305 [1, 2] 08:27:53.4 - 08:29:30.0 MARS [1, 2] 08:42:59.9 - 08:44:50.0 1411+522 [1, 2] 08:57:09.9 - 08:58:50.0 1331+305 [1, 2] 09:13:03.3 - 09:14:50.1 NGC7027 [1, 2] 09:26:59.9 - 09:28:40.0 1411+522 [1, 2] 09:40:33.4 - 09:42:09.9 1331+305 [1, 2] 09:56:19.9 - 09:58:10.0 NGC7027 [1, 2] 10:12:59.9 - 10:14:50.0 MARS [1, 2] 10:27:09.9 - 10:28:50.0 1411+522 [1, 2] 10:40:30.0 - 10:42:00.0 1331+305 [1, 2] 10:56:10.0 - 10:57:50.0 NGC7027 [1, 2] 11:28:30.0 - 11:35:30.0 NEPTUNE [1, 2] 11:48:20.0 - 11:50:10.0 1411+522 [1, 2] 12:01:36.7 - 12:03:10.0 1331+305 [1, 2] 12:35:33.3 - 12:37:40.0 URANUS [1, 2] 12:46:30.0 - 12:48:10.0 NEPTUNE [1, 2] 13:00:29.9 - 13:02:10.0 1411+522 [1, 2] 13:15:23.3 - 13:17:10.1 NGC7027 [1, 2] 13:33:43.3 - 13:35:40.0 URANUS [1, 2] 13:44:30.0 - 13:46:10.0 NEPTUNE [1, 2] 14:00:46.7 - 14:01:39.9 0137+331 [1, 2] 14:10:40.0 - 14:12:09.9 JUPITER [1, 2] 14:24:06.6 - 14:25:40.1 URANUS [1, 2] 14:34:30.0 - 14:36:10.1 NEPTUNE [1, 2] 14:59:13.4 - 15:00:00.0 0137+331 [1, 2] 15:09:03.3 - 15:10:40.1 JUPITER [1, 2] 15:24:30.0 - 15:26:20.1 NGC7027 [1, 2] 15:40:10.0 - 15:45:00.0 URANUS [1, 2] 15:53:50.0 - 15:55:20.0 NEPTUNE [1, 2] 16:18:53.4 - 16:19:49.9 0137+331 [1, 2] 16:29:10.1 - 16:30:49.9 JUPITER [1, 2] 16:42:53.4 - 16:44:30.0 URANUS [1, 2] 16:54:53.4 - 16:56:40.0 NGC7027 [1, 2] 17:23:06.6 - 17:30:40.0 0542+498 [1, 2] 17:41:50.0 - 17:43:20.0 0437+296 [1, 2] 17:55:36.7 - 17:57:39.9 VENUS [1, 2] 18:19:23.3 - 18:20:09.9 0137+331 [1, 2] 18:30:23.3 - 18:32:00.0 JUPITER [1, 2] 18:44:49.9 - 18:46:30.0 NGC7027 [1, 2] 18:59:13.3 - 19:00:59.9 0542+498 [1, 2] 19:19:10.0 - 19:21:20.1 0521+166 [1, 2] 19:32:50.1 - 19:34:29.9 0437+296 [1, 2] 19:39:03.3 - 19:40:40.1 VENUS [1, 2] 20:08:06.7 - 20:08:59.9 0137+331 [1, 2] 20:18:10.0 - 20:19:50.0 JUPITER [1, 2] 20:33:53.3 - 20:35:40.1 0813+482 [1, 2] 20:40:59.9 - 20:42:40.0 0542+498 [1, 2] 21:00:16.6 - 21:02:20.1 0521+166 [1, 2] 21:13:53.4 - 21:15:29.9 0437+296 [1, 2] 21:20:43.4 - 21:22:30.0 VENUS [1, 2] 21:47:26.7 - 21:48:20.1 0137+331 [1, 2] 21:57:30.0 - 21:59:10.0 JUPITER [1, 2] 22:12:13.3 - 22:14:00.1 0542+498 [1, 2] 22:28:33.3 - 22:30:19.9 VENUS [1, 2] 22:53:33.3 - 22:54:19.9 0137+331 [1, 2] Fields: 13 ID Name Right Ascension Declination Epoch 1 0137+331 01:37:41.30 +33.09.35.13 J2000 2 0813+482 08:13:36.05 +48.13.02.26 J2000 3 0542+498 05:42:36.14 +49.51.07.23 J2000 4 0437+296 04:37:04.17 +29.40.15.14 J2000 5 VENUS 04:06:54.11 +22.30.35.91 J2000 6 0521+166 05:21:09.89 +16.38.22.05 J2000 7 1411+522 14:11:20.65 +52.12.09.14 J2000 8 1331+305 13:31:08.29 +30.30.32.96 J2000 9 MARS 14:21:41.37 -12.21.49.45 J2000 10 NGC7027 21:07:01.59 +42.14.10.19 J2000 11 NEPTUNE 20:26:01.14 -18.54.54.21 J2000 12 URANUS 21:15:42.83 -16.35.05.59 J2000 13 JUPITER 00:55:34.04 +04.45.44.71 J2000 Data descriptions: 2 (2 spectral windows and 1 polarization setups) ID Ref.Freq #Chans Resolution TotalBW Correlations 1 4885.1 MHz 1 50000 kHz 50000 kHz RR RL LR LL 2 4835.1 MHz 1 50000 kHz 50000 kHz RR RL LR LL The sources of interest are JUPITER (13), 0137+331 (1), and 0521+166 (6). The times of interest are therefore 16-Apr-1999/14:00:46.7 to 21:59:10.0. myms.done() The DONE button in the lower right of the GUI gets rid of the tool (and the window) and should be used when you are absolutely done with that tool. Otherwise, use DISMISS or the close function if you want that tool later. You can see a list of the tools active using the Tools in use button on the Tool manager window. ------------------------------------------------------------------------------- 2. Examine and edit the data. The data may be plotted in various ways using the msplot tool. Interactive editing is also available in msplot. The flagger tool is available for basic data editing. >Follow 1.6 Data examination and inspection >>>BUG: The cookbook says "It can also be launched from the Toolmanager via the >>>Synthesis package, and msplot module." but it seems to be under General. >>>What gives here? [George says its on his list...] Create: general->ms->msplot Construct: msplot->mymsplot( edit=T, msfile='planets.ms') You can grab the msfile name from the Catalog GUI by using the "wrench" menu to the right of the msfile data entry window. Click on the file name in the window and then punch Send&dismiss to send that filename back to the msplot Constructor window. Also - be sure to set Edit to True! >>>BUG: I first did it without giving an msfile name and I had to wait about >>>30 seconds for it to come up with an empty msplot. It should be able to >>>detect the fact that you gave it no ms right away and not waste time! >>>I know of no use for bringing up msplot with no ms... >>>[submitted as defect AOCso03683: syntax/input checking needed for msfile] >>>NOTE(George): msplot can be started up without an input file, and you >>>can use the File menu to select an ms later. >>>BUG: The cookbook description of msplot is very thin. The link to the >>>msplot part helps (though the stable version of the doc points to the >>>wrong place. Note that you have to use the NEXT/PREVIOUS buttons at >>>the top to navigate the pages as there are no other links or index >>>to the help topics on this (I think there should be). >>>[not submitted, suggestion to/from NAUG?] You will first want to set the source names and polarizations. First look at the calibrator 0137+331 = FieldID 1 in RR,LL. Use the Data Selection button. >>>BUG: I tried the From MS in the wrench for Field ID and got the >>>Field Chooser which looked promising. However, it didnt seem to work >>>right and only put Id:1 in the window. How is it supposed to work? >>>[submitted as part of defect AOCso03684: various problems] >>>NOTE: after playing much later on, I see it inherited field 1 selection >>>from somwhere, but if you hit Refresh the plot comes up ok. First choose "Plot UV Coverage" and then Go. This shows the Fourier plane sampling of the uv-data. Baselines are in meters. >>>BUG: On UV Coverage wouldn't it be more useful to have the option to plot in >>>wavelengths (Klam, Mlam)? There might be a way to change this but I >>>don't see it. >>>[submitted as part of defect AOCso03684: various problems] Lets look at the configuration of antennas in the array. Choose 'Go' to the Data Selection, and select 'from MS' in the wrench to the right of the antennas field. Note the antennas in the center of the array for later use - we will pick one of these (5,7,11,23,24) as a reference. After done looking, Dismiss (we dont really want to select anything). >>>BUG: A plot of the antenna positions should be available from the >>>ms tool, this is so useful. >>>[not submitted, for NAUG consideration] Go back to the msplot tool window and select Plot X vs Y. Do not dismiss the plot window. The default X Axis is uv distance and Y Axis is amplitude. These are ok. GO. This is essentially an unresolved source and should have a flat amplitude versus uv distance. You see that there are some visibilities below the main group, particularly some near zero. The data is uncalibrated at this time but the VLA is normally set up so that the gains are about the same for all antennas (at about 1 Jy = 10). There is something wrong with this data. You can use the mouse cursor to draw a box around bad data and then use the Locate buttion to figure out what it is - locate says this is antenna 9. Leave msplot using the red Done button. >1.7.0.1 Direct editing We will use flagger to blow away bad antenna 9 Create: synthesis->flagger->flagger Construct: flagger->myflagger('planets.ms') myflagger.setantennas(ants=9) myflagger.setids(spectralwindowid=[1,2],fieldid=[1,6,13]) To actually flag things with flagger, you have to choose one of the functions that actually does flagging, e.g. timerange, query, etc. myflagger.timerange(starttime="15-Apr-1999/23:00:00.0", endtime="16-Apr-1999/23:00:00.0"); Go on the timerange should go away and flag it. While we are here, lets remove the first 10s of scans (when the antennas are often not quite there yet) using the "quack" function myflagger.quack(scaninterval='30s', delta='10s') Done with flagger, hit Done >>>BUG: The first time I did this, editing in msplot then in flagger then >>>back to msplot (starting the tool anew each time), even though I flagged >>>AN9 in flagger some integrations appeared in msplot. I have not tried >>>to repeat this problem, but it was strange. >>>[submitted as part of defect AOCso03684: various problems] Return to msplot and Plot X vs Y for Field 1 (0137+331). Go back to the msplot window and select time as the X Axis. GO. This has the full time range of data. Lets restrict it. Bring up the Data selection and hit the Range button to the right of the Time selection. Edit the ranges to say 1999/04/16/14:00:0.0 1999/04/16/22:00:0.0 then GO. You can see the scans on 3C48 (0137+331) and though most of the bad visibilities (down near zero) are gone, there are still low ones. Select them with the cursor by setting boxes. When all the bad points are covered, the hit the green Flag button. Do the same with Field ID 6 (0521+166) Now select Jupiter Field ID 13. GO. What happened? Jupiter is resolved so there are a spread of visibilities. Note that the last scan is very low - select and Flag these. Set the X Axis back to uvdist and GO. You now clearly see the envelope of the good data (and the ringing in the uvplane - what could this be due to?). Use the cursor to select the bad data under the good data and Flag it. When finished editing and looking, hit the red Done button. A window with a stop sign will come up and ask if you want to keep these edits - hit Yes. ------------------------------------------------------------------------------- 3. Set initial calibration models. By default, all sources have unit flux density point source models. Use the imager tool's setjy function to set (point-source) flux densities for amplitude calibrators. Before calibrating, make an uncalibrated image. Create: synthesis->imager->imager Construct: calibrater->myimager(filename='planets.ms') First thing, choose 'setup' from the Function group and click on setdata. Pick Jupiter (FieldID=13) and first IF. Note that some tools like imager have the list of functions broken into Function groups. This can often be confusing, as you have to guess which group functions are in and just what ones are necessary. It is easier to select the Fuction group 'all' to see all the available functions (though the list will be long!). I will assume you do so. myimager.setdata(fieldid=13, spwid=1); Need to get image parameters. Grab the 'advise' function and ask for a 30' field-of-view myimager.advise(fieldofview='30arcmin'); The output in logger says: Maximum uv distance = 16747.8 wavelengths Recommended cell size < 6.15795 arcsec Recommended number of pixels = 300 Dispersion in uv, w distance = 6146.7, 3491.4 wavelengths Best fitting plane is w = 0.272991 * u + -0.507509 * v Dispersion in fitted w = 2703.57 wavelengths Wide field cleaning is not necessary Choose 3arcsec pixels 512 square. Go back to setimage under the 'setup' functions myimager.setimage(nx=512, ny=512, cellx='3arcsec', celly='3arcsec', stokes='IQUV', fieldid=13, spwid=1) Set the weighting, pick weight function, choose 'natural' first myimager.weight(type='natural') >>>BUG: I first tried making the image without setting the weight as >>>in the script, but I got the error that the weights were zero! >>>[submitted as part of defect AOCso03685: weight and setjy problems] Now make the image, select function makeimage . We are not done with imager, so do not hit Done yet. (If you do you can make another anyway.) myimager.makeimage(type='observed', image='jupiter_uncal.im') Now bring up the default viewer. In the tool manager, click 'Tools in use' and select the 'dv viewer' then hit green Show button. Windows for the Data Manager and viewer should appear. If you cannot find an entry for dv in the list of tools in use, then choose the display->viewer->viewer construction and make on called dv. >>>BUG: maybe you should have a function in imager to bring up a viewer like >>>in the image package. On the other hand, maybe the viewer should be >>>separate. In any event, it would be nice if this was consistent across >>>the package. >>>[not submitted, for NAUG consideration] >>>NOTE: The following was pointed out by Kumar: A more subtle way to bring up the viewer is to click on the wrench in the tool entry for the image name, and choose the View option. This will bring up the viewer. You can also bring up the default viewer dv by clicking on the image name in the Catalog (dc), and hitting View. For a description of the viewer, see Volume 1 of Getting Results, Basic Tools, Chapter 3 on Display of Data. In the Data Manager, select the 'jupiter_uncal.im' File and hit Raster Image on the right. Look at the image in the viewer. Its pretty crappy, there is a blobby thing in the center. The problem is that there are significant phase errors in the data due to antenna gain errors. You can play with the viewer, for example, use the left mouse button which is assigned to Zoom to draw a box around the source in the center and double-click on it to zoom in. Use Unzoom button at bottom to unzoom, natch. Use right mouse assigned to rectangle drawing to draw a box in the lower right of image containing only noise. Go to Tools menu at top and select Image Analysis, which will bring up yet another window, hit the Statistics pulldown. Double-click with right-mouse button in the box you drew and the statistics should get filled in. Look at the rms, 3.128e-3 and compare it to the max in the image (how would you get that?) to see our crappy dynamic range right now! I get max=2.069e-2, use glish to compute '2.069e-2/3.128e-3' giving 6.61445013. Thus we have max=2.069e-2 rms (corner) = 3.128e-3 => dyn range=6.61 When done, do not bother to hit Done in the viewer, though you can Dismiss the Image Analysis window. Go back to the Mangager in the Tool Manager window. Time to calibrate... >1.8.1 Setting the calibration source model >1.8.2.1 continuum polarimetry case Our phase calibrator 0137+331 is also a recognized amplitude calibrator (3C48). We need to fill in the model column of the measurement set. In Tool manager, click on the 'imager' button which should be lurking after the Manager and Tools in use buttons. Grab 'setjy' from the 'calibrate' function group. myimager.setjy(fieldid=1,fluxdensity=[-1.0]) it will report in logger: Fourier transforming: replacing MODEL_DATA column 0137+331 spwid= 1 [I=5.489, Q=0, U=0, V=0] Jy, (Perley-Taylor 95) 0137+331 spwid= 2 [I=5.542, Q=0, U=0, V=0] Jy, (Perley-Taylor 95) >>>BUG: I thought we fixed this to say Perley-Taylor 99 long ago! >>>[submitted as part of defect AOCso03685: weight and setjy problems] >>>The source 0521+166 is also a known calibrator, but we will treat it as an >>>unknown target to see how good our calibration is. For this case, we >>>set its flux density to be 1 Jy and will solve for gains relative to this. myimager.setjy(fieldid=6,fluxdensity=[1.0,0.0,0.0,0.0]) Fourier transforming: replacing MODEL_DATA column 0521+166 spwid= 1 [I=1, Q=0, U=0, V=0] Jy, (user-specified) 0521+166 spwid= 2 [I=1, Q=0, U=0, V=0] Jy, (user-specified) ------------------------------------------------------------------------------- 4. Obtain standard "electronic gain" (G) calibration for calibrators. If planning to do instrumental polarization (D) calibration, application of parallactic angle (P) calibration is required. >1.8.2 Solving for complex gain >1.8.2.1 continuum polarimetry case Now we need to solve for gains Create: synthesis->calibrater->calibrater Construct: calibrater->mycalibrater(filename='planets.ms') Choose our two calibrators mycalibrater.setdata(msselect='FIELD_ID IN [1,6]') >>>BUG: this should be selected using a fieldid field (and spwid etc.) >>>[not submitted, Im pretty sure this is on the to-do list] Because we will be doing polarization calibration, we will need to rotate the RL and LR phases with the parallactic angle (explain this to students) in order to line them up for the leakage solution later on. (George's script says to do this every 5 seconds - this is probably overkill, and the integs are 3.3s so I will use 10sec instead.) mycalibrater.setapply(type='P', t=10.0) Now we want to solve for the G gain calibration table which we will output as 'planets.gcal', we do 1 minute solutions and use antenna 5 as the reference (e.g. its phase = 0) - this was one of the antennas at the center of the array. mycalibrater.setsolve(type='G', table='planets.gcal', t=60.0, refant=5) It is useful to see what we have set, use the state function which will print to the logger window mycalibrater.state() If we are happy with these, go ahead and solve mycalibrater.solve() ------------------------------------------------------------------------------- 5. Obtain bandpass calibration (B) (spectral-line observations only) for bandpass calibrator(s). This operation is very much like the G calibration. The G calibration obtained above (and P, if necessary) is applied on-the-fly to ensure that the B calibration is normalized. >1.8.3 Bandpass calibration This part is not needed here, as it is for spectral line channels only The continuum IFs (spectral windows spwid=1,2) are treated differently than channels. ------------------------------------------------------------------------------- 6. Transfer flux density scale from amplitude calibrator(s). Only the G calibration solutions obtained for the amplitude calibrators are properly scaled in amplitude. The solutions for the other calibrators are rescaled with the calibrator fluxscale function by demanding that the mean gain amplitude be the same for all calibrators. >1.8.4 Establishing the flux density scale >1.8.4.1 continuum polarimetry case Now lets go back and bootstrap 0521+166 from 0137+331 mycalibrater.fluxscale(tablein='planets.gcal',tableout='planets.gcal1', reference='0137+331', transfer='0521+166') Flux density for 0521+166 (spw=1) is: 3.732 +/- 0.001 Jy Flux density for 0521+166 (spw=2) is: 3.768 +/- 0.001 Jy The Perley-Taylor 99 fluxes are 3.774, 3.802 (you can run setjy on these with fluxdensity -1), so our calibration is about 1% low which is pretty darn close! Done for day. =============================================================================== Date - 2002.05.29-30 Version - Weekly (1.7 build 388) Directory - /home/sandrock/smyers/Testing/2002-05-23/ Data - /aips++/data/nrao/VLA/planets_6cm.fits =============================================================================== >>>NOTE: will have to rebuild imager and calibrater tools on restart, >>>but I will pretend this is one long session Plot up the solutions mycalibrater.plotcal(tablename='planets.gcal1',plottype='amp') You can see all the amplitudes of the gains are around the nominal value of 0.32-0.33 ~ sqrt(10), so g_i*g_j = 0.1 and thus 10 counts = 1 Jy. If you want you can plot the phases, selecting one antenna at a time in one spwid. What is the average change in phase between scans? How will these residuals after calibration affect the data? >1.8.5 Correct the data >1.8.5.1 continuum polarimetry case Apply the calibration tables we just made. Reset the setapply and setsolve, then setapply to the P and G, mycalibrater.reset() mycalibrater.setapply(type='P', t=10.0) mycalibrater.setapply(type='G', t=0.0, table='planets.gcal1') Be sure to set the data to include Jupiter=13 mycalibrater.setdata(msselect='FIELD_ID IN [1,6,13]') mycalibrater.state() mycalibrater.correct() >1.8.6 Calibration Progress We can now make another image to see how this first calibration went. >>>NOTE: I set these to what they were before quitting aips++ and restarting. >>>These will be already set when students are continuing session. In any >>>event, the settings used last time will still be there, but still have >>>to do GO on all these set functions myimager.setdata(fieldid=13, spwid=1); myimager.setimage(nx=512, ny=512, cellx='3arcsec', celly='3arcsec', stokes='IQUV', fieldid=13, spwid=1) myimager.weight(type='natural') myimager.makeimage(type='corrected', image='jupiter_cal1.im') Bring up the default viewer again, and plot up this image. You clearly see Jupiter now in the center, along with some "rays" extending out from the source - what are these? You can zoom in and see it is clearly not a point source, but still no real features are visible You can use the Data Manager to bring in the jupiter_uncal.im as a comparison. It will overlay the cal1.im - you have to go to the Tools->Register menu and deselect uncal.im to see the cal1.im What happened to the noise level and dynamic range? I get: max=1.548 rms (corner) = 4.944e-2 => dyn range=31.3 >>>BUG: As we have pointed out numerous times, there should be a way to >>>switch and blink between registered images like between poln planes >>>[not submitted - save for viewer meeting] You can also play with the imager weighting scheme for example, go back and select uniform weighting myimager.weight(type='uniform') myimager.makeimage(type='corrected', image='jupiter_cal1unif.im') Bring this into the viewer as a Raster, you will have to hit the Update button at lower left in the Data Manager to see the file in the list. What was the effect of changing the weighting, and why? I get: max=1.010 rms (corner) = 2.987e-2 => dyn range=33.81 the rms is better than the natural, probably due to reduced sidelobes What is the synthesized beam? we can make this in imager myimager.weight(type='natural') myimager.makeimage(type='psf', image='jupiter.psf') Look at it in viewer >>>BUG: it would be great to have the imagefitter available from the viewer >>>directly! [save for viewer discussion] Use the fitpsf function of imager to get the beamsize myimager.fitpsf(psf='jupiter.psf') Beam fit: 22.2519 by 15.9719 (arcsec) at pa -46.7583 (deg) You can contrast this with myimager.weight(type='uniform') myimager.makeimage(type='psf', image='jupiter_unif.psf') myimager.fitpsf(psf='jupiter_unif.psf') Beam fit: 17.7531 by 13.2426 (arcsec) at pa 135.048 (deg) We will go with uniform weighting for the better resolution, but you can switch this at any time during the imaging operations. >>>BUG: the cookbook does not mention the psf image. In general, its flow >>>should contain more suggestions for data exploration - it is very spare >>>in its current state [not submitted - save for cookbook upgrade with >>>Debra] We jump ahead and make a clean image of I right now you can look ahead in the cookbook at section 1.10 Imaging. To see the source structure, we need to remove the smearing caused by the psf (due to gaps in the uv coverage) We will use the interactive feature of clean to mask off regions See the User Reference manual entry for clean for explanation of the inputs myimager.clean(niter=1000, gain=0.1, model='jupiter.mod1', image='jupiter.cln1', residual='jupiter.resid1', mask='jupiter.mask1', interactive=T, npercycle=500) A display of the image pops up inviting you to define a mask region. First, hit the "Add region" green button, use the rectangle tool (right mouse) to define a small box around the source (zoom in using Zoom tool = left mouse if desired), and double-click (with right mouse) in it to define the region. Then hit "Done with masking" to start the cleaning. It will probably come back with a residual image that has the source cleaned from it, just hit "Done with masking" to finish the clean. >>>BUG: I had some problem with this at first, as it wasn't clear what you >>>had to do to define the region (e.g. hit Add, define box with rectangle >>>using right mouse, double click in it). It would be a great help if there >>>was some short statement sent to the logger about what you were being asked >>>to do, like happens in AIPS TVALL and TVFLG for example which says what >>>you are supposed to do with the buttons. This would be simple to add and >>>I think would help users alot. And it would apply to more than this >>>function, perhaps lots of viewer based things. You could also do it by >>>having a window at the bottom of the viewer for special dialog messages >>>like this. [not submitted - for NAUG/viewer discussion] Use the viewer to look at the restored (clean model + residuals) image jupiter.cln1: notice the much lower noise level - most of the noise was sidelobes from the source! Zoom in on Jupiter - its still a fuzzy blob, alas. You can compare the cleaned and uncleaned jupiter_cal1unif.im if you want. Image statistics: max=1.042 rms (corner) = 9.289e-3 => dynrange=112.2 ------------------------------------------------------------------------------- 7. Obtain instrumental polarization calibration (D) (polarization observations only). Using a full-polarization model for a calibrator (obtained by imaging or other means, and set with the imager tool's setjy function), the D calibration procedes in much the same way as the previous calibrations, with these previous calibration factors applied on-the-fly. >1.9 Instrumental Polarization Calibration >1.9.1 continuum polarimetry case (only) Make IQUV images of the polarization calibrator 0137+331 = field 1: myimager.setdata(fieldid=1, spwid=1); myimager.setimage(nx=512, ny=512, cellx='3arcsec', celly='3arcsec', stokes='IQUV', fieldid=1, spwid=1) >>>QUERY: why is there spwid and fieldid in both setdata and setimage??? >>>NOTE: From Kumar, the selection of image parameters in setimage can in >>>general be different than in the data selection, e.g. in the case of >>>multiple fields in a mosaic where you want an image center etc. myimager.weight(type='natural') myimager.makeimage(type='corrected', image='3c48_cal1.im') In the viewer, bring up 3c48_cal1.im and zoom in around the maximum in the I image. Pull up the Statistics panel from the Image Ananlysis Tool, then use the rectangle tool to box around the source, and double-click to fill in the information. Note the maximum. Step through the planes Q,U,V image planes and double click in the box, noting the maximum or minimum, whichever has the largest absolute value. The images will appear black, due to the color scale set to the intensity. Record these values. I get: I=5.491 Q=-0.115 U=-0.248 V=0.006 You can play with Adjusting the data range to look at Q,U,V if you wish, particularly useful is the little histogram widget next to the Data Range window in the Adjust panel Make the spwid=2 image and do the same. I get: I=5.545 Q=-0.220 U=0.186 V=0.006 >>>NOTE: There is a better way to get the I,Q,U,V at the center >>>of the image - average up the visibilities with zero phase rotation >>>(why does this work?). Will describe this if it is ready. >>>BUG: George's jupiter.g script says to use the glish function ptsrc.g >>>e.g. include 'ptsrc.g', but I cannot find that anywhere (aips++ does >>>not find it). George's values were : >>> >>>myimager.setjy(fieldid=1,spwid=1,fluxdensity=[5.408,-0.114,-0.243,0.0]); >>>myimager.setjy(fieldid=1,spwid=2,fluxdensity=[5.461,-0.216,0.183,0.0]); >>> >>> which agree with the numbers I got from the image. Use our values for Q and U, but set the I to the values we got from the Perley-Taylor 99 run of setjy above, and set V=0 rather than using the noisy value off the images: myimager.setjy(fieldid=1,spwid=1,fluxdensity=[5.489,-0.115,-0.248,0.0]); myimager.setjy(fieldid=1,spwid=2,fluxdensity=[5.542,-0.220,0.186,0.0]); Now go to the mycalibrater tool, and reset the setapply etc. for safety if desired (you can go with the old ones, you can check with the .state() function) mycalibrater.reset() mycalibrater.setdata(msselect='FIELD_ID==1') mycalibrater.setapply(type='P', t=10.0) mycalibrater.setapply(type='G', t=0.0, table='planets.gcal1') Set it up to solve for the D-terms with a single solution over the 1.5 days of data (really 1 day, but put a buffer). I will try a 30sec pre-avg mycalibrater.setsolve(type='D', t=1.5*86400.0, refant=5, preavg=30.0, table='planets.dcal'); mycalibrater.state() mycalibrater.solve() We can look at the dterm solutions mycalibrater.plotcal(plottype='DRI', tablename='planets.dcal') which plots the imaginary vs. real. All are less than 1% except one antenna with about a 4% leakage (AN15, IF2) >>>BUG: this should be in the cookbook also ------------------------------------------------------------------------------- 8. Apply all calibrations to target source data. The calibrater's correct function applies calibration to the observed data and writes the result in its own column in the dataset. >1.9.1 continuum polarimetry case (only) - continued Now we go ahead and apply all calibrations determined to this point: mycalibrater.setdata(msselect='FIELD_ID IN [1,6,13]' mycalibrater.reset() mycalibrater.setapply(type='P', t=10.0) mycalibrater.setapply(type='G', t=0.0, table='planets.gcal1') mycalibrater.setapply(type='D', t=0.0, table='planets.dcal') mycalibrater.state() mycalibrater.correct() ------------------------------------------------------------------------------- 9. Image the target source. The imager tool is used to form images from the calibrated data, and to deconvolve these images. >1.10 Imaging >1.10.0.1 continuum polarimetry case Now we make our images (and clean) again myimager.setdata(fieldid=13, spwid=1); myimager.setimage(nx=512, ny=512, cellx='3arcsec', celly='3arcsec', stokes='IQUV', fieldid=13, spwid=1) myimager.weight(type='uniform') myimager.clean(niter=1000, gain=0.1, model='jupiter.mod2', image='jupiter.cln2', residual='jupiter.resid2', mask='jupiter.mask2', interactive=T, npercycle=500) >>>BUG: The cookbook should be updated to reflect the fact that clean >>>now restores the full image automatically (I think). >>>[Im sure George knows this] Look at jupiter.cln2 in viewer. The I pol should be the same as before, but now Q and U should be better. Look at pols 2 and 3 (Q and U) - use the Adjust->Data range using the histogram widget. Notice that the images are sensible, though strange looking! Those are real polarized features around Jupiter! >>>BUG: the Data range and especially the histogram widget seem to always use >>>the I pol plane, and do not change to reflect the plane actually being >>>viewed. There should be a way to adjust each plane, by choosing which >>>plane the histogram is of. Also, there seems to be no button or control >>>to reset the contrast and colormap shift/slope (set with cursor) to >>>the original default, other than to click in the center. I would think >>>some sort of global reset (like the Unzoom button at the bottom) would >>>be useful. [submitted as defect AOCso03697] ------------------------------------------------------------------------------- 10. As necessary, self-calibrate the target source data and reimage. New and improved calibration may be obtained by repeating the initial calibration step(s), this time using the model image obtained in the first round of imaging. Re-image. >>>BUG: This also should be in the cookbook by the time of the school!!!! We will use our clean model to improve the calibration of the Jupiter data itself. This is called "self-calibration". First, we need to fill the MODEL column of the measurement set for Jupiter. We will only do spwid=1 at this point, all the following can be repeated for spwid=2 if desired. Grab the myimager tool. myimager.setdata(fieldid=13,spwid=1); # Jupiter, Spw 1 myimager.ft(model='jupiter.mod2') # Set model to current model Grab the mycalibrater tool. Set it up to apply the P and D and previous G corrections. mycalibrater.setdata(msselect='FIELD_ID==13 && SPECTRAL_WINDOW_ID==1') mycalibrater.reset() mycalibrater.setapply(type='P', t=10.0) mycalibrater.setapply(type='G', t=0.0, table='planets.gcal1) mycalibrater.setapply(type='D', t=0.0, table='planets.dcal) We want an incremental solution over the current G, so we will to an "atmospheric phase" T calibration (phase-only) mycalibrater.setsolve(type='T', table='jupiter.tcal1', t=60.0, phaseonly=T, refant=5) >>>NOTE: The Jupiter scans are 100 sec long. What we really want is scan-based >>>solution intervals. As it is we are getting failed solutions, probably >>>due to mismatch in the solution intervals and the scans. mycalibrater.state() mycalibrater.solve() Use mycalibrater.plotcal to look at the solutions. Note that there are some >10deg corrections. Now apply the solutions P,G,D and now T mycalibrater.setapply(type='T', t=0.0, table='planets.tcal1) mycalibrater.state() mycalibrater.correct() >>>BUG: The first time I did this, I gave it the T the wrong name of >>>planets.tcal1 in setapply. There was no warning or error at any point, >>>I only noticed that it didnt report a T application. There should be >>>some sort of checking here - it is way too easy to screw up without >>>noticing. [submitted as defect AOCso03698] Now we make our images and clean once more, so grab imager:myimager myimager.setdata(fieldid=13, spwid=1); myimager.setimage(nx=512, ny=512, cellx='3arcsec', celly='3arcsec', stokes='IQUV', fieldid=13, spwid=1) myimager.weight(type='uniform') myimager.clean(niter=1000, gain=0.1, model='jupiter.mod3', image='jupiter.cln3', residual='jupiter.resid3', mask='jupiter.mask3', interactive=T, npercycle=500) Look at jupiter.cln3 in viewer. Get stats: max = 1.066 rms (corner) = 6.689e-3 => dynrange = 159.4 Getting better! You can see the 'wings' in the I image now. Q and U are still looking good, though its hard to tell from Q and U what is going on. ------------------------------------------------------------------------------- 11. Repeat the self-calibration as necessary. Self-calibration need not be not limited to only the G calibration factor; any solvable calibration component (G, B, and/or D, in the present example) may be self-calibrated. We will skip extra self-cal iterations in this demo >1.12 Image analysis >>>BUG: This section of the cookbook should be expanded, preferably by the >>>time of the school. Note that what is there suggests using the image >>>tool, which though it has lots of functionality, it is more a toolkit >>>and it is not clear that this is the best choice for beginners. As for >>>the image.view() how is this different from just using the viewer directly? >>>The separate tools imagefitter, imagepol, and imageprofilfitter which are >>>built on image seem to be much more appropriate for the cookbook. >>>[just a rant] Now we will make specialized polarization vector images This is described in Volume 1 of Getting Results, Basic Tools, Chapter 4 on Image Analysis: >4.14 Polarimetry Create: general->images->imagepol Construct: imagepol->myimagepol('jupiter.cln3') Form a linear polarization intensity image, e.g. P = sqrt( Q^2 + U^2 ) myimagepol.linpolint(return='mylinpolint')->mylinpolint When you do 'Go' on this what it will actually do is go an make an image tool which will appear in the Tool Manager - this is very confusing at first, since this is not the normal behavior we have seen so far when you do Go on functions. Bear with this... Lets store this on disk using the subimage mylinpolint.subimage(outfile='jupiter.pint')->mysubim This makes yet another image tool, mysubim, but more importantly writes a file. Use the viewer to look at jupiter.pint - why does it look so strange in the noise-dominated regions? Now, hit the red 'Done' button on the image:mysubim tool. Also you are Done with the image:mylinpolint tool. >>>BUG: this behavior is annoying - you should be able, for example, >>>to the return in order to only make the file without the >>>extraneous tool. If you do the return, it does nothing. >>>There should be a way to not make these spurious image tools. >>>In Georges script, my guess is that he needs to do lots of cleanup >>>at the end that is not currently there. [not submitted, for discussion] >>>NOTE: Georges script makes a linpolposang image, which is not really >>>used for anything, so I will skip that. Now we want to make a vector image of the electric vectors of the polarization. Grab the imagepol:myimagepol tool from the manager window and choose the complexlinpol function and set output file to 'jupiter.cmplx' myimagepol.complexlinpol(outfile='jupiter.cmplx') which will form an image of Q + iU. >>>NOTE: Georges script made a complexfraclinpol image which I had trouble >>>dealing with. Instead, I follow the Getting Results and make a >>>complexlinpol image. >>>BUG: note that this tool does what is expected - only creates a file >>>and not an image tool! The others should do this, or at least have an >>>option to. Done with imagepol:myimagepol Now grab the viewer. Load the raster image of jupiter.pint you had previously. Load jupiter.cln3 as a raster image. Load jupiter.pint as a contour plot. Load jupiter.cmplx as a vector map. Zoom in on Jupiter in the center. We need to fiddle now with these plots which are all overlaid. Since we are going to be doing some heavy viewer gymnastics, you might want to look at the Volume 1 of Getting Results, Basic Tools, Chapter 3 on Display of Data which describes the Viewer in detail. Use the DisplayData->Adjust for jupiter.cmplx from top menu. Play with the Amplitude scale factor, x-increment, y-increment to get nice vector sizes and spacings. I found incs of 2 and scale of 2 ok. It would be nice to get rid of the vectors where there was low polarization or intensity signal. Grab the Tool->Image Analysis, select jupiter.pint, find the rms in a corner of the displayed region (I get 9.5e-4), so a reasonable mask is 4 sigma or 0.0038 Jy. In the Adjust panel for the vector image, set the Mask expression to: jupiter.cln3 > 0.0038 note - you cannot mask on the jupiter.cln3 image since it has a different shape (4 planes) to the cmplx image! We are not quite done yet. The electric-vector position angles (EVPA), which are defined as EVPA = 0.5*arctan(U/Q) are not fully calibrated. The ratio of U/Q that sets these is arbitrary as we have not specified that true ratio for our calibrator 0137+331. In our calibration above, we found Q=-0.115 U=-0.248 in spwid=1 and thus the R-L phase difference is arctan(-0.248/-0.115); in the glish window, type atan(-0.248/-0.115) giving 1.1366 or converting to degrees atan(-0.248/-0.115)*180/pi gives 65.12 degress. Since U and Q are both negative, this is the third quadrant, and thus the angle is really -180 deg from this (angles run from -180 to 180) or -114.88 deg. Looking on the web at the VLA calibration pages, we see that the true R-L phase difference for 3C48 at C-band is -148 deg. Thus, our R-L phases which are arg(Q+iU) should be rotated by adding -148 - (-115) = -33 deg to the values given. You can do this in the Adjust for the vector image, by entering a value for the Extra rotation. Since the Phase Type is set to polarimetric, the vectors shown have the correct factor of 0.5 to get EVPA from R-L phase. Thus we should add -33/2 = -26.5 deg, so put that in the Extra rotation (e.g. put in -33/2) and Apply. To see where the B-vectors like, add in an additonal 90 deg, e.g. set the Extra rotation to -33/2 + 90. >>>BUG: first time through I goofed and used arctan(Q/U), doh! This is because >>>the documentation on imagepol.linpolposang said "This function (short-hand >>>name lppa) returns the linearly polarized position angle image ( >>>0.5tan-1(Q/U)) in degrees", which is incorrect. >>>[submitted as defect AOCso03708] >>>NOTE: the cookbook should really link to the explanations like http://aips2.nrao.edu/weekly/docs/user/Display/node159.html#viewer:vddoptions >>>for all the adjust options Now we want a hardcopy of this. You can find a description (sort of) of this in the Volume 1 of Getting Results, Basic Tools, Chapter 3 on Display of Data, which describes the viewer, in the section: >3.12 A step-by-step recipe to make a post-script plot of your image First, bring up a DisplayData->Adust window for jupiter.cln3 (raster). Put the axis labels on, by using the Axis labels panel of Adjust, and setting 'Axis labelling & annotation" to True. You can add a title like 'Jupiter 6cm'. You can get a representation of the synthesized beam by using the Beam panel in the Adjust window, and setting 'Plot beam' to True. In the Basic function of Adjust, set the Colormap to Greyscale 2 to reverse the image from white on black to black on white versus Greyscale 1. Next, get (or bring up) the Adjust window for jupiter.cmplx (vector) set the Line color (at bottom of the basic functions) to black or blue. Finally, bring up a DisplayData->Adust window for jupiter.pint (contours) and set the Line color to black or red for contrast. Hit the Print... button at the bottom of the viewer display. Select the output file 'jupiter.ps', select the Output media 'A4' (why the heck is the default A4? Bloody Aussies...). Hit the 'Save PS' button. You can also print the image if you want, using the Print button (choose the right printer nearest your office, should be what the workstation default is). >>>BUG: The default plot comes out with black vectors and contours, even >>>though they are white on the display. It would be nice to have a quick >>>way to reset foreground, background and colormap to printer-friendly >>>values, perhaps though a button on the Canvas Print Manger window. >>>[not submitted, for discussion, and besides I think this was harped >>>on a while ago] This is a good place to stop the main tutorial. ------------------------------------------------------------------------------- Summary: After alot of pain, I got through this OK with a sort of Ok polarization image (modulo the angle problem). With some tweaking, I think the VLA cookbook chapter can be made to read smoothly enough to guide the students (at least for this dataset, I don't know about the spectral line one). In practice, I might skip some of the exploratory bits like weighting differences. This should work...