Common data reduction

Now the data processing can begin. We will start by running reduce_switch in order to subtract the off from the on and to split the data array of the raw data into separate components.

% reduce_switch 59
SURF: Opening apr8_dem_0059 in /scuba/observe/apr8/dem
SURF: run 59 was a MAP observation of object 3c279
SURF: file contains data for 2 switch(es) in 4 exposure(s) in 3 integration(s)
in 1 measurement(s)
OUT - Name of output file to contain reduced switch data /'o59'/ >

or with the full file specification:

% reduce_switch apr8_dem_0059
SURF: Opening apr8_dem_0059 in /scuba/observe/apr8/dem
SURF: run 59 was a MAP observation of object 3c279
SURF: file contains data for 2 switch(es) in 4 exposure(s) in 3 integration(s)
in 1 measurement(s)
OUT - Name of output file to contain reduced switch data /'o59'/ >

In this example the calibrator signal has not been used and any datum from which more than 5 spikes were removed by the transputers is marked bad (these are the default settings). The processed data are then written to file o59.sdf.

In this case we need to change the flatfield file (since the flatfield was updated after the data were taken) using change_flat:

% change_flat
IN - Name of input file containing demodulated map data /@o59/ > 
SURF: run 59 was a MAP observation of 3c279
NEW_FLAT - The name of the file containing the new flat-field > photflat1.dat

The next task is to flatfield the data:

% flatfield o59 o59_flat
SURF: run 59 was a MAP observation of 3c279
SURF: applying flatfield from photflat1.dat

If the input and output files are not specified on the command line they will be requested.

The data can now be corrected for airmass (elevation) and sky opacity by using extinction. According to the skydip observation taken prior to the map, the tau at 850 $\mu$m is 0.220 and a skydip taken after the map shows it was 0.187. If extinction is given two $\tau$ values from different times then the actual $\tau$ for each jiggle will be calculated by linear interpolation - obviously this assumes that the $\tau$ varied linearly with time. For this example we will assume the $\tau$ variations are correct in order to demonstrate the principle:

% extinction
IN - Name of NDF containing demodulated data /@o59_flat/ > 
SURF: run 59 was a MAP observation with JIGGLE sampling of object 3c279
SURF: file contains data for 4 exposure(s) in 3 integration(s) in 1
measurement(s)
SURF: observation started at sidereal time 12 19 59 and ended at 12 28 40
SURF: file contains data for the following sub-instrument(s)
 - SHORT with filter 450
 - LONG with filter 850
SUB_INSTRUMENT - Name of sub-instrument to be extinction corrected /'SHORT'/ > l
FIRST_TAU - First zenith sky opacity measured /0/ > 0.22
FIRST_LST - Sidereal time of first opacity measurement; hh mm ss.ss /'0.0'/ > 11 54
SECOND_TAU - Second zenith sky opacity measured /0.22/ > 0.187
SECOND_LST - Sidereal time of second opacity measurement; hh mm ss.ss /'11 54'/ > 16 10
OUT - Name of output NDF /'o59_lon_ext'/ >

The arrays are separated at this point (since the extinction correction would be different). In this case the LONG-wave array was selected; extinction would have to be re-run to select the SHORT-wave array (the question is not asked if only one sub-instrument is present).

Some comment is probably required for the use of LST for the tau measurements - this is not as bad as it sounds. The skydip task prints the LST of the skydip and extinction prints the LST of the observation; in many cases it is known that the tau value was taken a certain time before or after the observation so this value can simply be added. In addition, most of the time a constant tau value is used and for the case of a constant tau the LST is irrelevant.

If desired, sky noise can now be removed. Sky signal can be identified in two ways: firstly using bolometers that are known to be looking at sky (implemented in remsky) and secondly, using a model of the source structure to enable the sky signal to be calculated with the source subtracted from the data (implemented in calcsky). In this section we will examine the first method since this is the simplest and can be used for jiggle observations. The more complex approach will be dealt with in § where sky removal from scan map data is discussed. remsky works in a very simplistic way: Sky bolometers are specified, each jiggle is then analysed in turn, the average value for the sky bolometers (either MEDIAN or MEAN) is then removed from the entire array. At present it is not possible to specify sky regions, only sky bolometers can be specified. This may cause problems with extended sources (rebin the map in NA coordinates initially to find the sky bolometers). remsky should normally be run after rebin in order to choose sky bolometers that are really looking at sky, for this example we will skip that step. remsky is sufficient for the mapping of compact sources in jiggle mode and for photometry; if the source structure is complex then calcsky should be considered (§).

% remsky
IN - Name of input file containing demodulated map data /@w48_newrebin/ > o59_lon_ext
SURF: run 59 was a MAP observation with JIGGLE sampling of object 3c279
OUT - Name of output file /'o59_lon_sky'/ > 
BOLOMETERS - The Sky bolometers, [a1,a2] for an array /['all','-h7']/ > 
SURF: Using 36 sky bolometers
MODE - Sky removal mode /'median'/ > 
Adding mean background level onto data (value=1.5316721E-6)

In this example we have used all the bolometers except for the central pixel (H7) and then used median sky removal for each jiggle. The average background level has also been added back onto the data.

Figure: The 3C279 data after processing through extinction and remsky. The next stage is to regrid the data using rebin. The source can clearly be seen in bolometer 19 (H7). The negative stripes are indicating that the chop throw was smaller than the array.

The output of extinction or remsky can be displayed using, say, KAPPA display to see whether some integrations or bolometers should be marked bad.⁸ Sometimes bad bolometers can only be identified after a rebin/scuover phase. The output so far can be seen in figure - the axes are labeled with bolometer number along the X-axis and integration number up the Y-axis.

Now that the data have been extinction corrected and, optionally, processed with remsky, the data reduction path diverges according to the type of observation. Map making⁹ (jiggle and scan) and photometry will be dealt with separately. Note that scuquick can be used to automate some of the more repetitive tasks during this stage of the data reduction process.

Before diverging though, we should first take a diversion into the question of despiking.