Reduction of Keck-NIRSPEC observations of Pluto

Section 1 Introduction

This web page describes the reduction of spectroscopic observations of Pluto from 2.8 to 3.7 mm recorded on Keck II (www2.keck.hawaii.edu) with the NIRSPEC instrument (www2.keck.hawaii.edu/inst/nirspec/nirspec.html). The data were recorded on August 12, 2001 by Eliot Young and Will Grundy.

Section 2 Wavelength Calibration

The first step is the wavelength calibration to determine the mapping from pixel space to wavelength. We used an image of a Neon lamp for this calibration and a list of Neon lines from the NIST database. We visually identified lines on this NIST list in the 2-D image of the Neon lamp. In order to have an initial guess as to the location of the lines in our image, we used the NIRSPEC Format Simulator. This is how the simulator looked for our setup. The purple text is our estimate of line locations:

Figure 1. The NIRSPEC simulator

Once we knew the approximate spectral range, we could visually identify the lines on the Neon lamp image. Figure 3 shows the image of the Neon lamp (on the left) and the list of lines from the NIST database. Three times, multiple lamp lines were not resolved. For two of the occasions, the two line centers were averaged – a reasonable solution as the lines had similar strengths. For the other case, the middle line of the three blended lines used.

Figure 2. The list of Neon lines and our Neon image.

Now that we have identified the lamp lines, we can establish the transformation from pixel space to wavelength space. Ultimately we would like our rectified images to have the wavelength axis correspond to increasing column numbers. Therefore, before doing the wavelength calibration or later the rectification of the raw images, we rotate the image by 270°, using the command rotate(data, 3) in IDL. This is automatically done in the rectification procedure used later.

For the polynomial relating x and y to wavelength, we used a cubic:

l (x, y) = a + b x + c y + d x y + e x² + f y² + gx²y + hxy² + ix³ + jy³

We used the IDL procedure GetLambdas.pro to generate of set of (wavelength, x, y) triplets for least-squares fitting. This procedure allows you to enter the wavelength for a line, the click on the line in another window displaying the lamp image. It uses a lower level function called GetClix.Pro which returns the coordinates of a pixel that you click on. We identified ~5 positions along each line. The x-coordinate of this list (remember that we rotated the image before doing GetLambdas) were refined by a 4-parameter Gaussian fit to improve the center position (using refinelines.pro). All IDL procedures written for this analysis can be found in an adjacent directory.

The list of wavelengths and refined pixel coordinates generated is input into a least squares routine to fit the parameters of the bivariate quadratic (or cubic, or quadratic) equation described above. The procedure is called WaveCal.Pro. It uses the method of least-squares with a single value decomposition to fit the parameters (a, b, c, d, e, and f). The procedure returns the fitted parameters, residuals, chi-squared and the wavelength map – a 2-D image describing the mapping from image coordinates to wavelength space. The residuals from our fit are shown in Figure 3. For each of the lamp lines, you see five different points corresponding to different slit positions. There is no apparent trend in wavelength or slit position that would indicate the need for a higher polynomial fit.

Figure 3. The residual from the cubic fit of the Neon lines.

Section 3 Rectification

Next we would like to rectify the image so that it has orthogonal coordinate of wavelength and slit position instead of x and y from the raw image. The slit position is the location of the trace peak at row=512 in the raw data. We denote slit position by the variable ‘p’.

In order to carry out the transformation from column and row to wavelength and slit position, we must know the mapping from column and row to wavelength and slit position. We have one piece from the wavelength calibration. Next we need to determine the transformation from x and y to slit position. We construct triplets of (slit position, x and y) to fit for the parameters using waveCal.pro. We do this using the trace of the star (BS6060) using A-B pair subtracted and flat-fielded images. The trace center is the column with the maximum signal in each row. The list of columns of maximum signal is robustly fit to a quadratic to determine a profile of the trace which is not sensitive to outliers (the yfit parameters from goodpoly). This is done in FitTrace.pro. For each star image, we have 1024 points that describe the trace. We chose 8 of the BS6060 images and constructed an array of 8192 (8x1024) triplets describing the slit position as a function of row and column. The slit position is simply the column from the goodpoly fit at row=512 (see “makePmap5”). The quadratic polynomial was not sufficient to characterize the transformation from row and column to slit position. The residuals showed a pattern that would be removed with a higher order polynomial hence we increased the order of the polynomial to a cubic and the fit was improved.

Figure 4. The residuals of the fit of slit position.

Now that we have wavelength and slit position as a function of x and y, we can turn it around to get x as a function of wavelength and slit position and y as a function of wavelength and slit position. The details are described in ‘convertToLambdaP8’. We construct a grid of x and y coordinates from 0 to 1024 with 10,000 points in each axis. Then using the fitted parameters from the wavelength calibration compute the wavelength for each grid point. Similarly, using the fitted parameters from the slit position cubic equation, compute the slit position for each grid point. Next construct a triplet of points to determine the mapping: (x, wavelength, slit position). You will have 10,000 triplets. With these we can do a least-squares fit to find the parameters describing the bivariate polynomial. We found that x as a function of wavelength and slit was best described by a cubic polynomial. The same procedure was repeated to derive a quartic polynomial relationship for y as a function of wavelength and slit position. A detailed script of the steps is given in convertToLambdaP8.

Now that we have the map parameters to convert from wavelength and slit position to x and y coordinates, we are able to construct our rectified image (see Rectify.pro). First, we flatten the data. Next, we construct a grid of wavelengths and slit positions that covers the range of our data and is oversampled by a factor of 8 compared to the data. The grid of wavelength and slit position is transformed into a grid in x and y using the parameters found above. We construct the rectified image using the grid x and y coordinates and then re-bin the new array back to the original size.

To be conservative, the wavelength and slit position ranges extended a little beyond our data. This led to extrapolated data outside the region of the image. Hence, we only used columns 38 to 996 of the rectified image.

The rectification procedure was repeated on the wavelength map image generated back in section 1. This gave us an indication of the quality of our rectification. We find the wavelength scale to vary by ~8x10^-5 microns (or about 1/10^th of a pixel) across the range of slit positions observed.

Figure 5. The difference in wavelength scale for slit position of 225 and 60.

Section 4 Removing the Background

Because of the fluctuation in the background signal with wavelength, the background was fit and subtracted on a column by column basis. For each column we consider two background “boxes” (boxes with a width of 1 pixel in column) equidistant from the center of the trace. The box location was chosen such that no object signal was present, only background. This was done by visually inspecting the Gaussian profile of the brightest star observations. An inner limit of 15 pixels and an outer limit of 25 pixels was chosen. A linearly sloping background was fit to the median pixel values in the background boxes of the AB subtracted images. The fitted column-wise linearly sloping background was subtracted from the AB subtracted images. The resulting background subtracted images were used for the reduction from this point on.

We also need a value for the background for each individual image to calculate the noise per pixel for the optimal extraction. These backgrounds were determined using the same method, just on an individual image instead of a AB subtracted image. The results were stored in a directory for later use.

Section 5 Constructing the Normalized Gaussians

The next step is the construct normalized Gaussians for the optimal extraction of the Pluto data. The procedure gaussFitByCol.Pro fits a 3-parameter Gaussian model to each column of the data given the center of the trace and the half-width of the number of rows to extract around the center for the fit. The trace center was identified by summing the first 300 columns (where the Pluto signal is the strongest) and locating the maximum row in the sum. The fitted Gaussian parameters are saved in a FITS file.

The Gaussian model parameters were slowly varying over the wavelength range. Figure 6 shows the fitted Gaussian centers and widths for four different images. This was a series of ABBA. Note in the A and B nod positions the trace of the center positions has a different shape. For the optimal extraction of the Pluto spectra, the normalized Gaussian is derived from a BS6060 image taken with a similar slit position and just prior to the set of Pluto images. We could not use the exact slit position because the star images were taken with slit positions of 90 and 225 pixels and the Pluto images have a slit positions near 60 and 198 pixels.

Figure 6. The Gaussian center and width for 4 different images.

We construct a smoothly-varying, normalized-Gaussian model for each BS6060 image by doing a robust quadratic fit (using the function goodpoly) to the parameters in Figure 6. The robust fit excluded data that was more than 3-sigma from the model. The model was normalized by the area under the Gaussian. This was done in MakeNormalizedGauss.pro.

For the optimal extraction, we also need an estimate of the noise in each pixel. We consider three sources of noise: source noise, background noise and read noise. In particular, the noise of a single image is:

s = SQRT(gain * DN_S + gain * DN_B + RN^2)/gain

The gain of the NIRSPEC system is 5 e^–/DN and the read noise (RN) is 33 e^– (from page 4 of www.astro.ucla.edu/~irlab/nirspec/nirspec.htm). DN_S is the counts of the object without the background (as calculated above) and DN_B is the counts of the background (also calculated previously). Hence the noise of a AB-subtracted pair is:

s = SQRT(s _A^2 + s _B^2)

Section 6 Optimally Extracting the Spectra

For each wavelength, a fit was performed (using madscalew.pro) to minimize the average weighted deviation of a function: y=a*x. For us, the normalized guassian model was the “x” and the (AB and background subtracted) data is “y” with weights on a per pixel basis given by 1/s^2, where s was described in the previous section. The resulting scaling factors, a, (one for each wavelength) is the optimally extracted spectrum.

To extract the BS6060 spectra, the choice of normalized Gaussian model was easy (remember, we have a normalized Gaussian model for each BS6060 image). We simply use the model derived from the data we want to extract. The fitting of the scale factors was done only on pixels within 22 pixels of the trace center. The optimal extraction is performed by the IDL procedure OptimalExtraction.pro.

To extract the Pluto spectra, (as mentioned previously) we chose the normalized Gaussian model derived from the BS6060 image that was recorded just before the series of Pluto images and has the closest slit position. For the Pluto extraction, we need to align the center of the normalized Gaussian model and the Pluto image. This is done by cross-correlation. We started by reading in a 45 pixel wide swath of the (background and nod subtracted) Pluto image centered on the Pluto trace center (as determined by the row with the max counts in a sum of ~300 columns). The same size area was subtracted around the BS6060 trace center.

Next we resample the data (both Pluto and BS6060) expanding the scale by a factor of 10. We compute the product of the two images shifting the star position incrementally by 1 pixel in row (corresponding to 1/10^th pixel in the raw scale). We find the relative centers of the two images by the shift that has the largest product. We apply the shift reduced by a factor of 10 to the normalized Gaussian model prior to performing the optimal extraction (described above). The procedure that carries this out is MakeShiftedNG.pro.

Section 7 Remove extinction and the solar contribution

Our standard star is an excellent solar analog and we will use it to remove extinction and solar lines. Using ATRAN, we identified wavelengths where the atmospheric transmission was greater than 0.97 over the span of observed airmass. These “clear atmosphere” windows were used to normalize the observed BS6060 spectra.

For each wavelength, we determined the extinction coefficients from a linear fit. The star observations bracketed the Pluto observations in airmass. Using the extinction coefficients, we computed a spectrum containing the solar analog signal with the extinction corresponding to the airmass of every Pluto observations. The Pluto data were normalized by this spectrum to remove the effect of extinction and solar lines. Details of this procedure are given in the document “extinction”.

The final result is the mean of 37 individual Pluto spectra binned in wavelength by 10 pixels.

Section 8 Error Estimate

Last we propagate the error. To do this, we normalize all the spectra to the level of the mean Pluto spectrum using madscale (unweighted) to determine the scaling factor. For each wavelength, we calculate the standard deviation of the 37 scaled spectra. The square root of this set of standard deviations summed in quadrature divided by the sqrt(N), where N=37, is the error in the mean spectrum. Details of the error propogation are given in errorProp4.

The function, bin.pro, bins the wavelength scale, spectrum and error array to the desired resolution. The resulting spectrum for Pluto in the 3 um region (binned in wavelength by 10) is shown below. There is a small gap in the spectrum due to a strong telluric feature at this wavelength that was difficult to correct for sufficiently.