diff --git a/skyreader/constants.py b/skyreader/constants.py
index 37c79ef..e7b9463 100644
--- a/skyreader/constants.py
+++ b/skyreader/constants.py
@@ -1,5 +1,10 @@
 # LAT 14-year Source Catalog (4FGL-DR4)
 # fermi.gsfc.nasa.gov/ssc/data/access/lat
-CATALOG_NAME = "gll_psc_v35.fit"
+# CATALOG_NAME = "gll_psc_v35.fit"
+
+# LAT 16-year Source List (FL16Y)
+# https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl16y/
+CATALOG_NAME = "gll_psc_v41.fit"
+
 BASE_DIR = "/cvmfs/icecube.opensciencegrid.org/users/"
 CATALOG_PATH = f"{BASE_DIR}azegarelli/realtime/catalogs/{CATALOG_NAME}"
diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py
index 51dedda..d65a78c 100644
--- a/skyreader/plot/plot.py
+++ b/skyreader/plot/plot.py
@@ -26,7 +26,7 @@
     plot_catalog
 )
 
-# LAT 14-year Source Catalog (4FGL-DR4 in FITS format) ; https://fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/
+# Fermi catalog
 from skyreader.constants import CATALOG_PATH
 
 from ..utils.areas import calculate_area, get_contour_areas
@@ -276,7 +276,7 @@ def create_plot_zoomed(
         extra_radius=np.nan,
         systematics=False,
         plot_bounding_box=False,
-        plot_4fgl=False,
+        plot_fermi_sources=False,
         circular=False,
         circular_err50=0.2,
         circular_err90=0.7,
@@ -544,8 +544,8 @@ def bounding_box(ra, dec, theta, phi):
             rot=(lon, lat, 0),
             bounds=(lower_lon, upper_lon, lower_lat, upper_lat)
         )
-        if plot_4fgl:
-            # Overlay 4FGL sources
+        if plot_fermi_sources:
+            # Overlay Fermi sources
             plot_catalog(
                 equatorial_map, cmap, lower_ra, upper_ra, lower_dec, upper_dec, catalog_path
             )
diff --git a/skyreader/plot/plotting_tools.py b/skyreader/plot/plotting_tools.py
index bbd5c41..28aedf6 100644
--- a/skyreader/plot/plotting_tools.py
+++ b/skyreader/plot/plotting_tools.py
@@ -128,9 +128,9 @@ def hp_ticklabels(zoom=False, lonra=None, latra=None, rot=None, bounds=None):
 
 
 def plot_catalog(master_map, cmap, lower_ra, upper_ra, lower_dec, upper_dec, catalog_path, cmap_min=0., cmap_max=250.):
-    """"Plots the 4FGL catalog in a color that contrasts with the background
+    """"Plots the Fermi catalog in a color that contrasts with the background
     healpix map."""
-    hdu = pyfits.open(catalog_path)  # LAT 14-year from skyreader.constants or user-specified
+    hdu = pyfits.open(catalog_path)  # Fermi catalog defined in skyreader.constants 
     fgl = hdu[1]
     pe = [path_effects.Stroke(linewidth=0.5, foreground=cmap(0.0)),
         path_effects.Normal()]
diff --git a/skyreader/utils/handle_map_data.py b/skyreader/utils/handle_map_data.py
index f6cee0a..7815293 100644
--- a/skyreader/utils/handle_map_data.py
+++ b/skyreader/utils/handle_map_data.py
@@ -331,9 +331,9 @@ def find_filled_pixels(uniqs: np.ndarray):
             nside, next_nside, uniqs
         )
         already_filled_uniqs.append(already_filled_uniqs_nside)
-    already_filled_uniqs = np.concatenate(already_filled_uniqs)
+    already_filled_uniqs_array = np.concatenate(already_filled_uniqs)
     already_filled_indeces = np.array(
-        [np.where(uniqs == uni)[0][0] for uni in already_filled_uniqs]
+        [np.where(uniqs == uni)[0][0] for uni in already_filled_uniqs_array]
     )
     return already_filled_indeces