Date - 2000.04.14 to 2000.04.18 Tester - S.T. Myers (AOC) Platform - Linux (kernow) Version - Stable (1.4 build 1) Note: bugs are denoted by lines with prefix: >>>BUG and queries are denoted by lines with prefix: >>>QUERY while suggestions for improvement are prefixed by >>>ENHANCEMENT ------------------------------------------------------------------------------- GOAL: Read CBI data into aips++ and try to do some basic things 1. Read data into measurement set include 'synthesis.g' dowait:=T m:=fitstoms(msfile='J0739+016.ms',fitsfile='J0739+016.uvf'); m.close() Had trouble with some of the FITS keywords, but seems to have gone ahead anyway. 2. Get UV range plot of data Try GUI first. Could not figure out how to make it plot all the data on one plot. It seemed to be iterating over something. >>>ENHANCEMENT: VISPLOT AS IT STANDS IS BARELY USEABLE. DEVELOPMENT OF A MORE COMPREHENSIVE AND INTUITIVE VISIBILITY VISUALIZATION TOOL WOULD BE PREFERABLE 3. Now try imager. Use script imgr:=imager('J0739+016.ms') # make imager object ERROR: unhandled mount type (eventually ended with 1711 errors!) I knew this was going to be a problem, as Tim Pearson put the mount type in the FITS header as imgr.summary() =============================================================================== MeasurementSet Name: J0739+016.ms MS Version 1 =============================================================================== Observer: Project: Array: CBI Fields: 1 ID Name Right Ascension Declination Type 1 J0739+016 07:39:18.03 +01.37.04.62 J2000 Data records: 210100 Total integration time = 1089.25 seconds Observed from 01-Mar-2000/23:58:22 to 02-Mar-2000/00:16:31 Spectral Windows: 10 Ref.Freq RestFreq #Chans Resolution TotalBW Correlations 35500 MHz continuum 1 1e+06 kHz 1e+06 kHz LL 34500 MHz continuum 1 1e+06 kHz 1e+06 kHz LL 33500 MHz continuum 1 1e+06 kHz 1e+06 kHz LL 32500 MHz continuum 1 1e+06 kHz 1e+06 kHz LL 31500 MHz continuum 1 1e+06 kHz 1e+06 kHz LL 30500 MHz continuum 1 1e+06 kHz 1e+06 kHz LL 29500 MHz continuum 1 1e+06 kHz 1e+06 kHz LL 28500 MHz continuum 1 1e+06 kHz 1e+06 kHz LL 27500 MHz continuum 1 1e+06 kHz 1e+06 kHz LL 26500 MHz continuum 1 1e+06 kHz 1e+06 kHz LL Feeds = 12: printing first row only Antenna Spectral Window # Receptors Polarizations 1 -1 2 [ R, L] Antennas: 13 1-13: RX0 RX1 RX2 RX3 RX4 RX5 RX6 RX7 RX8 RX9 RX10 RX11 RX12 ------------------------------------------------------------------------------- Tables: MAIN 210100 rows ANTENNA 13 rows ARRAY 1 row FEED 12 rows FIELD 1 row OBSERVATION 1 row OBSLOG 25 rows SOURCE (see FIELD) SPECTRAL_WINDOW 10 rows SYSCAL WEATHER =============================================================================== General: MeasurementSet is J0739+016.ms Beam fit is not valid Image: parameters not yet set (use setimage) Data selection settings: (use setdata to change) All data selected Options settings: (use setoptions to change) Gridding cache has 4194304 complex pixels, in tiles of 16 pixels on a side Gridding convolution function is Spheroidal wave function No primary beam correction will be made Image plane padding : 1 Seems to have gone ahead and made the imager object anyway! Continue... There seem to be 10 spectral windows (the 10 CBI channels) each of 1 channel 1 GHz wide. imgr.setimage(cellx='1arcmin', celly='1arcmin', nx=128, ny=128, stokes='LL', spwid=[1,2,3,4,5,6,7,8,9,10], fieldid=1) ERROR: Illegal Stokes string LL Manual sez "Allowed: 'I'|'IV'|'IQU'|'IQUV' " WHAT, NO RR OR LL???? >>>ENHANCEMENT: OPTION FOR RR OR LL (OR SINGLE LINEAR) POLARIZATION SHOULD BE ALLOWED. I guess "I" will do what I want, however one of the 13 CBI receivers is cross-polarized (to R) versus the 12 others (in L) in order to do cmbr polarization. Have to watch that this gets treated right. imgr.setimage(cellx='1arcmin', celly='1arcmin', nx=128, ny=128, stokes='I', spwid=[1,2,3,4,5,6,7,8,9,10], fieldid=1) imgr.setdata(spwid=[1,2,3,4,5,6,7,8,9,10], mode="none") imgr.uvrange() imgr.weight('robust') Sum of weights = 168 What does this mean?? bmaj:=F; bmin:=F; bpa:=F; # Set up some return variables imgr.image(type='psf',image='J0739+016.psf') # Make the PSF imgr.fitpsf(psf='J0739+016.psf', bmaj=bmaj, bmin=bmin, bpa=bpa) # Fit beam Beam fit: 256.005 by 239.323 (arcsec) at pa 27.4288 (deg) imgr.image(type='observed',image='J0739+016.dirty') # Make the dirty image Make an image object im := image('J0739+016.dirty') View it im.view(raster=T) Wow! You can see the source (calibrator J0739+016). Use Statistics function Max = 1.363 Jy/beam on peak Rms = 67.5 mJy/beam off source (does include sidelobes) im.done() Now try to clean the image. Do inner 1/4 around center. imgr.make('J0739+016.clean') blc_max := 32 trc_max := 96 imgr.boxmask(mask='J0739+016.mask',blc=blc_max,trc=trc_max) IN HINDSIGHT, THESE BLC AND TRC ARE INCORRECT... imgr.clean(mask='J0739+016.mask',model='J0739+016.clean',niter=100,gain=0.03) Clean used 100 iterations to get to a max residual of 0.140504 imgr.clean(mask='J0739+016.mask',model='J0739+016.clean',niter=100,gain=0.03) Initial maximum residual: 1.01709 Clean used 100 iterations to get to a max residual of 0.107537 I have assumed that this will have restarted the clean correctly from where it left off... Look at the residual so far: imgr.residual(model='J0739+016.clean',image='J0739+016.resid') im := image('J0739+016.resid') im.view(raster=T) im.done() Peak = 1.058 mJy/bm rms outside peak = 55.4 mJy Something is not cleaning right - check current clean model im := image('J0739+016.clean') im.view(raster=T) im.done() This shows a +0.73 Jy pixel next to a -0.49 pixel below it. Clearly something has gone wrong. >>>NOTE: LATER ON I FIND OUT WHY THIS BOMBED - I SPECIFIED BLC,TRC INCORRECTLY TO THE MASK End of 04-14 session. 4. Pick up again 2000.04.17 with looking at the ms generated. Use GUI version of ms, equivalent to myms := ms('J0739+016.ms') We will try retrieving some data, similar to the example in the User's Reference myms.selectinit(arrayid=1,spectralwindowid=5) myms.select([antenna1=[1,2]) Now get a record with the desired stuff. Note it is unclear now just what you can get. In the URM an example is "uvw amplitude corrected_amplitude" so it is clearly more complex than just "corrected_data". Note that myms.summary() gives no help here. In the ms.range description we read "Possible items include most scalar columns, ifrnumber (1000*antenna1+antenna2), uvdist(ance), u, v, w, amplitude, phase, real and imaginary of the data (and corrected and model versions of these). See the table at the top of the document to find the exact list." Sure enough, there is a big table in the ms description :-) The GUI is slowing me down, so I will revert to the glish level... rec := myms.getdata("weight corrected_data model_data"); print rec.weight[1], rec.corrected_data[1], rec.model_data[1] 0.0998466462 3.45497704+6.03565454i 0.213887364+0.101477325i or rec := myms.getdata("uvdist corrected_amplitude corrected_phase"); print rec.uvdist[1], rec.corrected_amplitude[1], rec.corrected_phase[1] 4.5825758 6.95456648 1.05089724 Following the plotter trail to the Glish/PGPLOT bindings, we try f := frame (); pg := pgplot (f); We get an itty bitty frame, then a slightly bigger black background, myms.range("uvdist corrected_amplitude corrected_phase"); [uvdist=[0.999999977 5.00000013] , correctedamplitude=[0.118513152 38.8557816] , correctedphase=[-3.14056754 3.14126396] ] Note the uvdist must be in either meters or hundreds of wavelengths, so I would guess meters. I think Tim puts out the visibilities in Jy (after calibration was applied). pg->env ( 0.0, 6.0, 0.0, 40.0, 0, 0 ); Whoa. That really is a tiny plot. Try to drag it bigger. I guess I should redraw or refresh somehow. After reading the stuff under "Pgplot and Window Resizing" I decide to just do it over... pg->eras (); pg->env ( 0.0, 6.0, 0.0, 40.0, 0, 0 ); pg->pt ( rec.uvdist, rec.corrected_amplitude, -1 ); Good, that worked. This seems easier than trying to drive the visplot function from the GUI, and is scriptable. pg->end (); I wonder how you get rid of the frame. Does "frame" have a done function? Well, after reading under the Glish/Tk in the Glish manual about the frame, I dont see anything. Do it manually... myms.done(); 5. Browse the ms checking the tables under J0739+016.ms in the dc GUI. Will eventually need to get the FEED info for the peculiar CBI setup with a fixed deck with rotation FEED table has some questionable entries: INTERVAL(s) ... TIME(UTC) 1.79769313e+308 ... 1858/11/17/00:00:00 ANTENNA table has some strange defaults including the usual garbled text in the POSITION columns DISH_DIAMETER(m) MOUNT 25 bizarre I doubt if these get used for cleaning so I doubt this is a big problem. 6. Check the psf made earlier in imager im := image('J0739+016.psf') im.view(raster=T) im.done() Beam looks fine, about like the dirty image as expected. Why did clean fail before? 7. 2000.04.18 : Try clean again, make ms from scratch include 'synthesis.g' dowait:=T m:=fitstoms(msfile='J0739+016.ms',fitsfile='J0739+016.uvf'); m.close() imgr:=imager('J0739+016.ms') As before, gave a whole lot of errors "unhandled mount type". Noted strange behavior that it made 3 small SEVERE ERROR windows in addition to the one it was printing in, but the extras went away when the errors stopped being generated. Also gave errors (reported in the glish window) error, self.vsb is not an agent WheneverStmt ref count error, self.text is not an agent WheneverStmt ref count >>>BUG: IS THIS BEHAVIOR EXPECTED? imgr.setimage(cellx='1arcmin', celly='1arcmin', nx=128, ny=128, stokes='I', spwid=[1,2,3,4,5,6,7,8,9,10], fieldid=1) imgr.setdata(spwid=[1,2,3,4,5,6,7,8,9,10], mode="none") imgr.uvrange() imgr.weight('robust') bmaj:=F; bmin:=F; bpa:=F; imgr.image(type='psf',image='J0739+016.psf') imgr.fitpsf(psf='J0739+016.psf', bmaj=bmaj, bmin=bmin, bpa=bpa) Beam fit: 256.005 by 239.323 (arcsec) at pa 27.4288 (deg) imgr.image(type='observed',image='J0739+016.dirty') Frequency = 26.5 GHz, synthesized continuum bandwidth = 10 GHz im := image('J0739+016.dirty'); im.view(raster=T); im.done(); im := image('J0739+016.psf'); im.view(raster=T); im.done(); These look OK. imgr.make('J0739+016.clean'); blc_max := [ 32, 32, 0, 0 ]; trc_max := [ 96, 96, 0, 0 ]; This could be where I went wrong before, I just used blc=32 trc=96 imgr.boxmask(mask='J0739+016.mask',blc=blc_max,trc=trc_max); imgr.clean(mask='J0739+016.mask',model='J0739+016.clean',niter=300,gain=0.03); Clean used 300 iterations to get to a max residual of 0.06064 This seems much more sensible. im := image('J0739+016.clean'); im.view(raster=T); im.done() Model consists of one pixel of 0.92 Jy and one of 0.44 right below imgr.restore(model='J0739+016.clean', image='J0739+016.clean.restored', bmaj=bmaj, bmin=bmin, bpa=bpa) THIS GIVES ERROR : named arg "bmaj" does not match any formal in call to function (val model = , val complist = , val image = , val residual = , val async = !dowait) { { { { { (self.restoreRec.model := model); (self.restoreRec.complist := complist);} (self.restoreRec.image := image);} (self.restoreRec.residual := residual);} (returnval := defaultservers.run(self.agent, self.restoreRec, async));} return returnval} File: imager.g, Line 239 Stack: .() EVEN THOUGH - print bmaj,bmin,bpa [value=256.004608, unit=arcsec] [value=239.322708, unit=arcsec] [value=27.4288387, unit=deg] Lets try imgr.restore(model='J0739+016.clean', image='J0739+016.clean.restored', bmaj='247arcsec', bmin='247arcsec') >>>BUG: STILL GIVES SAME ERROR, EVEN THOUGH THE DOCUMENTATION SAYS TO DO IT THIS WAY. WHAT IS WRONG? Try imgr.restore(model='J0739+016.clean', image='J0739+016.clean.restored', usefit=T) SAME ERROR. Connect to imgr tool with GUI, and show the restore function. NOTE: THERE ARE NO INPUTS OTHER THAN model,complist,image,residual >>>BUG: FIX DOCUMENTATION Lets try using clean itself to make the restored image imgr.clean(mask='J0739+016.mask',model='J0739+016.clean', image='J0739+016.clean.restored',niter=0); Initial maximum residual: 1.24842 Clean used 0 iterations to get to a max residual of 1.24842 DID THIS MAKE THE PROPER RESTORED IMAGE? (Note: should I have specified fixed=T?) im := image('J0739+016.clean.restored'); im.view(raster=T); im.done() This looks more like a residual image. Lets try this from the top of clean again (remove previous files J0739+016.clean* ) imgr.make('J0739+016.clean'); imgr.clean(mask='J0739+016.mask',model='J0739+016.clean', image='J0739+016.clean.restored',residual='J0739+016.clean.residual', niter=300,gain=0.03); Clean used 300 iterations to get to a max residual of 0.06064 Good - as before. im := image('J0739+016.clean.restored'); im.view(raster=T); im.done() This worked. A nice point source of 1.355 Jy/bm max, with rms=30 mJy/bm in the cleaned region vs. 60 mJy/bm outside the cleaned region. 8. Now try a self-cal cycle on the data Make calibrater object and then in imager Fourier transform the model and insert into the MODEL_DATA column cl:=calibrater('J0739+016.ms') imgr.ft(model='J0739+016.clean') Set selfcal parameters for phase-only gain solution, then solve cl.setsolve (type="T",t=10.0,dophase=T,table='J0739+016.phasecal',append=F); : named arg "dophase" does not match any formal in call to function (val type = , val t = 60, val phaseonly = F, val table = , val append = F) { { { { { (self.setsolveRec.type := type); (self.setsolveRec.t := t);} (self.setsolveRec.phaseonly := phaseonly);} (self.setsolveRec.table := table);} (self.setsolveRec.append := append);} return defaultservers.run(self.agent, self.setsolveRec)} File: calibrater.g, Line 109 Stack: .() BOOM. LOOKS LIKE THERE HAVE BEEN SOME CHANGES NOT PROPAGATED THROUGH CORRECTLY EVEN THOUGH THIS IS THE "STABLE" VERSION. Fire up the Tool Manager for cl object. Aha, the dophase input has been changed to phaseonly. >>>BUG: THERE SHOULD BE A MORE ELEGANT WAY THAT A FUNCTION CAN LET YOU KNOW YOU GAVE AN UNKNOWN INPUT PARAMETER!!!! >>>ENHANCEMENT: PERHAPS THERE SHOULD BE A .usage METHOD TO RETURN A LIST OF THE ALLOWED INPUTS. cl.setsolve (type="T",t=10.0,phaseonly=T,table='J0739+016.phasecal',append=F); Bingo. Note that the daily documentation is correct. cl.solve(); T Jones Final fit per unit weight = 7.16277 Jy, after 12 iterations T Jones Did not converge for 30 time intervals Seems OK. cl.setapply (type="T", t=10.0, table='J0739+016.phasecal'); cl.correct(); This should have corrected the visibilities if I am not mistaken. cl.close() Lets look at the solutions gtable := table('J0739+016.phasecal') gains := gtable.getcolslice(columnname='GAIN',blc=[1,1,1],trc=[1,1,2]) goods := gtable.getcol(columnname='SOLUTION_OK') gtable.done() print shape(gains) [1 1 2 6500] print shape(goods) [10 6500] gmask := goods[1,] ggain := gains[1,1,1,gmask] print shape(ggain) 550 ngain := len(ggain) phases := arg( ggain )*180.0/pi pharms := stddev( phases ) phamin := min( phases ) phamax := max( phases ) Make glish record of phase cal stats gstats := [ rms=dq.quantity(pharms,'deg'), min=dq.quantity(phamin,'deg'), max=dq.quantity(phamax,'deg'), number=ngain ] print gstats [rms=[value=45.2178759, unit=deg], min=[value=-171.198024, unit=deg], max=[value=172.500943, unit=deg], number=550]; This seems reasonable, though the rms phase correction of 45deg seems a little high. 9. Lets go back and clean again imgr.image(type='corrected',image='J0739+016.corr.dirty'); im := image('J0739+016.corr.dirty'); im.view(raster=T); im.done(); Looks about the same as before, but with Peak = 1.357 Jy/bm Rms = 66 mJy/bm (slightly different than before) imgr.make('J0739+016.corr.clean'); imgr.clean(mask='J0739+016.mask',model='J0739+016.corr.clean', image='J0739+016.corr.clean.restored', residual='J0739+016.corr.clean.residual', niter=300,gain=0.03); Clean used 300 iterations to get to a max residual of 0.0624659 im := image('J0739+016.corr.clean.restored'); im.view(raster=T); im.done() New Peak = 1.348 Jy/bm (cf. old 1.355 Jy/bm) Rms = 26 mJy/bm in cleaned region (cf. 30 mJy/bm) Rms = 60 mJy/bm outside mask (cf. 60 mJy/bm) imgr.done() This seems to have been successful. Next step is to look at some of the "blank" field data which should be dominated by the CMBR! POST-MORTEM: The imager/calibrater (and image) routines were successful even given the strange mount type that Tim put into the CBI header. However, this will eventually need to be fixed as this data was taken with the deck at fixed angle zero (so the parallactic angle is just the normal parallactic angle) but in general we will use deck rotations (effectively rotating the parallactic angle) to help improve the uv-coverage and diagnose systematics (eg. ground pickup). Session complete.