exotic/api/colab.py

@@ -370,7 +370,7 @@ def make_inits_file(planetary_params, image_dir, output_dir, first_image, targ_c
             "Observing Notes": "%s",
 
             "Plate Solution? (y/n)": "y",
-            "Add Comparison Stars from AAVSO? (y/n)": "n",
+            "Add Comparison Stars from AAVSO? (y/n)": "y",
 
             "Target Star X & Y Pixel": %s,
             "Comparison Star(s) X & Y Pixel": %s,

exotic/exotic.py

@@ -902,18 +902,18 @@ def simbad_query(sample):
 
 # Apply calibrations if applicable
 def apply_cals(image_data, gen_dark, gen_bias, gen_flat, i):
-    if gen_dark.size != 0:
+    if gen_dark is not None and gen_dark.size != 0:
         if i == 0:
             log_info("Dark subtracting images.")
         image_data = image_data - gen_dark
-    elif gen_bias.size != 0:  # if a dark is not available, then at least subtract off the pedestal via the bias
+    elif gen_bias is not None and gen_bias.size != 0:  # if a dark is not available, then at least subtract off the pedestal via the bias
         if i == 0:
             log_info("Bias-correcting images.")
         image_data = image_data - gen_bias
     else:
         pass
 
-    if gen_flat.size != 0:
+    if gen_flat is not None and gen_flat.size != 0:
         if i == 0:
             log_info("Flattening images.")
         gen_flat[gen_flat == 0] = 1
@@ -1272,7 +1272,9 @@ def aperPhot(data, starIndex, xc, yc, r=5, dr=5):
 
 def skybg_phot(data, starIndex, xc, yc, r=10, dr=5, ptol=99, debug=False):
     # create a crude annulus to mask out bright background pixels
-    xv, yv = mesh_box([xc, yc], np.round(r + dr))
+    # the box will not extend beyond the borders of the image
+    image_height, image_width = data.shape
+    xv, yv = mesh_box([xc, yc], np.round(r + dr), maxx=image_width, maxy=image_height)
     rv = ((xv - xc) ** 2 + (yv - yc) ** 2) ** 0.5
     mask = (rv > r) & (rv < (r + dr))
     try:
@@ -1327,6 +1329,55 @@ def skybg_phot(data, starIndex, xc, yc, r=10, dr=5, ptol=99, debug=False):
         plt.show()
     return mode(dat.flatten(), nan_policy='omit', keepdims=True).mode[0], np.nanstd(dat.flatten()), np.sum(mask)
 
+def process_dark_frames(dark_files):
+    """Process dark frames and return the master dark."""
+    if not dark_files:
+        return None
+    # Dark files whose median is much higher than the overall dark files median will be filtered
+    # e.g. to discard saturated dark files that may negatively affect the master dark used to calibrate the science frames
+    # First pass: collect all dark frame medians
+    darks_medians = [(dark_file, np.nanmedian(fits.getdata(dark_file))) for dark_file in dark_files]
+
+    d_median = np.median([median for _, median in darks_medians])
+    threshold = 1.7  # 70% higher than overall median
+
+    # Second pass: collect valid dark frames
+    darks_img_list = []
+    for dark_file, dark_median in darks_medians:
+        median_ratio = dark_median / d_median
+        if median_ratio > threshold:
+            log_info(
+                f"\nWarning: Skipping suspicious dark frame {dark_file}: "
+                f"median/overall_median = {median_ratio:.2f}\n",
+                warn=True
+            )
+            continue
+        dark_data = fits.getdata(dark_file)
+        darks_img_list.append(dark_data)
+
+    return np.median(darks_img_list, axis=0) if darks_img_list else None
+
+def process_bias_frames(bias_files):
+    """Process bias frames and return the master bias."""
+    if not bias_files:
+        return None
+
+    biases_img_list = [fits.getdata(bias_file) for bias_file in bias_files]  
+    return np.median(biases_img_list, axis=0) if biases_img_list else None
+
+def process_flat_frames(flat_files, master_bias=None):
+    """Process flat frames and return the normalized master flat."""
+    if not flat_files:
+        return None
+
+    flats_img_list = [fits.getdata(flat_file) for flat_file in flat_files]      
+    master_flat = np.median(flats_img_list, axis=0)
+    # Bias subtract after creating master flat
+    if master_bias is not None:
+        master_flat = master_flat - master_bias
+    # Normalize
+    medi = np.median(master_flat)
+    return master_flat / medi
 
 def convert_jd_to_bjd(non_bjd, p_dict, info_dict):
     try:
@@ -1905,34 +1956,9 @@ def main():
 
         if fitsortext == 1:
             # Only do the dark correction if user selects this option
-            if exotic_infoDict['darks']:
-                darksImgList = []
-                for darkFile in exotic_infoDict['darks']:
-                    darkData = fits.getdata(darkFile)
-                    darksImgList.append(darkData)
-                generalDark = np.median(darksImgList, axis=0)
-
-            if exotic_infoDict['biases']:
-                biasesImgList = []
-                for biasFile in exotic_infoDict['biases']:
-                    biasData = fits.getdata(biasFile)
-                    biasesImgList.append(biasData)
-                generalBias = np.median(biasesImgList, axis=0)
-
-            if exotic_infoDict['flats']:
-                flatsImgList = []
-                for flatFile in exotic_infoDict['flats']:
-                    flatData = fits.getdata(flatFile)
-                    flatsImgList.append(flatData)
-                notNormFlat = np.median(flatsImgList, axis=0)
-
-                # if the bias exists, bias subtract the flatfield
-                if exotic_infoDict['biases']:
-                    notNormFlat = notNormFlat - generalBias
-
-                # NORMALIZE
-                medi = np.median(notNormFlat)
-                generalFlat = notNormFlat / medi
+            generalDark = process_dark_frames(exotic_infoDict['darks'])
+            generalBias = process_bias_frames(exotic_infoDict['biases'])
+            generalFlat = process_flat_frames(exotic_infoDict['flats'], generalBias)
 
             if exotic_infoDict['demosaic_fmt']:
                 demosaic_fmt = exotic_infoDict['demosaic_fmt'].upper()

exotic/plots.py

@@ -50,6 +50,7 @@ def plot_centroids(x_targ, y_targ, x_ref, y_ref, times, target_name, save, date)
     plt.close()
 
 def plot_fov(aper, annulus, sigma, x_targ, y_targ, x_ref, y_ref, image, image_scale, targ_name, save, date, opt_method, min_aper_fov, min_annulus_fov):
+
     ref_circle, ref_circle_sky = None, None
     picframe = 10. * (aper + 15. * sigma)
 
@@ -65,12 +66,14 @@ def plot_fov(aper, annulus, sigma, x_targ, y_targ, x_ref, y_ref, image, image_sc
         else:
             outer_circle_color = 'lime'
 
-        # Create the target and reference circles
-        target_circle = plt.Circle((x_targ, y_targ), aper, color='r', fill=False, ls='-')
-        target_circle_sky = plt.Circle((x_targ, y_targ), aper + annulus, color=outer_circle_color, fill=False, ls='-')
+        # Create the target circles
+        # We are using abs(aper) to account for a negative aperture in case EXOTIC is not using a comparison star
+        target_circle = plt.Circle((x_targ, y_targ), abs(aper), color=outer_circle_color, fill=False, ls='-')
+        target_circle_sky = plt.Circle((x_targ, y_targ), abs(aper) + annulus, color=outer_circle_color, fill=False, ls='-')
 
+        # IF EXOTIC is using a comparison star, create its circles
         if aper >= 0:
-            ref_circle = plt.Circle((x_ref, y_ref), aper, color='r', fill=False, ls='-')
+            ref_circle = plt.Circle((x_ref, y_ref), aper, color=outer_circle_color, fill=False, ls='-')
             ref_circle_sky = plt.Circle((x_ref, y_ref), aper + annulus, color=outer_circle_color, fill=False, ls='-')
 
         interval = ZScaleInterval()
@@ -83,23 +86,23 @@ def plot_fov(aper, annulus, sigma, x_targ, y_targ, x_ref, y_ref, image, image_sc
 
         ax.add_artist(target_circle)
         ax.add_artist(target_circle_sky)
-        ax.text(x_targ + aper + annulus + 5, y_targ, targ_name, color='w', fontsize=10,
+        ax.text(x_targ + abs(aper) + annulus + 5, y_targ, targ_name, color='w', fontsize=10,
                 path_effects=[path_effects.withStroke(linewidth=2, foreground='black')])
 
-        if aper >= 0:
+        if aper >= 0: #EXOTIC is using a comparison star
             ax.add_artist(ref_circle)
             ax.add_artist(ref_circle_sky)
             ax.text(x_ref + aper + annulus + 5, y_ref, 'Comp Star', color='w', fontsize=10,
                     path_effects=[path_effects.withStroke(linewidth=2, foreground='black')])
 
         handles = []
-        label_aper = f"{opt_method} Photometry\n(Min Aper: {min_aper_fov:.2f} px)\n(Min Annulus: {min_annulus_fov:.2f} px)"
+        label_aper = f"{opt_method} Photometry\n(Min Aper: {abs(min_aper_fov):.2f} px)\n(Min Annulus: {min_annulus_fov:.2f} px)"
 
         if opt_method == "Aperture":
-            aperture_line = Line2D([], [], color='r', linestyle='-', label=label_aper)
+            aperture_line = Line2D([], [], color=outer_circle_color, linestyle='-', label=label_aper)
             handles.append(aperture_line)
         elif opt_method == "PSF":
-            psf_line = Line2D([], [], color='lime', linestyle='-', label=label_aper)
+            psf_line = Line2D([], [], color=outer_circle_color, linestyle='-', label=label_aper)
             handles.append(psf_line)
 
         plt.title(f"FOV for {targ_name}\n({image_scale})")

exotic/version.py

@@ -1 +1 @@
-__version__ = '4.2.3'
+__version__ = '4.3.1'

requirements.txt

@@ -10,11 +10,12 @@ imreg_dft~=2.0.0
 LDTk~=1.8.4
 lmfit~=1.3.2
 matplotlib~=3.9.2
+numba~=0.59.1
 numpy~=1.26.4
 pandas~=2.2.3
 panel~=1.5.2
 photutils~=2.0.0
-pylightcurve~=4.0.1
+pylightcurve>=4.0.1,<5
 python_dateutil~=2.9
 python_version
 pyvo~=1.5.3