Adding orderStacks as a default coadd.py output for Echelle spectrographs

pypeit/core/coadd.py

@@ -2401,7 +2401,7 @@ def ech_combspec(waves_arr_setup, fluxes_arr_setup, ivars_arr_setup, gpms_arr_se
         nbests=None, which will just use one fourth of the orders for a given
         setup.
     wave_method : str, optional
-        method for generating new wavelength grid with get_wave_grid. Deafult is
+        method for generating new wavelength grid with get_wave_grid. Default is
         'log10' which creates a uniformly space grid in log10(lambda), which is
         typically the best for echelle spectrographs
     dwave : float, optional
@@ -2557,6 +2557,8 @@ def ech_combspec(waves_arr_setup, fluxes_arr_setup, ivars_arr_setup, gpms_arr_se
     # data shape
     nsetups=len(waves_arr_setup)
 
+    msgs.info(f'Number of setups to cycle through is: {nsetups}')
+
     if setup_ids is None:
         setup_ids = list(string.ascii_uppercase[:nsetups])
 
@@ -2605,7 +2607,7 @@ def ech_combspec(waves_arr_setup, fluxes_arr_setup, ivars_arr_setup, gpms_arr_se
                                             wave_grid_min=wave_grid_min,
                                             wave_grid_max=wave_grid_max, dwave=dwave, dv=dv,
                                             dloglam=dloglam, spec_samp_fact=spec_samp_fact)
-
+    msgs.info(f'The shape of the giant wave grid here is: {np.shape(wave_grid)}')
     # Evaluate the sn_weights. This is done once at the beginning
     weights = []
     rms_sn_setup_list = []
@@ -2733,14 +2735,28 @@ def ech_combspec(waves_arr_setup, fluxes_arr_setup, ivars_arr_setup, gpms_arr_se
         for iord in range(norders[isetup]):
             ind_start = iord*nexps[isetup]
             ind_end = (iord+1)*nexps[isetup]
+            wave_grid_ord = waves_setup_list[isetup][ind_start]
+            # if the wavelength grid is non-monotonic, resample onto a loglam grid
+            wave_grid_diff_ord = np.diff(wave_grid_ord)
+            if np.any(wave_grid_diff_ord < 0):
+                msgs.warn(f'This order ({iord}) has a non-monotonic wavelength solution. Resampling now: ')
+                wave_grid_ord = np.linspace(np.min(wave_grid_ord), np.max(wave_grid_ord), len(wave_grid_ord))
+                wave_grid_diff_ord = np.diff(wave_grid_ord)
+
+            wave_grid_diff_ord = np.append(wave_grid_diff_ord, wave_grid_diff_ord[-1])
+            wave_grid_mid_ord = wave_grid_ord + wave_grid_diff_ord / 2.0
+            # removing the last bin since the midpoint now falls outside of wave_grid rightmost bin. This matches
+            # the convention in wavegrid above
+            wave_grid_mid_ord = wave_grid_mid_ord[:-1]
+            #trim off first and last pixel in case of edge effects in wavelength calibration
             wave_order_stack_iord, flux_order_stack_iord, ivar_order_stack_iord, gpm_order_stack_iord, \
                 nused_order_stack_iord, outgpms_order_stack_iord = spec_reject_comb(
-                wave_grid, wave_grid_mid, waves_setup_list[isetup][ind_start:ind_end],
+                wave_grid_ord, wave_grid_mid_ord, waves_setup_list[isetup][ind_start:ind_end],
                 fluxes_scale_setup_list[isetup][ind_start:ind_end], ivars_scale_setup_list[isetup][ind_start:ind_end],
                 gpms_setup_list[isetup][ind_start:ind_end], weights_setup_list[isetup][ind_start:ind_end],
                 sn_clip=sn_clip, lower=lower, upper=upper, maxrej=maxrej, maxiter_reject=maxiter_reject, debug=debug_order_stack,
                 title='order_stacks')
-            waves_order_stack.append(wave_order_stack_iord)
+            waves_order_stack.append(wave_grid_mid_ord)
             fluxes_order_stack.append(flux_order_stack_iord)
             ivars_order_stack.append(ivar_order_stack_iord)
             gpms_order_stack.append(gpm_order_stack_iord)

pypeit/core/flux_calib.py

@@ -31,6 +31,7 @@
 from pypeit import dataPaths
 
 
+
 # TODO: Put these in the relevant functions
 TINY = 1e-15
 SN2_MAX = (20.0) ** 2

pypeit/core/telluric.py

@@ -774,6 +774,7 @@ def general_spec_reader(specfile, ret_flam=False, chk_version=False, ret_order_s
             raise ValueError("This is an ugly hack until the DataContainer bug is fixed")
         head = sobjs.header
         wave, counts, counts_ivar, counts_gpm = unpack_orders(sobjs, ret_flam=ret_flam)
+        wave_grid_mid = None
         # Made a change to the if statement to account for unpack_orders now squeezing returned arrays
         #if (head['PYPELINE'] !='Echelle') and (wave.shape[1]>1)
         if (head['PYPELINE'] !='Echelle') and (wave.ndim>1):

pypeit/core/wavecal/autoid.py

@@ -1111,8 +1111,8 @@ def full_template(spec, lamps, par, ok_mask, det, binspectral, nsnippet=2, slit_
         lines_wav = template_dict['lines_wav'] 
         lines_fit_ord = template_dict['lines_fit_ord']
 
-    temp_wv = temp_wv_og
-    temp_spec = temp_spec_og
+    temp_wv = np.copy(temp_wv_og)
+    temp_spec = np.copy(temp_spec_og)
 
     # Deal with binning (not yet tested)
     if binspectral != temp_bin:
@@ -1145,7 +1145,7 @@ def full_template(spec, lamps, par, ok_mask, det, binspectral, nsnippet=2, slit_
         obs_spec_i = spec[:,slit]
         # get FWHM for this slit
         fwhm = set_fwhm(par, measured_fwhm=measured_fwhms[slit], verbose=True)
-
+        
         # Find the shift
         ncomb = temp_spec.size
         # Remove the continuum before adding the padding to obs_spec_i

pypeit/core/wavecal/wvutils.py

@@ -285,6 +285,9 @@ def get_wave_grid(waves=None, gpms=None, wave_method='linear', iref=0, wave_grid
             wave_grid = np.power(10.0,newloglam)
 
         elif wave_method == 'iref': # Use the iref index wavelength array
+            msgs.info(f'iref for the list is set to {iref}')
+            msgs.info(f'The shape of the list is: {np.shape(waves)}')
+            msgs.info(f'shape of the first wave_grid in the list is: {np.shape(waves[iref])}')
             wave_tmp = waves[iref]
             wave_grid = wave_tmp[wave_tmp > 1.0]
             if spec_samp_fact != 1: # adjust sampling via internal interpolation

pypeit/sensfunc.py

@@ -741,6 +741,8 @@ def write_QA(self):
             wave_gpm = self.sens['SENS_FLUXED_STD_WAVE'][iorddet] > 1.0
             model_flux_sav[iorddet][wave_gpm] = model_interp_func(self.sens['SENS_FLUXED_STD_WAVE'][iorddet][wave_gpm])
 
+        self.sens['SENS_STD_MODEL_FLAM'] = model_flux_sav
+