implement GNIRSSpect._apply_wavelength_model_bounds() and remove _get…

…_cenwave_accuracy from code (except ATRAN stuff: TODO!)

geminidr/core/primitives_spect.py

@@ -5285,10 +5285,14 @@ def _get_actual_cenwave(self, ext, asMicrometers=False,
 
         return actual_cenwave
 
-
     def _get_cenwave_accuracy(self, ext):
-        # Accuracy of central wavelength (nm) for a given instrument/setup.
-        return 10
+        # TODO: remove this
+        dispaxis = 2 - ext.dispersion_axis()
+        npix = ext.shape[dispaxis]
+        w1, w2 = am.get_named_submodel(ext.wcs.forward_transform, "WAVE").copy()([0, npix-1])
+        m_wave = models.Chebyshev1D(degree=1, c0=0.5*(w1+w2), c1=0.5*(w2-w1))
+        self._apply_wavelength_model_bounds(m_wave, ext)
+        return 0.5 * abs(np.diff(m_wave.c0.bounds)[0])
 
     def _apply_wavelength_model_bounds(self, model=None, ext=None):
         # Apply bounds to an astropy.modeling.models.Chebyshev1D to indicate

geminidr/gmos/primitives_gmos_spect.py

@@ -591,10 +591,6 @@ def _get_arc_linelist(self, waves=None, ext=None):
                                 'CuAr_GMOS{}.dat'.format('_mixord' if use_second_order else ''))
         return wavecal.LineList(filename)
 
-    def _get_cenwave_accuracy(self, ext=None):
-        # Assumed accuracy of central wavelength in nm for a given instrument/setup.
-        return 10
-
     def _apply_wavelength_model_bounds(self, model=None, ext=None):
         # Apply bounds to an astropy.modeling.models.Chebyshev1D to indicate
         # the range of parameter space to explore

geminidr/gnirs/primitives_gnirs_spect.py

@@ -427,20 +427,19 @@ def _get_arc_linelist(self, ext, waves=None):
 
         return wavecal.LineList(filename)
 
-
-    def _get_cenwave_accuracy(self, ext):
-        # Accuracy of central wavelength (nm) for a given instrument/setup.
-        # According to GNIRS instrument pages "wavelength settings are accurate
-        # to better than 5 percent of the wavelength coverage".
-        # However using 7% covers more cases. For the arcs dc0=10 works just fine for all modes.
-
-        mband = ext.filter_name(pretty=True).startswith('M')
-        lband = ext.filter_name(pretty=True).startswith('L')
-        dispaxis = 2 - ext.dispersion_axis()  # python sense
+    def _apply_wavelength_model_bounds(self, model=None, ext=None):
+        # Apply bounds to an astropy.modeling.models.Chebyshev1D to indicate
+        # the range of parameter space to explore
+        dispaxis = 2 - ext.dispersion_axis()
         npix = ext.shape[dispaxis]
-
-        if 'ARC' in ext.tags or not (mband or lband):
-            dcenwave = 10
-        else:
-            dcenwave = abs(ext.dispersion(asNanometers=True)) * npix * 0.07
-        return dcenwave
+        for i, (pname, pvalue) in enumerate(zip(model.param_names, model.parameters)):
+            if i == 0:  # central wavelength
+                if 'ARC' in ext.tags or ext.filter_name(pretty=True)[0] in 'LM':
+                    prange = 10
+                else:
+                    prange = abs(ext.dispersion(asNanometers=True)) * npix * 0.07
+            elif i == 1:  # half the wavelength extent (~dispersion)
+                prange = 0.02 * abs(pvalue)
+            else:  # higher-order terms
+                prange = 1
+            getattr(model, pname).bounds = (pvalue - prange, pvalue + prange)

gempy/library/wavecal.py

@@ -439,8 +439,7 @@ class holding parameters for the UI, passed from the primitive's Config
         init_models, ui_params.toDict(), peaks=peaks,
         peak_weights=weights[ui_params.weighting],
         linelist=input_data["linelist"], fwidth=fwidth, kdsigma=kdsigma, k=k,
-        dcenwave = input_data["cenwave_accuracy"],
-        bounds_setter = input_data["bounds_setter"], filename = ext.filename)
+        bounds_setter=input_data["bounds_setter"], filename=ext.filename)
 
     input_data["fit"] = fit1d
     return input_data, fit1d, acceptable_fit
@@ -485,7 +484,6 @@ def get_all_input_data(ext, p, config, linelist=None, bad_bits=0,
     "linelist" : LineList object
     "fwidth" : feature width (pixels)
     "location" : extraction location (if 2D spectrum)
-    "cenwave_accuracy" : accuracy of the central wavelength
      "bounds_setter" : a callable to set the uncertainty on polynomial parameters
     """
     cenwave = config["central_wavelength"]
@@ -597,18 +595,15 @@ def get_all_input_data(ext, p, config, linelist=None, bad_bits=0,
                     log.warning(f"{i}. Offset {m.right.offset_0.value} "
                                 f"scale {m.right.factor_1.value}")
 
-    # Get the accuracy of the central wavelength
-    dcenwave = p._get_cenwave_accuracy(ext=ext)
-
     return {"spectrum": np.ma.masked_array(data, mask=mask),
             "init_models": m_init, "peaks": peaks, "weights": weights,
             "linelist": linelist, "fwidth": fwidth, "location": location,
-            "cenwave_accuracy" : dcenwave, "refplot_data": refplot_dict,
+            "refplot_data": refplot_dict,
             "bounds_setter": partial(p._apply_wavelength_model_bounds, ext=ext)}
 
 def find_solution(init_models, config, peaks=None, peak_weights=None,
                   linelist=None, fwidth=4,
-                  kdsigma=1, k=1, filename=None, dcenwave=10, bounds_setter=None):
+                  kdsigma=1, k=1, filename=None, bounds_setter=None):
     """
     Find the best wavelength solution from the set of initial models.
 
@@ -676,8 +671,7 @@ def find_solution(init_models, config, peaks=None, peak_weights=None,
         matches = perform_piecewise_fit(model, peaks, arc_lines, pixel_start,
                                         kdsigma, order=config["order"],
                                         min_lines_per_fit=min_lines_per_fit,
-                                        k=k, dcenwave=dcenwave,
-                                        bounds_setter=bounds_setter)
+                                        k=k, bounds_setter=bounds_setter)
 
         # We perform a regular least-squares fit to all the matches
         # we've made. This allows a high polynomial order to be
@@ -742,7 +736,7 @@ def find_solution(init_models, config, peaks=None, peak_weights=None,
 
 def perform_piecewise_fit(model, peaks, arc_lines, pixel_start, kdsigma,
                           order=3, min_lines_per_fit=15, k=1,
-                          arc_weights=None, dcenwave=10, bounds_setter=None):
+                          arc_weights=None, bounds_setter=None):
     """
     This function performs fits in multiple regions of the 1D arc spectrum.
     Given a starting location, a suitable fitting region is "grown" outwards
@@ -788,7 +782,6 @@ def perform_piecewise_fit(model, peaks, arc_lines, pixel_start, kdsigma,
     wave_start = model(pixel_start)
     dw_start = np.diff(model([pixel_start - 0.5, pixel_start + 0.5]))[0]
     match_radius = 2 * abs(dw_start)
-    dc0 = dcenwave
     fits_to_do = [(pixel_start, wave_start, dw_start)]
 
     first = True