Coadding data, or how to deal with long integrations

Reducing multiple observations of the same source is about as easy as reducing a single map, it is just a bit more time consuming. Some aspects of the reduction process, like despiking, now becomes simpler and more reliable.

All the basic steps up to and including sky removal are done the same way as in the example above, but now we can despike differently, because we have much more redundant data. After we have extinction corrected, calibrated, done sky removal and pointing corrections we can now try despike on all the data sets.

To demonstrate despike we have collected 10 jiggle maps with two or three integrations of HL Tau, one of our secondary calibrators. All the maps are taken in good weather conditions, but some suffer from sharp variations in sky noise. We have set the really noisy bolometers to bad for all maps, but in some cases we have still left bolometers with 1.5 - 2 times the average noise level. For all maps the chop throw was 120" in Azimuth. Note that if you chop on the array this despiking technique will not work very well, unless you use a fixed chop angle, or regrid in az which works very well for a point source or a centrally condensed spherically symmetric source. To make the reduction easier, and to enable us to go back and redo the despiking, we make a small ASCII file for the maps, which includes the name, weight and pointing shift - one line per map. Here we start without applying pointing corrections because the pointing errors are negligible. Neither do we apply weighting, because we will need despiked maps in order to compute the weights. In this example we call the file hl.inp. The first line in this file is hl39_lon_sky 1 0 0, where hl39_lon_sky is the file name and 1 is the weight. The two zeros that follow the weight are the pointing shifts in the coordinate frame that the data are rebinned to, which in this case RJ. Since HL Tau is a point source with negligible extended emission one could equally well have rebinned the whole data set in az. We are now ready to try despike.

% despike noloop
OUT_COORDS - Coordinate sys of output map; PL,AZ,NA,RB,RJ,RD or GA 
/'RJ'/ > 
SURF: output coordinates are FK5 J2000.0
REF - Name of first data file to be rebinned /'hltau_lon'/ > hl.inp
SURF: Reading file hl39_lon_sky
SURF: run 39 was a MAP observation of HLTau with JIGGLE sampling
SURF: file contains data for 4 exposure(s) in 3 integrations(s) in 1
measurement(s)

   ...............   list continues until it has read in all maps, 
then
 
SURF Input data: (name, weight, dx, dy)
   -- 1: hl39_lon_sky (1, 0, 0)

     ..............
   -- 10: hl42_lon_sky (1, 0, 0)
 
NSIGMA - Sigma levels at which to despike /3/ >
SMODE - Smoothing mode for clipping envelope (none,hann) /'hann'/ > 
DEVICE - Name of graphics device /@xwindows/ > 
DMODE - Display mode (Spiral, Xlinear, Ylinear, Diag1, Diag2) 
/'spiral'/ > 
XRANGE - X range of plot (! to continue) /[1,5112]/ > 1,2000
XRANGE - X range of plot (! to continue) /[1,2000]/ > !
SURF: 8 spikes detected in file hl39_lon_sky
...........................
SURF: 196 spikes detected in file hl42_lon_sky
OUT - Name of despiked data file /'hl39_lon_dsp'/ > 
...........................
OUT - Name of despiked data file /'hl42_lon_dsp'/ >

Figure: These are the first 2000 data points analyzed by despike. The observed data are plotted in green, median in black and the Hanning smoothed 3-$\sigma$ envelope is plotted in red. Any data points outside this limit will be flagged as bad.

The display for the first 2000 data points are shown in Fig. . Note that after you run despike the display may get messed up. It is normally sufficient to clear the display with gdclear and re-issue the the color table, e.g. lutbgyrw. If there are still problems, delete the AGI-file from your home directory (agi_xxx.sdf, where xxx is the name of your work station).

We now run rebin on each data set to determine the offsets we need to center the map, which will ensure that we get the sharpest possible image. We also determine the rms noise level for each map. Here it is a clear advantage to use GAIA. After displaying the map with your preferred color scheme and contrast level using the View menu, go to the Image-Analysis menu and click on image regions. Choose the polygon option of ARD regions and outline the region you want to use with the cursor. After it is done, click on Stats selected and you are done. For determining pointing offsets it is easier to use the KAPPA command centroid i.e. centroid search=19 cerror, where we adjust the search parameter to minimize the fit errors. For 850 $\mu$m and a 1 arcsec cell setting search to 17 or 19 generally gives a good result. After we have determined the noise level and position offsets for each map we compute the weights for each map. An easy way to do this is to set the weight for the first map equal to one. All other maps are then weighted relative to this map. If we call the integration times t$_1$ and t$_2$, and the corresponding noise values n$_1$ and n$_2$ for map 1 and map 2 respectively, then we can compute the weighting factor from the following equation:

$$weight = (n_1/n_2)^2\times(t_1/t_2)$$ (3)

We now edit in the weights and position offset in hl.inp and run despike again on the original data set. The end result appears good enough.

Now we can form the final co-add of the data very easily, since we have already computed the weights and position shifts that we need to get the best S/N and sharpest image from our data. The only thing we need to do is to change the file extensions from sky to dsp, the name we used for despiked data files. We rename the input file to hl.dsp in case we would like to redo the despiking. We now run this file through rebin, i.e.,

% rebin noloop
REBIN_METHOD - Rebinning method to be used /'LINEAR'/ > 
SURF: Initializing LINEAR weighting functions
OUT_COORDS - Coordinate sys of output map; PL,AZ,NA,RB,RJ,RD or GA 
/'RJ'/ > 
SURF: output coordinates are FK5 J2000.0
REF - Name of first data file to be rebinned /'hl42_lon_dsp'/ > 
hl.inp
SURF: Reading file hl39_lon_dsp
SURF: run 39 was a MAP observation of HLTau with JIGGLE sampling
SURF: file contains data for 4 exposure(s) in 3 integrations(s) in 1
measurement(s)
..................
SURF: Reading file hl42_lon_dsp
SURF: run 42 was a MAP observation of HLTau with JIGGLE sampling
SURF: file contains data for 4 exposure(s) in 3 integrations(s) in 1
measurement(s)
 
SURF Input data: (name, weight, dx, dy)
   -- 1: hl39_lon_dsp (1, 0.08, 1.19)
   -- 2: hl64_lon_dsp (0.4, 0.64, 0.75)
   -- 3: hl41_lon_dsp (0.9, -0.23, 2.7)
   -- 4: hl55_lon_dsp (1, -0.86, 0.38)
   -- 5: hl61_lon_dsp (0.8, 0.14, 1.04)
   -- 6: hl75_lon_dsp (0.75, 0.4, -0.13)
   -- 7: hl84_lon_dsp (0.5, 0.42, 0.83)
   -- 8: hl60_lon_dsp (0.9, 0.11, -1.41)
   -- 9: hl23_lon_dsp (0.75, 0.1, 0.46)
   -- 10: hl42_lon_dsp (0.3, 0.75, 1.06)

LONG_OUT - Longitude of of output map centre in hh (or dd) mm ss.ss 
format
/'+04 31 38.46'/ > 
LAT_OUT - Latitude of output map centre in dd mm ss.ss format 
/'+ 18 13 57.8'/ > 
OUT_OBJECT - Object name for output map /'HLTau'/ > 
PIXSIZE_OUT - Size of pixels in output map (arcsec) /3/ > 1
SIZE - Number of pixels in output map (NX,NY) /[194,192]/ > 
OUT - Name of file to contain rebinned map /'HLTau'/ > hltau_lon
WTFN_REGRID: Beginning regrid process
WTFN_REGRID: Entering second rebin phase (T = 3.113513 seconds)
WTFN_REGRID: Entering third rebin phase (T = 14.55025 seconds)
WTFN_REGRID: Regrid complete. Elapsed time = 15.23921 seconds

We now have our final co-added map, which is almost as good as it can get from the data we had available.

If you have extremely large data sets, then you may find that your computer runs out of memory in rebin. The recommended solution, if you cannot find a more powerful workstation, is to split the data set in two parts, run each separately through rebin, and co-add the two final maps using the KAPPA command add or the task makemos in CCDPACK.