Merge pull request #38 from radicamc/dev

dev

changelog.md

@@ -1,6 +1,11 @@
 # Changelog
 All notable changes to this project will be documented in this file.
 
+### [1.3.1] -- 2024-01-19
+#### Added
+- Fixed implementation of STScI up-the-ramp JumpStep which would not run consistently.
+- Made it easier to select stellar models to use for limb-darkening calculation in Stage 4. 
+
 ### [1.3.0] -- 2023-12-21
 #### Added
 - Compatibility with v1.12.5 of the STScI jwst package.

setup.py

@@ -4,7 +4,7 @@
 from setuptools import setup
 
 setup(name='supreme_spoon',
-      version='1.3.0',
+      version='1.3.1',
       license='MIT',
       author='Michael Radica',
       author_email='[email protected]',

supreme_spoon/fit_lightcurves.py

@@ -207,7 +207,8 @@
             assert np.all(np.array([m_h, logg, teff]) != None), msg
             c1, c2 = stage4.gen_ld_coefs(config['spectrace_ref'], wave_low,
                                          wave_up, order, m_h, logg, teff,
-                                         config['ld_data_path'])
+                                         config['ld_data_path'],
+                                         model_type=config['ld_model_type'])
             q1, q2 = juliet.reverse_q_coeffs('quadratic', c1, c2)
             # Save calculated coefficients.
             target = fits.getheader(config['infile'], 0)['TARGET']

supreme_spoon/fit_lightcurves.yaml

@@ -63,6 +63,9 @@ teff : None
 # Path to ExoTiC-LD Reference Files
 # See ExoTiC-LD documentation for more info: https://exotic-ld.readthedocs.io/en/latest/views/installation.html.
 ld_data_path :
+# Stellar models to use for limb-darkening calculation.
+# See ExoTiC-LD documentation for more info.
+ld_model_type : 'stagger'
 # Path to JWST spectrace Reference File
 # Will be in crds_cache with file name like jwst_niriss_spectrace_XXXX.fits
 spectrace_ref : './crds_cache/references/jwst/niriss/jwst_niriss_spectrace_0023.fits'

supreme_spoon/stage1.py

@@ -355,8 +355,15 @@ def run(self, save_results=True, force_redo=False, do_plot=False,
                         '.fits', '_2.pdf')
                 else:
                     plot_file1, plot_file2 = None, None
+
+                # For scale-chromatic correction, collapse the 2D timeseries to
+                # 1D for plotting purposes.
+                if self.method == 'scale-chromatic':
+                    this_ts = np.nanmedian(self.timeseries, axis=1)
+                else:
+                    this_ts = self.timeseries
                 plotting.make_oneoverf_plot(results,
-                                            timeseries=self.timeseries,
+                                            timeseries=this_ts,
                                             baseline_ints=self.baseline_ints,
                                             outfile=plot_file1,
                                             show_plot=show_plot)
@@ -366,7 +373,7 @@ def run(self, save_results=True, force_redo=False, do_plot=False,
                 else:
                     window = False
                 plotting.make_oneoverf_psd(results, self.datafiles,
-                                           timeseries=self.timeseries,
+                                           timeseries=this_ts,
                                            baseline_ints=self.baseline_ints,
                                            pixel_masks=self.pixel_masks,
                                            outfile=plot_file2,
@@ -465,16 +472,24 @@ def run(self, save_results=True, force_redo=False, flag_up_ramp=False,
                 # Get number of groups in the observation - ngroup=2 must be
                 # treated in a special way as the default pipeline JumpStep
                 # will fail.
+                # Also need to set minimum_sigclip_groups to something >nints,
+                # else the up-the-ramp jump detection will be replaced by a
+                # time-domain sigma clipping.
                 testfile = datamodels.open(self.datafiles[0])
                 ngroups = testfile.meta.exposure.ngroups
+                i_start = testfile.meta.exposure.integration_start
+                i_end = testfile.meta.exposure.integration_end
+                nints = i_end - i_start
                 testfile.close()
                 # For ngroup > 2, default JumpStep can be used.
                 if ngroups > 2 and flag_up_ramp is True:
                     step = calwebb_detector1.jump_step.JumpStep()
                     res = step.call(segment, output_dir=self.output_dir,
                                     save_results=save_results,
                                     rejection_threshold=rejection_threshold,
-                                    maximum_cores='quarter', **kwargs)
+                                    maximum_cores='quarter',
+                                    minimum_sigclip_groups=nints+100,
+                                    **kwargs)
                 # Time domain jump step must be run for ngroup=2.
                 else:
                     res = segment

supreme_spoon/stage4.py

@@ -344,8 +344,8 @@ def fit_lightcurves(data_dict, prior_dict, order, output_dir, fit_suffix,
     return results
 
 
-def gen_ld_coefs(spectrace_ref, wavebin_low, wavebin_up, order, m_h, logg, teff,
-                 ld_data_path, model_type='stagger'):
+def gen_ld_coefs(spectrace_ref, wavebin_low, wavebin_up, order, m_h, logg,
+                 teff, ld_data_path, model_type='stagger'):
     """Generate estimates of quadratic limb-darkening coefficients using the
     ExoTiC-LD package.
 

supreme_spoon/utils.py

@@ -310,10 +310,10 @@ def get_filename_root(datafiles):
     # Now assuming everything is in chronological order, just increment the
     # segment number.
     split = fileroot.split('seg')
-    for segment in range(seg_start+1, len(datafiles)+1):
+    for segment in range(seg_start+1, seg_start+len(datafiles)):
         if segment < 10:
             seg_no = 'seg00{}'.format(segment)
-        elif 10 <= segment < 99:
+        elif 10 <= segment <= 99:
             seg_no = 'seg0{}'.format(segment)
         else:
             seg_no = 'seg{}'.format(segment)
@@ -709,6 +709,7 @@ def outlier_resistant_variance(data):
     return var
 
 
+# TODO: Include something to check that inputs are correct format/not nonsensical.
 def parse_config(config_file):
     """Parse a yaml config file.
 
@@ -793,11 +794,11 @@ def save_extracted_spectra(filename, wl1, wu1, f1, e1, wl2, wu2, f2, e2, t,
     hdu3 = fits.ImageHDU(wu1, header=hdr)
     hdr = fits.Header()
     hdr['EXTNAME'] = "Flux O1"
-    hdr['UNITS'] = "Electrons"
+    hdr['UNITS'] = "DN/s"
     hdu4 = fits.ImageHDU(f1, header=hdr)
     hdr = fits.Header()
     hdr['EXTNAME'] = "Flux Err O1"
-    hdr['UNITS'] = "Electrons"
+    hdr['UNITS'] = "DN/s"
     hdu5 = fits.ImageHDU(e1, header=hdr)
 
     # Pack order 2 values.
@@ -811,17 +812,17 @@ def save_extracted_spectra(filename, wl1, wu1, f1, e1, wl2, wu2, f2, e2, t,
     hdu7 = fits.ImageHDU(wu2, header=hdr)
     hdr = fits.Header()
     hdr['EXTNAME'] = "Flux O2"
-    hdr['UNITS'] = "Electrons"
+    hdr['UNITS'] = "DN/s"
     hdu8 = fits.ImageHDU(f2, header=hdr)
     hdr = fits.Header()
     hdr['EXTNAME'] = "Flux Err O2"
-    hdr['UNITS'] = "Electrons"
+    hdr['UNITS'] = "DN/s"
     hdu9 = fits.ImageHDU(e2, header=hdr)
 
     # Pack time axis.
     hdr = fits.Header()
     hdr['EXTNAME'] = "Time"
-    hdr['UNITS'] = "BJD"
+    hdr['UNITS'] = "MJD"
     hdu10 = fits.ImageHDU(t, header=hdr)
 
     if save_results is True: