telluric -- remove telluric features from 1D spectra
Telluric calibration spectra are shifted and scaled to best divide out telluric features from data spectra. This may be done non-interactively to minimize the RMS in some region or regions of the data spectra and interactively with a graphically search.
telluric input output cal
Input one dimensional spectra are corrected to remove telluric features by dividing by shifted and scaled calibration spectra. The calibration spectra are generally of hot, nearly featureless stars; hence this procedure is sometimes referred to as a B-star correction. The shifting allows for possible small shifts or errors in the dispersion zeropoints. The intensity scaling allows for differences in the airmass and variations in the abundance of the telluric species. The intensity scaling uses Beer's law which is the approximation that the change in absorption with abundance is an exponential relation.
The following describes the correction. Let J(x_i) be the calibration spectrum at a set of pixels x_i. An interpolation function is fit to this spectrum to give J(x). The shifted and scaled calibration function is then
(1) J'(x) = max (threshold, J(x+dx)) ** (A / A_cal * scale)
where dx is the pixel shift parameter, A is the airmass of the input spectrum, A_cal is the airmass of the calibration spectrum, and scale is the scale parameter. The operator "**" is exponentiation. The max operation limits the calibration spectrum to be greater than or equal to the specified threshold value. If the calibration value is ever less than or equal to zero then the task will quit with a warning error.
The output corrected spectrum is then computed as
(2) I'(x_i) = I(x_i) / (J'(x_i) / <J'>)
where I' is the corrected spectrum, I is the input spectrum, and <J'> is the mean of the shifted and scaled calibration spectrum to keep the output intensities comparable to the input spectrum. The value of <J'> is printed in the output as the "normalization". If the spectra are dispersion calibrated, possibly with different dispersion parameters, then the x values in (2) from the input spectrum are converted to matching pixels in the calibration spectrum using the dispersion functions of the two spectra.
The purpose of this task is to determine the best values of the shift and scale parameters dx and scale. There are automatic and interactive methods provided. The automatic methods are cross-correlation of the calibration and input spectra to find a shift and an iterative search for the in both shift and scale that minimizes the RMS of I' in some region. The automatic methods are performed first, if selected, followed by the interactive, graphical step. The following describes the steps in the order in which they occur.
The initial values of the shift and scale are set by the parameters shift and scale for the first spectrum. After that the values determined for the previous spectrum, those actually applied to correcting that spectrum, are used as the initial values for the next spectrum. The search steps and sample regions are also initialized by task parameters but may be modified during the interactive step and the modified values apply to subsequent spectra.
If the xcorr parameter is yes and the lag parameter is not zero the calibration spectrum is cross-correlated against the input spectrum. Each spectrum is prepared as follows. A large scale continuum is fit by a quadratic chebyshev using 5 iterations of sigma clipping with a clipping factor of 3 sigma below the fit and 1 sigma above the fit and rejecting the deviant points along with one pixel on either side. This attempts to eliminate the effects of absorption lines. The continuum fit is subtracted from the spectrum and the spectrum is extended and tapered by a cosine function of length given by the lag parameter.
The prepared spectra are then cross-correlated by shifting the calibration spectrum plus and minus the specified lag amount about the current shift value. Only the regions in the input spectrum specified by the sample regions parameter are used in the correlation. This produces a correlation profile whose peak defines the relative shift between the two spectra. The current shift value is updated. This method assumes the common telluric features dominate within the specified sample regions. The lag size should be roughly the profile widths of the telluric features.
If the tweakrms parameter is yes and dshift is greater than zero trial corrections at the current shift value and plus and minus one shift step with the scale value fixed at its current value are made and the RMS in the sample regions computed. If the RMS is smallest at the current shift value the shift step is divided in half otherwise the current shift value is set to the shift with the lowest RMS. The process is then repeated with the new shift and shift step values. This continues until either the shift step is less than 0.01 pixels or the shift is more than two pixels from the initial shift. In the latter case the final shift is reset to the original shift.
The scale factor is then varied if dscale is greater than zero by the scale step at a fixed shift in the same way as above to search for a smaller RMS in the sample regions. This search terminates when the scale step is less than 0.01 or if the scale value has departed by 100% of the initial value. In the latter case the scale value is left unchanged.
The search over the shifts and scales is repeated a second time after which the tweak algorithm terminates.
After the optional cross-correlation and tweak steps the interactive search mode may be entered. This occurs if interactive = yes. A query is asking whether to search interactively. The answers may be "no", "yes", "NO", or "YES". The lower case answers apply to the current spectrum and the upper case answers apply to all subsequent spectra. This means that if an answer of "NO" or "YES" is given then there will be no further queries for the remaining input spectra.
If the interactive step is selected a graph of three candidate corrections for the input spectrum is displayed. There also may be a graph of the calibration or input spectrum shown for reference. Initially the calibration spectrum is displayed. The additional graph may be toggled off and on and between the input and calibration spectra with the c and d keys. The three candidate corrected spectra will be with the current shift and scale in the middle and plus or minus one step in either the shift or scale. Initially the spectra will be at different scale values. Information about the current shift and scale and the step used is given in the graph title.
One may toggle between shift steps and scale steps with the x (for shift) or y (for scale) keys. The RMS in the title is the RMS within the currently defined sample regions. If one of the step values is zero then a display of different values of that parameter will not be selected. The step size will need to be set with a colon command to search in that parameter.
If x is typed when the three spectra are at different shifts then the nearest spectrum to the y cursor at the x cursor position will be selected. If the central spectrum is selected the step size is divided in half otherwise the current shift is changed and the selected spectrum becomes the middle spectrum. Three new spectra are then shown. The same applies if y is typed when the three spectra are at different scales. This allows an interactive search similar to the iterative tweakrms method described previously except the user can use whatever criteria is desired to search for the best scale and shift.
There are additional keystrokes and colon commands to set or change sample regions, reset the current shift, scale, and step sizes, expand the step size in the current mode, adjust the offsets between the spectra, and get help. The w key and GTOOLS colon commands are available to window the graphs. Any changes in the x limits apply to both graphs while y limit adjustments apply to the graph pointed to by the cursor.
Two other commands require a short explanation. The a key may be used to run the tweakrms algorithm starting from the current shift, scale, and steps and the current sample regions. This allows one to graphically set or reset the sample regions before doing the RMS minimization. The ":smooth" command and associated smooth task parameter allow the corrected spectra to be displayed with a boxcar smoothing to better see faint features in noise. It is important to realize that the smoothing is only done on the displayed spectra. The telluric correction and computed RMS are done in the unsmoothed data.
After the interactive step is quit with q or if the interactive step is not done then the final output spectrum is computed and written to the output image. A brief log output is printed for each spectrum.
? - print help a - automatic RMS minimization within sample regions c - toggle calibration spectrum display d - toggle data spectrum display e - expand (double) the step for the current selection q - quit r - redraw the graphs s - add or reset sample regions w - window commands (see :/help for additional information) x - graph and select from corrected shifted candidates y - graph and select from corrected scaled candidates :help - print help :shift [value] - print or reset the current shift :scale [value] - print or reset the current scale :dshift [value] - print or reset the current shift step :dscale [value] - print or reset the current scale step :offset [value] - print or reset the current offset between spectra :sample [value] - print or reset the sample regions :smooth [value] - print or reset the smoothing box size
1. To interactively search for a best correction with the default cross-correlation and tweak steps:
cl> telluric spec001.ms telspec001.ms spec005.ms
2. To search only for a scale factor:
cl> telluric spec001.ms telspec001.ms spec005.ms xcorr- dshift=0.
3. To processes a set of spectra non-interactively with the same calibration spectrum and to replace the input spectra with the corrected spectra and log the processing:
cl> telluric spec* "" calspec inter- > log
4. To apply the simplest scaling by the ratio of the airmasses alone:
cl> telluric spec* tel//spec* calspec inter- xcorr- tweak- inter- \ >>> scale=1. shift=0.