@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@ First Pass Reduction: IC10 mc7 27 June 2001 @%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@ by steve The following notes refer to my BIMA summer school reduction of an observation of molecular cloud complex 7 in IC10 (Wilson 1995). The goal here was to get as high a S/N as possible to augment Adam Leroy's map. This is my first ever use of MIRIAD, so things aren't so slick. Also, the data was taken (mostly) in the rain. For the record, it was a 3mm observation, with CO centered in the upper side band. The resulting map should be 120" square with a spatial resolution of 15". The spectral bandwidth was XXX MHz with a velocity resoltution of 1 km/s. The combined grade for the observation was a C. The reduction notes loosely follow the suggestion of Mel Wright in lecture on 27 June 2001. ------------------------------------------ Stuff You Should Know Before You Start ------------------------------------------ 1. To get help on a miriad task: > mirhelp taskname 2. I'm running miriad at the command line. In Berkeley, I'll have to set the appropriate paths so that the commands are known. 3. Miriad deals with data structures. These are somtime directories, sometimes not. They're sometimes ASCII, but more often not. To remove something (like "gains", for instance), use > \rm -r bad_thing It might not get killed otherwise. 4. If a miriad keyword value is to contain a parenthesis, but the keyword=value string in single quotes. 5. "options" is plural. This is bound to screw you up. -------------------------- Baby's First Reduction -------------------------- Beware: I'm guessing about a lot of this! 1. Plot amplitude vs. time for calibrator to look for bad data (e.g. too high/low; umps). Interested in upper side band (USB) average, as that's where the CO line is. My phase calibrator was 0228+673. > uvplt vis=0228+673 line=wide,1,2 device=/xw axis=time,amp They all look pretty smooth, hovering around unity. There is one catastrophically hight point -- notice the y-axis scale. If this were the data, we'd need to take that out. The inverse fourier transform would smear that all over the map. Since this is just the calibrator, we can ignore it. 2. Plot phase vs. time for calibrator. > uvplt vis=0228+673 line=wide,1,2 device=/xw axis=time,phase These should be somewhat tight. Mine, on the other hand are really noisy. This is why you need to fit the phase gains and correct the data! 3. Fit the phase gains of the calibrator. The interval should be sufficient to interpolate between calibrator observations -- 20 or 30 minutes is recommended for a script that has 20 minutes on the source and 5 minutes on the calibrator. > selfcal vis=0228+673 line=wide,1,2 interval=20 refant=1 Here, we set the reference antenna (that with arbitrarily zero phase) to be antenna 1. This step made a data structure called "gains" under the phase calibrator directory. 4. Check your phase gain calibration by plotting phase vs. time for the calibrator again. UVPLT automatically applies these gains. > uvplt vis=0228+673 line=wide,1,2 device=/xw axis=time,phase It only looks a little better. Ugh. However, there are so many points that it's hard to tell what's going on. One way to check is by binning up the points. Check out: > uvplt vis=0228+673 line=wide,1,2 device=/xw axis=time,phase average=5 These should hover around zero, and they do, sort of. Keep in mind that you have a signal of about 1 Jy being added vectorially to an error of something like 5 Jy with random phase. Still, this ain't bad. And there isn't much else you can do, anyway. 5. Another check on the phase gain calibration is to make a map out of the calibrator. You'll make a map out of the beam as well. > invert vis=0228+673 map=0228+673.mp beam=0228+673.bm line=wide,1,2 Now look at the map and beam with: > cgdisp in=0228+673.mp device=/xw > cgdisp in=0228+673.bm device=/xw This looks okay -- it's a point source, after all. And Dick says that the beam looks okay. For comparison, do the inversion without applying the phase calibration. The beam should look the same -- that has nothing to do with the phase calibration. The data, on the other hand, should look like random noise. It's a bunch of out-of-phase fourier components! > invert vis=0228+673 map=0228+673.mp_nocal beam=0228+673.bm_nocal line=wide,1,2 options=nocal > cgdisp in=0228+673.mp_nocal device=/xw > cgdisp in=0228+673.bm_nocal device=/xw Notice that I appended "nocal". Otherwise, invert would have crashed out, not wanting to overwrite the previous map and beam. 6. Yet another thing to check is a plot of the gains with gpplot. I think the only difference between this and what we did in (4) is that the phases don't wrap, and you're looking at the relative phases between the antennes, rather than look at 36 out 45 baselines. > gpplt vis=0228+673 device=/xw yaxis=phase Notice that antenna 1 is flat. That's 'cause it's the reference! 7. In any case, it's time to apply these phase gains to the data. Do this by copying the gains subdirectory under the phase calibrator directory (created by the selfcal up in step 3) over to the data directory. Doing this the miriad way invokes gpcopy: > gpcopy vis=0228+673 out=IC10 Notice that where there wasn't a gains subdirectory under the source directory (IC10/), there now is. 8. We now want to make a dirty map of the source. But we want to see the spectral line of interest, so we should pre-determine the interesting velocity... Do this with uvlist with the "spectra" option, the idea being to gain information about the spectral windows. > uvlist vis=IC10 options=spectra This gives: rest frequency : 115.27120 115.27120 starting channel : 1 257 number of channels : 256 256 starting frequency : 112.11579 115.35551 frequency interval : -0.00039 0.00039 starting velocity : 8229.617 -195.441 ending velocity : 8488.657 -454.480 velocity interval : 1.016 -1.016 invert vis=IC10 options=systemp 'select=window(1)' map=junk.mp beam=beam.mp invert vis=IC10 options=systemp 'select=window(1)' map=junk.mp beam=junk.bm line=wide,1,1 cgdisp in=junk.mp device=/xw invert vis=IC10 options=systemp line=velocity,10,-400,20,20 map=junk.mp beam=beam.mp \rm -r junk.mp \rm -r junk.bm invert vis=IC10 options=systemp line=velocity,10,-400,20,20 map=junk.mp beam=junk.bm cgdisp in=junk.mp device=/xw options=3value csize=1,1,1,1 This makes maps 20 km/s wide, spanning -400 km/s to -200 km/s. Looks like there's some flux in the 4th map (-340 to -320? which is right for v=-330???). Looks wavy though, which could mean that some hight point is corrupting things... Let's try making just that swath centered on 330. \rm -r junk.mp \rm -r junk.bm invert vis=IC10 options=systemp line=velocity,1,-340,20,20 map=junk.mp beam=junk.bm cgdisp in=junk.mp device=/xw options=3value csize=1,1,1,1 And now let's trying cutting out the high T_sys data at the beginning of the track. \rm -r junk.mp \rm -r junk.bm invert vis=IC10 options=systemp line=velocity,1,-340,20,20 map=junk.mp beam=junk.bm 'select=time(09:26:01.0,16:24:54.0)' cgdisp in=junk.mp device=/xw options=3value csize=1,1,1,1 *** USEFUL STUFF WITH RICK *** uvflag vis=IC10 flagval=f 'select=time(09:26:01.0,16:24:54.0),ant(4)' ###flag that starred data corset freq=115.271204 iffreq=400 cormode=1 corf=447.632 corbw=100 device=/ps ###make the corset plot uvjima IC10 ###shortcut for amp vs. time (4 to a page) for checking the flagging invert vis=IC10 options=systemp line=velocity,1,-340,20,20 map=junk.mp beam=junk.bm \ cell=2 imsize=128 ### choose better cell size / total angular extent, use better display xmtv & \rm -r junk.mp ; \rm -r junk.bm invert vis=IC10 options=systemp line=velocity,1,-340,20,20 map=junk.mp beam=junk.bm cell=2 imsize=128 tvdisp in=junk.bm server=xmtv@bima2 tvdisp in=junk.mp server=xmtv@bima2 ### do same, but make data cube \rm -r junk.mp ; \rm -r junk.bm invert vis=IC10 options=systemp line=velocity,10,-400,20,20 map=junk.mp beam=junk.bm cell=2 imsize=128 tvdisp in=junk.mp server=xmtv@bima2 'region=images(1,1)' cgdisp in=junk.mp device=/xw options=3value csize=1,1,1,1 ### let's check out a spectrum. first make a cube with better resolution, ### then use velplot. Actually, first test velplot on the coarse cube above. cgdisp in=junk.mp device=/xw options=3value csize=1,1,1,1 'region=images(4,4)' velplot in=junk.mp device=/xw ### This is not straightforward. Start it up, as above. You want to pick your spectrum position interactively off of a velocity integrated map, so hit V. You don't want help, so hit return. You want a straight average, so hit return. You don't want blanking, so hit return. You want to just see your lines, so enter L for list. Now tell it 1,-340,-320 for "1 map (#1) over that velocity" . You don't want any more maps, so hit return. That's a fine display, so hit return. "L" gives a list of your options in the interactive map. "S" picks off the points for spectra. "E" exits when you're done. We don't want another plot device, so hit return. At the main menu, we now want to see our spectra. Hit "S". No help, so hit return. We already have positions, so don't enter a grid -- hit return. We don't want to edit the positions, so hit return. No convolution, so hit return. No smoothing, so hit return. No Gaussians, no scaling -- return! return! Fine windowing -return! There's your spectrum. You can make an ASCII file if you like. Now make your finer data cube: \rm -r junk.mp ; \rm -r junk.bm invert vis=IC10 options=systemp line=velocity,40,-350,1,1 map=junk.mp beam=junk.bm cell=2 imsize=128 cgdisp in=junk.mp device=/xw options=3value csize=1,1,1,1 velplot in=junk.mp device=/xw ### Alright -- what you'd like to do here is flag bad data. But I need some guidance with that, so what I'm going to do instead is press blindly on. The next thing to do is clean, I guess. clean map=junk.mp beam=junk.bm out=junk_c.mp cgdisp in=junk.mp device=/xw options=3value csize=1,1,1,1 cgdisp in=junk.mp device=/xw options=3value csize=1,1,1,1 'region=images(18,18)' ### Wow -- it's grey with a couple of pixels. Let's see what restore does. I'm going to use the standard mode, and I'm going to let the fwhm be chosen from the beam header. You want the gaussian to be about the size of the primary lobe of the beam; this will fit a gaussian to the beam, which ought to do the trick. Check it: restor model=junk_c.mp beam=junk.bm map=junk.mp mode=clean out=junk_cr.mp cgdisp in=junk.mp device=/xw options=3value csize=1,1,1,1 cgdisp in=junk.mp device=/xw options=3value csize=1,1,1,1 'region=images(18,18)' ### Well, that kinda looks like what we started with. Might as well clean and restore the map of the calibrator, while I'm here. Of course, this is gonna have those fatty pixels... clean map=0228+673.mp beam=0228+673.bm out=0228+673_c.mp cgdisp in=0228+673_c.mp device=/xw restor model=0228+673_c.mp beam=0228+673.bm map=0228+673.mp mode=clean out=0228+673_cr.mp cgdisp in=0228+673_cr.mp device=/xw ### Let's remake this with the better pixels: \rm -r 0228+673.bm \rm -r 0228+673.mp \rm -r 0228+673_c.mp \rm -r 0228+673_cr.mp \rm -r 0228+673_cr.mpcgdisp invert vis=0228+673 map=0228+673.mp beam=0228+673.bm line=wide,1,2 cell=2 imsize=128 cgdisp in=0228+673.mp device=/xw clean map=0228+673.mp beam=0228+673.bm out=0228+673_c.mp cgdisp in=0228+673_c.mp device=/xw restor model=0228+673_c.mp beam=0228+673.bm map=0228+673.mp mode=clean out=0228+673_cr.mp cgdisp in=0228+673_cr.mp device=/xw velplot in=0228+673_cr.mp device=/xw