diff --git a/tutorials/euclid/euclid_clusters_tutorial.md b/tutorials/euclid/euclid_clusters_tutorial.md
index 81652b95..d6f1f905 100644
--- a/tutorials/euclid/euclid_clusters_tutorial.md
+++ b/tutorials/euclid/euclid_clusters_tutorial.md
@@ -188,6 +188,7 @@ To explore a different control field, change the seed value or remove the seed e
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def check_mer_tile_requirement(coord, search_radius=2.0):
     """Check whether a sky coordinate is covered by exactly 4 Euclid MER science images.
@@ -230,6 +231,7 @@ def check_mer_tile_requirement(coord, search_radius=2.0):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def find_control_field(cluster_df, cluster_ra, cluster_dec, min_distance_arcmin=15, max_attempts=100):
     """Find a random control field offset from the known cluster catalog.
@@ -342,6 +344,7 @@ Because the four photometric bands are independent of each other, `download_and_
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def _download_band(band, mer_images, field_coord, field_id, cache_dir, s3):
     """Download one photometric band from S3 and write a cutout FITS to the local cache.
@@ -401,6 +404,7 @@ def _download_band(band, mer_images, field_coord, field_id, cache_dir, s3):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def download_and_cache_field(mer_images, field_name, field_coord, field_id):
     """Stream cutout FITS images from S3 and cache them locally.
@@ -482,6 +486,7 @@ Displaying the cluster and control fields side by side at the same stretch gives
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def normalize_with_consistent_stretch(cluster_cutouts, control_cutouts, lower_percentile=1, upper_percentile=99):
     """Normalize cluster and control cutouts using a shared percentile stretch.
@@ -538,6 +543,7 @@ def normalize_with_consistent_stretch(cluster_cutouts, control_cutouts, lower_pe
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def downsample(arr, factor=4):
     """Block-average a 2-D or 3-D image array by an integer factor for faster display.
@@ -646,6 +652,7 @@ Then we show part of the cluster dataframe to see what information we have avail
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def query_galaxies_for_field(ra, dec, field_name, redshift_center, redshift_width=0.1):
     """Query galaxies within a redshift slice around a sky position.
@@ -781,6 +788,7 @@ The function prints the number of galaxies within bounds as a diagnostic, if a l
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def apply_dbscan_clustering(galaxy_df, wcs, rgb_image, field_name, eps=500, min_samples=18):
     """Apply DBSCAN to detect galaxy overdensities in a redshift-selected sample. DBSCAN operates on the 2-D projected
@@ -938,6 +946,7 @@ Combining field galaxies from both the cluster and control fields gives us a lar
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def identify_cluster_members(galaxy_df, labels, galaxy_coords, field_name):
     """Separate a galaxy sample into cluster members and field galaxies.
@@ -1058,6 +1067,7 @@ We convert the uniform-aperture fluxes in the photo-z catalog to AB magnitudes a
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def calculate_color_magnitude(df):
     """Convert uniform-aperture fluxes to AB magnitudes and compute Y-H color.
@@ -1102,6 +1112,7 @@ field_cmd = calculate_color_magnitude(all_field_galaxies[['flux_y_unif', 'flux_h
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def remove_outliers_bounds(df, h_min=17, h_max=25, yh_min=-0.5, yh_max=1.5):
     """Filter a colour-magnitude table to a physically motivated range.
@@ -1242,6 +1253,7 @@ Ten spectra per group is sufficient to show whether the cluster and field popula
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def get_n_spectra(obj_ids, n=10):
     """
@@ -1434,6 +1446,7 @@ The three helper functions below handle this: `preprocess_spectrum` continuum-su
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def preprocess_spectrum(spec):
     """Continuum-remove + robust-normalize a spectrum for visualization."""
@@ -1487,6 +1500,7 @@ def preprocess_spectrum(spec):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def build_stack(spectra_dict, w_grid):
     """Interpolate processed spectra onto a common grid and return matrix [nobj, ngrid]."""
@@ -1514,6 +1528,7 @@ def build_stack(spectra_dict, w_grid):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def lines_in_window(wmin, wmax, z):
     """Return dict of lines whose observed wavelength falls in [wmin, wmax]."""
@@ -1553,6 +1568,7 @@ field_stack   = build_stack(field_spectra,   w_grid)
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def label_emission_line(ax, x, label, y=0.92, rotation=90, color='k'):
     """Place a vertical emission-line label at wavelength x."""
@@ -1639,6 +1655,7 @@ This provides an independent check of cluster membership using spectroscopic red
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def search_ned_field(ra, dec, z_min, z_max, label, radius_arcmin=3, max_retries=3):
     """Search NED within radius_arcmin of (ra, dec), filter to z_min–z_max, and print a summary.
@@ -1731,6 +1748,7 @@ We now overlay the NED-catalogued objects on the cluster and control field image
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def overlay_ned_sources(
     cluster_objects, control_objects,
diff --git a/tutorials/euclid/merged-objects-hats-catalog/4-euclid-q1-hats-magnitudes.md b/tutorials/euclid/merged-objects-hats-catalog/4-euclid-q1-hats-magnitudes.md
index 90bbc33b..3acb3f0b 100644
--- a/tutorials/euclid/merged-objects-hats-catalog/4-euclid-q1-hats-magnitudes.md
+++ b/tutorials/euclid/merged-objects-hats-catalog/4-euclid-q1-hats-magnitudes.md
@@ -232,6 +232,7 @@ Since the template-fit photometry is recommended for extended objects, we'll sep
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 # Galaxy + any. Star + galaxy. QSO + galaxy.
 classes = {"Galaxy": (2, 3, 6, 7), "Star": (1, 3), "QSO": (4, 6)}
@@ -304,6 +305,7 @@ This figure is inspired by Romelli Fig. 6 (top panel).
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 # Only consider objects within these mag and mag difference limits.
 mag_limits, mag_diff_limits = (16, 24), (-1, 1)
diff --git a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md
index b2f1f94a..00977df2 100644
--- a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md
+++ b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_SED_fit.md
@@ -110,6 +110,7 @@ list out some basic information including column names
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def inspect_parquet_columns(s3_path, *, region='us-east-1', max_rows=0):
     """
@@ -176,6 +177,7 @@ Next we merge the SN sample with the host galaxy fluxes and physical properties.
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def assemble_SN_data(sn_flux_file, galaxy_flux_file, galaxy_info_file, *, region='us-east-1'):
     """
@@ -251,6 +253,7 @@ A tight correlation along the 1:1 line indicates successful matching.
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_redshift_comparison(df):
     """
@@ -349,6 +352,7 @@ roman_bands = {
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def scale_rubin_to_roman(df, rubin_bands, roman_bands, match_wave=0.87, verbose=False):
     """
@@ -445,6 +449,7 @@ def scale_rubin_to_roman(df, rubin_bands, roman_bands, match_wave=0.87, verbose=
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def photons_to_maggies_with_filter(photon_flux, filt):
     """
@@ -493,6 +498,7 @@ def photons_to_maggies_with_filter(photon_flux, filt):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def add_maggie_columns(df):
     """
@@ -545,6 +551,7 @@ def add_maggie_columns(df):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_single_sed(df, rubin_bands, roman_bands, galaxy_index=0, loglog=False):
     """
@@ -636,6 +643,7 @@ plot_single_sed(df_scaled, rubin_bands, roman_bands, 1, loglog=True)
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_many_seds(df, rubin_bands, roman_bands, n_galaxies=10, loglog=False):
     """
@@ -744,6 +752,7 @@ Prospector has its own built in functions to do writing and reading, but
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def df_to_all_obs(df, flux_err_frac=0.1):
     """
@@ -801,6 +810,7 @@ def df_to_all_obs(df, flux_err_frac=0.1):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def process_output(output, model, obs, sps):
     """
@@ -903,6 +913,7 @@ def process_output(output, model, obs, sps):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def run_fit(obs, model_params, sps, lnprobfn, noise_model, fitting_kwargs):
     """
@@ -965,6 +976,7 @@ def run_fit(obs, model_params, sps, lnprobfn, noise_model, fitting_kwargs):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def run_fit_star(args):
     """
@@ -1007,6 +1019,7 @@ def run_fit_star(args):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def fit_sample_parallel(obs_list, model_params, sps, lnprobfn, noise_model=(None,None),
                         fitting_kwargs=None, nproc=None):
@@ -1057,6 +1070,7 @@ def fit_sample_parallel(obs_list, model_params, sps, lnprobfn, noise_model=(None
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def save_outputs_hdf5(outputs, obs_list, filename="all_galaxies.h5"):
     """
@@ -1240,6 +1254,7 @@ use this section to quickly load the fitted results instead of re-running Prospe
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def download_results_hdf5_from_gdrive(file_id, dest_path="./sn_fits.h5"):
     """
@@ -1272,6 +1287,7 @@ def download_results_hdf5_from_gdrive(file_id, dest_path="./sn_fits.h5"):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def load_outputs_hdf5(filename="sn_fits.h5",
                       file_id="1BuymtEXJsxbd8PqwkIEmff-76JGGa650"):
@@ -1353,6 +1369,7 @@ After fitting, we can visualize the results by comparing observed and modeled SE
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_SED_fit(outputs, obs_list, galaxy_id):
     """
@@ -1464,6 +1481,7 @@ for gid in galaxy_ids[:5]:
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_corner(outputs, obs_list, galaxy_id):
     """
@@ -1556,6 +1574,7 @@ Agreement along the 1:1 line indicates successful recovery of physical parameter
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_stellar_mass_verification(df, outputs):
     """
@@ -1642,6 +1661,7 @@ We compare the distribution of Prospector-derived stellar masses for the two SN
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_1a_vs_cc_mass(df_sn, outputs):
     """
diff --git a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md
index 3fc43b21..bd18bf99 100644
--- a/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md
+++ b/tutorials/simulated-data/OpenUniverse2024/openuniverse2024_TDE_light_curve.md
@@ -143,6 +143,7 @@ Let's take a look at a few images to see what we are dealing with.
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def show_gallery(files, max_images=9):
     """
@@ -292,6 +293,7 @@ OU_RUBIN_SIA_COLLECTION = 'simulated_rubin_openuniverse2024'
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def get_s3_fpath(cloud_access):
     cloud_info = json.loads(cloud_access) # converts str to dict
@@ -307,6 +309,7 @@ First, we find the filenames of the images in the Roman TDS survey which include
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def Roman_TDS_image_search(host_galaxy, radius, bandname):
     """
@@ -360,6 +363,7 @@ Since there are > 100 TDS images per band, we will:
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def select_images_by_mjd_quantiles(images, n_select=9):
     """
@@ -426,6 +430,7 @@ The two plotting functions then compile these measurements into time-ordered plo
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def run_aperture_photometry(host_galaxy, bandname, image_column="image_filenames", aperture_radius=1.0):
     """
@@ -525,6 +530,7 @@ host_galaxy
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_single_band_light_curve(df, bandname, start_mjd):
     """
@@ -595,6 +601,7 @@ def plot_single_band_light_curve(df, bandname, start_mjd):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_multiband_light_curve(df, bands, start_mjd):
     """
@@ -685,6 +692,7 @@ We follow the example in this [tutorial](https://caltech-ipac.github.io/irsa-tut
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def make_cutout(fname, ra, dec, size=100):
     """
@@ -739,6 +747,7 @@ def make_cutout(fname, ra, dec, size=100):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def cutout_gallery(image_filenames, mjd_list, ra, dec, aperture_radius_pix_list, size=100, ncols=4,
                    galaxy_id=None):
diff --git a/tutorials/simulated-data/roman_hlss_number_density.md b/tutorials/simulated-data/roman_hlss_number_density.md
index a727fc46..6734cbb7 100644
--- a/tutorials/simulated-data/roman_hlss_number_density.md
+++ b/tutorials/simulated-data/roman_hlss_number_density.md
@@ -80,6 +80,7 @@ import requests
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def download_file(file_url, save_path):
     """
@@ -118,6 +119,7 @@ def download_file(file_url, save_path):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def download_files_in_parallel(file_url_list, save_dir):
     """
@@ -154,6 +156,7 @@ def download_files_in_parallel(file_url_list, save_dir):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def download_simulations(download_dir, download_all=True):
     """
@@ -236,6 +239,7 @@ download_simulations(download_dir, download_all=False)
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def fetch_column_names():
     """
@@ -274,6 +278,7 @@ def fetch_column_names():
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def read_hdf5_to_pandas(file_path, columns_to_keep, columns_to_convert):
     """
@@ -335,6 +340,7 @@ def read_hdf5_to_pandas(file_path, columns_to_keep, columns_to_convert):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def assign_redshift_bins(df, z_min, z_max, dz=0.1):
     """
@@ -464,6 +470,7 @@ Do the numbers we are getting and their distributions make sense?
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def get_bin_area(dec_edges, ra_width):
     """
@@ -492,6 +499,7 @@ def get_bin_area(dec_edges, ra_width):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def calculate_counts_per_deg2(df):
     """
@@ -545,6 +553,7 @@ def calculate_counts_per_deg2(df):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_dec_bins(dec_bins, galaxy_counts):
     """
@@ -583,6 +592,7 @@ def plot_dec_bins(dec_bins, galaxy_counts):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_dec_bins_per_deg2(dec_bins, counts_per_deg2):
     """
@@ -635,6 +645,7 @@ This is a first pass at making a number density plot with just a single data fil
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def jackknife_galaxy_count(df, grid_cells=10):
     """
@@ -715,6 +726,7 @@ def jackknife_galaxy_count(df, grid_cells=10):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_galaxy_counts_vs_redshift_with_jackknife(counts, std_dev_counts, z_bin_centers):
     """
@@ -797,6 +809,7 @@ Originally this section attempted to use dask instead of pandas to hold the data
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def calculate_area(df):
     """
@@ -830,6 +843,7 @@ def calculate_area(df):
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def galaxy_counts_per_hdf5_binned(
     df,
@@ -903,6 +917,7 @@ def galaxy_counts_per_hdf5_binned(
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_binned_galaxy_counts_vs_redshift_with_jackknife(
     counts_summary, z_bin_centers, halpha_flux_thresholds, halpha_summary,
@@ -995,6 +1010,7 @@ def plot_binned_galaxy_counts_vs_redshift_with_jackknife(
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def jackknife_wrapper(file_list,  halpha_flux_thresholds, find_halpha_bins=False):
     """
@@ -1151,6 +1167,7 @@ The uncertainties are plotted as error bars on the data points above.  There are
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def plot_jackknife_fractional_uncertainty(z_bin_centers, counts, std_dev_counts):
     """
diff --git a/tutorials/spherex/spherex_mosaic.md b/tutorials/spherex/spherex_mosaic.md
index 65f9c382..26f61c63 100644
--- a/tutorials/spherex/spherex_mosaic.md
+++ b/tutorials/spherex/spherex_mosaic.md
@@ -81,6 +81,7 @@ import matplotlib as mpl
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 ## Define some plotting formats
 def load_plotting_defaults():
@@ -108,10 +109,10 @@ We first show how to open the SPHEREx mosaic cube and how to assign wavelengths
 
 The mosaic cube was obtained from the [IRSA SPHEREx mosaic tool](https://irsa.ipac.caltech.edu/applications/spherex/tool-mosaic) GUI. The mosaic tool computes a cube (x,y,$\lambda$) from SPHEREx data at a given sky position provided by the user. In brief, the tool gathers all the SPHEREx spectral images at that position, extracts the pixels at similar wavelenghts, and reprojects them into cube layers at the respective wavelengths. The final cube has 102 wavelength layers from $0.75$ to $5\,{\rm \mu m}$. The user can also specify the spatial size of the cube.
 
-Because the IRSA SPHEREx mosaic tool does not have an API, we have created and downloaded the mosaic already and added it to the `./data/` directory. 
+Because the IRSA SPHEREx mosaic tool does not have an API, we have created and downloaded the mosaic already and added it to the `./data/` directory.
 
 ```{tip}
-To retrieve the same mosaic as provided here, go to the [IRSA mosaic tool](https://irsa.ipac.caltech.edu/applications/spherex/tool-mosaic) and type in *M101* in the `Output Mosaic Center` text field. For the size of the mosaic, choose 30 arcminutes for both axis. For the output mosaic pixel scale choose 9 arcseconds. 
+To retrieve the same mosaic as provided here, go to the [IRSA mosaic tool](https://irsa.ipac.caltech.edu/applications/spherex/tool-mosaic) and type in *M101* in the `Output Mosaic Center` text field. For the size of the mosaic, choose 30 arcminutes for both axis. For the output mosaic pixel scale choose 9 arcseconds.
 ```
 
 We first define the path to the mosaic cube.
@@ -152,7 +153,7 @@ with fits.open(fn_spherex) as hdul:
     # adjust image units and update header BUNIT
     cube_img = hdul["IMAGE"].data * 23.5045 * (pixscale_mosaic)**2 / 1e3 # Mjy/sr -> mjy
     cube_hdr['BUNIT'] = "mjy"
-    
+
     # extract the 2D spatial WCS (dropping the spectral axis)
     cube_wcs = WCS(cube_hdr, fobj=hdul).celestial
     print(f"Loaded SPHEREx cube with pixel scale {pixscale_mosaic} arcsec/px")
@@ -202,13 +203,14 @@ Furthermore, the function creates a plot showing the collapsed cube and the aper
 
 
 ```{tip}
-This function provides a very simple aperture photometry with background subtraction. The method can be extended. For example the function does not calculate photometric errors.  
+This function provides a very simple aperture photometry with background subtraction. The method can be extended. For example the function does not calculate photometric errors.
 ```
 
 ```{code-cell} ipython3
 ---
 jupyter:
   source_hidden: true
+tags: [hide-cell]
 ---
 def measure_aperture_photometry_cube(cube, *, r_aperture_px=20, r_inner_px=60, r_outer_px=100, makeplot=True):
     '''
@@ -233,7 +235,7 @@ def measure_aperture_photometry_cube(cube, *, r_aperture_px=20, r_inner_px=60, r
     -------
     astropy.QTable
         A table including the photometry results (sum of aperture flux, plane ID, etc).
-    
+
     '''
 
     ## Define helper function to compute the photometry efficiently
@@ -252,34 +254,34 @@ def measure_aperture_photometry_cube(cube, *, r_aperture_px=20, r_inner_px=60, r
             Aperture for photometry extraction.
         annulus_aperture : photutils.CircularAnnulus
             Annulus aperture for background estimation.
-        
+
         Returns
         -------
         astropy.QTable
             A table including the photometry results (sum of aperture flux, plane ID, etc).
-        
+
         '''
-        
+
         ## Calculate background
         sigclip = SigmaClip(sigma=3.0, maxiters=10)
         aperstats = ApertureStats(img, annulus_aperture, sigma_clip=sigclip)
         bkg_mean = aperstats.mean
         aperture_area = aperture.area_overlap(img)
         total_bkg = bkg_mean * aperture_area
-        
+
         ## Get photometry and subtract background
         phot_table = aperture_photometry(img, aperture)
         phot_table["aperture_sum_bkgsub"] = phot_table["aperture_sum"] - total_bkg
-    
+
         return(phot_table)
 
-    
+
     ## Define position and apertures
     position = (cube.shape[1]//2,cube.shape[2]//2)
     aperture = CircularAperture(position, r=r_aperture_px)
     annulus_aperture = CircularAnnulus(position, r_in=r_inner_px, r_out=r_outer_px)
 
-    
+
     ## Compute photometry (use helper function to iterate over planes)
     phot_all = vstack( [apphot_helper(cube[ii,:,:] , position, aperture, annulus_aperture) for ii in range(cube.shape[0])] ) # run all planes
     phot_all["id"] = np.arange(cube.shape[0])+1 ## add plane numbers back
@@ -442,7 +444,7 @@ def make_map(cube,
     -------
     numpy.ndarray
         Two-dimensional map.
-    
+
     '''
     img_map = np.nansum( cube[np.asarray(planes_feature)-1 , :,:], axis=0 ) - np.nanmedian( cube[np.asarray(planes_continuum)-1 , :,:], axis=0 )
     return(img_map)