diff --git a/sotrplib/config/blind_search.py b/sotrplib/config/blind_search.py
index 8d114370..7693359e 100644
--- a/sotrplib/config/blind_search.py
+++ b/sotrplib/config/blind_search.py
@@ -1,6 +1,8 @@
 from abc import ABC, abstractmethod
 from typing import Literal
 
+import astropy.units as u
+from astropydantic import AstroPydanticQuantity
 from numpydantic import NDArray
 from pydantic import BaseModel
 from structlog.types import FilteringBoundLogger
@@ -36,6 +38,7 @@ class PhotutilsBlindSearchConfig(BlindSearchConfig):
     search_type: Literal["photutils"] = "photutils"
     parameters: BlindSearchParameters | None = None
     pixel_mask: NDArray | None = None
+    thumbnail_half_width: AstroPydanticQuantity[u.deg] | None = None
 
     def to_search_provider(
         self, log: FilteringBoundLogger | None = None
@@ -43,6 +46,7 @@ def to_search_provider(
         return SigmaClipBlindSearch(
             parameters=self.parameters,
             pixel_mask=self.pixel_mask,
+            thumbnail_half_width=self.thumbnail_half_width,
             log=log,
         )
 
diff --git a/sotrplib/filters/filters.py b/sotrplib/filters/filters.py
index a89df61b..a3b1e972 100644
--- a/sotrplib/filters/filters.py
+++ b/sotrplib/filters/filters.py
@@ -350,7 +350,6 @@ def matched_filter_depth1_map(
     ):
         # Get band-local versions of the map etc.
         bmap, bivar, bhit = [a[..., r.i1 : r.i2, :] for a in [imap, ivarmap, hit]]
-
         bny, _ = bmap.shape[-2:]
         if shift > 0:
             bS = ShiftMatrix(bmap.shape, bmap.wcs, info.profile)
diff --git a/sotrplib/maps/maps.py b/sotrplib/maps/maps.py
index ca755105..762ebce6 100644
--- a/sotrplib/maps/maps.py
+++ b/sotrplib/maps/maps.py
@@ -83,7 +83,11 @@ def subtract_sources(
         src_model = make_model_source_map(
             input_map.flux,
             sources,
-            nominal_fwhm=get_fwhm(input_map.frequency),
+            nominal_fwhm=get_fwhm(
+                input_map.frequency,
+                arr=input_map.array,
+                instrument=input_map.instrument,
+            ),
             matched_filtered=True,
             verbose=verbose,
             cuts=cuts,
diff --git a/sotrplib/maps/preprocessor.py b/sotrplib/maps/preprocessor.py
index 3d550515..a5e351f2 100644
--- a/sotrplib/maps/preprocessor.py
+++ b/sotrplib/maps/preprocessor.py
@@ -131,11 +131,13 @@ def preprocess(self, input_map: ProcessableMap) -> ProcessableMap:
         rho, kappa = matched_filter_depth1_map(
             imap=input_map.intensity * input_map.intensity_units.to(u.K),
             ivarmap=input_map.inverse_variance / input_map.intensity_units.to(u.K) ** 2,
-            band_center=get_frequency(input_map.frequency),
+            band_center=get_frequency(
+                input_map.frequency, instrument=input_map.instrument
+            ),
             infofile=self.infofile,
             maskfile=self.maskfile,
             source_mask=input_map.mask,
-            beam_fwhm=get_fwhm(input_map.frequency),
+            beam_fwhm=get_fwhm(input_map.frequency, instrument=input_map.instrument),
             beam1d=self.beam1d,
             shrink_holes=self.shrink_holes,
             apod_edge=self.apod_edge,
diff --git a/sotrplib/sifter/core.py b/sotrplib/sifter/core.py
index 4e370221..496059ef 100644
--- a/sotrplib/sifter/core.py
+++ b/sotrplib/sifter/core.py
@@ -153,8 +153,7 @@ def __init__(
         self.ra_jitter = ra_jitter
         self.dec_jitter = dec_jitter
         self.cuts = cuts or {
-            "fwhm_ra": [0.5 * u.arcmin, 5.0 * u.arcmin],
-            "fwhm_dec": [0.5 * u.arcmin, 5.0 * u.arcmin],
+            "fwhm": [0.2, 5.0],
             "snr": [5.0, np.inf],
         }
         self.crossmatch_with_gaia = crossmatch_with_gaia
diff --git a/sotrplib/sifter/crossmatch.py b/sotrplib/sifter/crossmatch.py
index a478860f..243e7867 100644
--- a/sotrplib/sifter/crossmatch.py
+++ b/sotrplib/sifter/crossmatch.py
@@ -1,16 +1,8 @@
-import math
-import os
-import os.path as op
-
 import numpy as np
-import pandas as pd
 from astropy import units as u
-from astropy.coordinates import SkyCoord, angular_separation
-from astropy.table import Table
+from astropy.coordinates import angular_separation
 from astropydantic import AstroPydanticQuantity
-from pixell import enmap
 from pixell import utils as pixell_utils
-from scipy import spatial
 from structlog import get_logger
 from structlog.types import FilteringBoundLogger
 
@@ -187,8 +179,7 @@ def sift(
     map_freq: str | None = None,
     arr: str | None = None,
     cuts: dict = {
-        "fwhm_ra": [0.2 * u.arcmin, 5.0 * u.arcmin],
-        "fwhm_dec": [0.2 * u.arcmin, 5.0 * u.arcmin],
+        "fwhm": [0.2, 5.0],
         "snr": [5.0, np.inf],
     },
     crossmatch_with_gaia: bool = True,
@@ -217,8 +208,8 @@ def sift(
            minimum matching radius, i.e. for a zero flux source, arcmin
        source_fluxes:list = None,
            a list of the fluxes of the extracted sources, if None will pull it from `extracted_sources` dict. 
-       fwhm_cut = 5.0
-           a simple cut on fwhm, in arcmin, above which something is considered noise.
+       cuts: dict
+           snr cut and a simple cut on fwhm, in units of fwhm, outside of range is considered noise.
        ra_jitter:float=0.0
               jitter in the ra direction, in arcmin, to add to the uncertainty of the source position.
        dec_jitter:float=0.0
@@ -239,6 +230,15 @@ def sift(
     """
     log = log if log else get_logger()
     log = log.bind(func_name="sift")
+
+    fwhm = get_fwhm(
+        map_freq, arr=arr, instrument=input_map.instrument if input_map else None
+    )
+    ## apply cuts in units of fwhm
+    if "fwhm" in cuts:
+        cuts["fwhm_ra"] = [cuts["fwhm"][0] * fwhm, cuts["fwhm"][1] * fwhm]
+        cuts["fwhm_dec"] = [cuts["fwhm"][0] * fwhm, cuts["fwhm"][1] * fwhm]
+
     ## log the parameters
     log.info(
         "sift.start",
@@ -251,8 +251,6 @@ def sift(
         crossmatch_with_million_quasar=crossmatch_with_million_quasar,
     )
 
-    fwhm = get_fwhm(map_freq, arr=arr)
-
     angle_coords = u.deg
 
     extracted_ra = np.asarray(
@@ -412,7 +410,7 @@ def sift(
     transient_candidates, new_noise_candidates = recalculate_local_snr(
         transient_candidates,
         input_map,
-        thumb_size=0.25 * u.deg,
+        thumb_size=np.maximum(0.25 * u.deg, 5 * fwhm),
         fwhm=fwhm,
         snr_cut=cuts["snr"][0],
     )
@@ -430,7 +428,12 @@ def sift(
     return source_candidates, transient_candidates, noise_candidates
 
 
-def get_cut_decision(candidate: MeasuredSource, cuts: dict = {}, debug=False) -> bool:
+def get_cut_decision(
+    candidate: MeasuredSource,
+    cuts: dict = {},
+    debug: bool = False,
+    log: FilteringBoundLogger | None = None,
+) -> bool:
     """
     cuts contains dictionary with key values describing the cuts.
 
@@ -439,7 +442,7 @@ def get_cut_decision(candidate: MeasuredSource, cuts: dict = {}, debug=False) ->
     each cut key has [min,max] value on which to make cuts.
 
     for example,
-    if i want a source iwth measured fwhm >0.5 arcmin and less than 5 arcmin
+    if i want a source iwth measured fwhm >0.5  and less than 5 times nominal fwhm,
     cuts={'fwhm':[0.5,5.0]}
     cut = (candidate['fwhm']<cuts['fwhm'][0]) | (candidate['fwhm']>cuts['fwhm'][1])
 
@@ -448,14 +451,22 @@ def get_cut_decision(candidate: MeasuredSource, cuts: dict = {}, debug=False) ->
     True : cut source
     False : within ranges
     """
+    log = log if log else get_logger()
+    log = log.bind(func_name="get_cut_decision")
     cut = False
     for c in cuts.keys():
-        val = getattr(candidate, c)
+        val = getattr(candidate, c, None)
+        if val is None:
+            log.debug("get_cut_decision.missing_value", cut_name=c)
+            continue
         cut |= (val < cuts[c][0]) | (val > cuts[c][1])
         if debug and (val < cuts[c][0]) | (val > cuts[c][1]):
-            print(
-                "cut %s: %.2f < %.2f or %.2f > %.2f"
-                % (c, val, cuts[c][0], val, cuts[c][1])
+            log.debug(
+                "get_cut_decision.cut_made",
+                measured_value=val,
+                cut_name=c,
+                cut_min=cuts[c][0],
+                cut_max=cuts[c][1],
             )
     return cut
 
@@ -675,633 +686,3 @@ def crossmatch_position_and_flux(
         raise ValueError("Some injected sources have mismatched fluxes.")
 
     return matched_mask, similar_fluxes
-
-
-####################
-## Act crossmatching between wafers / bands etc.
-####################
-
-
-def get_cats(
-    path, ctime, return_fnames=False, mfcut=True, renromcut=False, sourcecut=True
-):
-    keys = [
-        "pa4_f220",
-        "pa4_f150",
-        "pa5_f150",
-        "pa5_f090",
-        "pa6_f150",
-        "pa6_f090",
-        "pa7_f040",
-        "pa7_f030",
-    ]
-    cats = []
-    fnames = []
-    for i, key in enumerate(keys):
-        f = op.join(path, f"depth1_{ctime}_{key}_ocat.fits")
-
-        # check if file exists. If not, append empty array
-        if not op.exists(f):
-            cats.append(
-                pd.DataFrame({"ra": [], "dec": [], "flux": [], "snr": [], "ctime": []})
-            )
-            fnames.append(np.nan)
-            continue
-
-        fname = op.join(path, f)
-        t = Table.read(fname)
-        t = t.to_pandas()
-        if len(t) != 0:
-            if mfcut:
-                t = t[t["mf"].astype(bool)]
-            if renromcut:
-                t = t[t["renorm"].astype(bool)]
-            if sourcecut:
-                t = t[~t["source"].astype(bool)]
-        # reindex
-        t.index = np.arange(len(t))
-        cats.append(t)
-        fnames.append(fname)
-
-    cats = dict(zip(keys, cats))
-
-    if return_fnames:
-        return cats, fnames
-
-    else:
-        return cats
-
-
-def crossmatch_array(cats_in_order, N, radius):
-    """crossmatch btw 6 catalogs to get a matched catalog with inds from individual cats, the matched candidate should be detected by at least N arrays
-
-    Args:
-        cats_in_order: a list of 2D np.array of source, each with [[dec,ra]] in deg, order is pa4_f220,pa4_f150, pa5_f150,pa5f090,pa6f150,pa6f090. For non-existin cats, put an empty 2D array in the position in list
-        N: integer that candidate should be detected by at least N arrays
-        radius: radius to search for matches in arcmin
-    """
-    arr_names = ["pa4", "pa4", "pa5", "pa5", "pa6", "pa6", "pa7", "pa7"]
-    # keys = [
-    #     "pa4_f220",
-    #     "pa4_f150",
-    #     "pa5_f150",
-    #     "pa5_f090",
-    #     "pa6_f150",
-    #     "pa6_f090",
-    #     "pa7_f040",
-    #     "pa7_f030",
-    # ]
-    rad_deg = radius / 60.0
-    skip = 0
-    cat_full = []
-    for i, cat in enumerate(cats_in_order):
-        if cat.shape[0] == 0:
-            skip += 1
-            continue
-        else:
-            patch1 = np.arange(0, cat.shape[0], dtype=int)[np.newaxis].T
-            patch2 = np.ones(shape=(cat.shape[0], 1)) * i  # i for index of the cats
-            cat = np.concatenate((cat, patch1), axis=1)
-            cat = np.concatenate((cat, patch2), axis=1)
-            if i - skip == 0:
-                cat_full = cat
-            else:
-                cat_full = np.concatenate((cat_full, cat), axis=0)
-    del cat
-    if len(cat_full) == 0:
-        new_cat = np.zeros((0, 8))
-    else:
-        tree = spatial.cKDTree(cat_full[:, 0:2])
-        idups_raw = tree.query_ball_tree(tree, rad_deg)
-        new_cat = []
-        idups_merge = []
-        for i, srcs in enumerate(idups_raw):
-            if i == 0:
-                idups_merge.append(srcs)
-            else:
-                overlap = 0
-                for srcs2 in idups_merge:
-                    if bool(set(srcs) & set(srcs2)):
-                        srcs_m = list(set(srcs + srcs2))
-                        idups_merge.remove(srcs2)
-                        idups_merge.append(srcs_m)
-                        overlap = 1
-                if overlap == 0:
-                    idups_merge.append(srcs)
-        idups_redu = []
-        for i, srcs in enumerate(idups_merge):
-            inds = cat_full[srcs, -1].astype(int)
-            arrs = np.unique(np.array(arr_names)[inds])
-            if arrs.shape[0] > N - 1 and srcs not in idups_redu:
-                row = np.ones(8) * float("nan")
-                for src in srcs:
-                    ind = int(cat_full[src, -1])
-                    row[ind] = cat_full[src, -2]
-                new_cat.append(row)
-                idups_redu.append(srcs)
-
-    new_cat = np.array(new_cat)
-    if len(new_cat) == 0:
-        new_cat = np.zeros((0, 8))
-    ocat_dtype = [
-        ("pa4_f220", "1f"),
-        ("pa4_f150", "1f"),
-        ("pa5_f150", "1f"),
-        ("pa5_f090", "1f"),
-        ("pa6_f150", "1f"),
-        ("pa6_f090", "1f"),
-        ("pa7_f040", "1f"),
-        ("pa7_f030", "1f"),
-    ]
-    ocat = np.zeros(len(new_cat), ocat_dtype).view(np.recarray)
-    ocat.pa4_f220 = new_cat[:, 0]
-    ocat.pa4_f150 = new_cat[:, 1]
-    ocat.pa5_f150 = new_cat[:, 2]
-    ocat.pa5_f090 = new_cat[:, 3]
-    ocat.pa6_f150 = new_cat[:, 4]
-    ocat.pa6_f090 = new_cat[:, 5]
-    ocat.pa7_f040 = new_cat[:, 6]
-    ocat.pa7_f030 = new_cat[:, 7]
-    return ocat
-
-
-def crossmatch_array_new(cats_in_order, N, radius, ctime):
-    """crossmatch btw 8 catalogs to get a matched catalog with inds from individual cats, the matched candidate should be detected by at least N arrays
-       but it would still keep the catalog crossmatched by two bands if there's only 1 array's data available
-
-    Args:
-        cats_in_order: a list of 2D np.array of source, each with [[dec,ra]] in deg, order is pa4_f220,pa4_f150, pa5_f150,pa5f090,pa6f150,pa6f090. For non-existin cats, put an empty 2D array in the position in list
-        N: integer that candidate should be detected by at least N arrays
-        radius: radius to search for matches in arcmin
-    """
-    from ..utils.utils import obj_dist
-
-    arr_names = ["pa4", "pa4", "pa5", "pa5", "pa6", "pa6", "pa7", "pa7"]
-    keys = [
-        "pa4_f220",
-        "pa4_f150",
-        "pa5_f150",
-        "pa5_f090",
-        "pa6_f150",
-        "pa6_f090",
-        "pa7_f040",
-        "pa7_f030",
-    ]
-    rad_deg = radius / 60.0
-    skip = 0
-    cat_full = []
-    for i, cat in enumerate(cats_in_order):
-        if cat.shape[0] == 0:
-            skip += 1
-            continue
-        else:
-            patch1 = np.arange(0, cat.shape[0], dtype=int)[np.newaxis].T
-            patch2 = np.ones(shape=(cat.shape[0], 1)) * i  # i for index of the cats
-            cat = np.concatenate((cat, patch1), axis=1)
-            cat = np.concatenate((cat, patch2), axis=1)
-            if i - skip == 0:
-                cat_full = cat
-            else:
-                cat_full = np.concatenate((cat_full, cat), axis=0)
-    del cat
-    if len(cat_full) == 0:
-        new_cat = np.zeros((0, 8))
-    else:
-        tree = spatial.cKDTree(cat_full[:, 0:2])
-        idups_raw = tree.query_ball_tree(tree, rad_deg)
-        new_cat = []
-        idups_merge = []
-        for i, srcs in enumerate(idups_raw):
-            if i == 0:
-                idups_merge.append(srcs)
-            else:
-                overlap = 0
-                for srcs2 in idups_merge:
-                    if bool(set(srcs) & set(srcs2)):
-                        srcs_m = list(set(srcs + srcs2))
-                        idups_merge.remove(srcs2)
-                        idups_merge.append(srcs_m)
-                        overlap = 1
-                if overlap == 0:
-                    idups_merge.append(srcs)
-        idups_redu = []
-        for i, srcs in enumerate(idups_merge):
-            inds = cat_full[srcs, -1].astype(int)
-            arrs = np.unique(np.array(arr_names)[inds])
-            keys_srcs = np.unique(np.array(keys)[inds])
-            if arrs.shape[0] > N - 1 and srcs not in idups_redu:
-                row = np.ones(8) * float("nan")
-                for src in srcs:
-                    ind = int(cat_full[src, -1])
-                    row[ind] = cat_full[src, -2]
-                new_cat.append(row)
-                idups_redu.append(srcs)
-            # check if all other freq arr data is not available
-            elif arrs.shape[0] == 1 and keys_srcs.shape[0] == 2:
-                src = srcs[0]
-                ra_deg = cat_full[src, 0]
-                dec_deg = cat_full[src, 1]
-                keys_left = [key for key in keys if key not in keys_srcs]
-                count_ava = 0
-                for key in keys_left:
-                    ivar_file = (
-                        "/scratch/gpfs/snaess/actpol/maps/depth1/release/%s/depth1_%s_%s_ivar.fits"
-                        % (ctime[:5], ctime, key)
-                    )
-                    if os.path.exists(ivar_file):
-                        ivar = enmap.read_map(ivar_file)
-                        if enmap.contains(
-                            ivar.shape,
-                            ivar.wcs,
-                            [
-                                dec_deg * pixell_utils.degree,
-                                ra_deg * pixell_utils.degree,
-                            ],
-                        ):
-                            dist = obj_dist(ivar, np.array([[dec_deg, ra_deg]]))[0]
-                            dec_pix, ra_pix = ivar.sky2pix(
-                                [
-                                    dec_deg * pixell_utils.degree,
-                                    ra_deg * pixell_utils.degree,
-                                ]
-                            )
-                            ivar_src = ivar[int(dec_pix), int(ra_pix)]
-                            if ivar_src != 0.0 and dist > 10.0:
-                                count_ava += 1
-                if count_ava == 0 and srcs not in idups_redu:
-                    row = np.ones(8) * float("nan")
-                    for src in srcs:
-                        ind = int(cat_full[src, -1])
-                        row[ind] = cat_full[src, -2]
-                    new_cat.append(row)
-                    idups_redu.append(srcs)
-    new_cat = np.array(new_cat)
-    if len(new_cat) == 0:
-        new_cat = np.zeros((0, 8))
-    ocat_dtype = [
-        ("pa4_f220", "1f"),
-        ("pa4_f150", "1f"),
-        ("pa5_f150", "1f"),
-        ("pa5_f090", "1f"),
-        ("pa6_f150", "1f"),
-        ("pa6_f090", "1f"),
-        ("pa7_f040", "1f"),
-        ("pa7_f030", "1f"),
-    ]
-    ocat = np.zeros(len(new_cat), ocat_dtype).view(np.recarray)
-    ocat.pa4_f220 = new_cat[:, 0]
-    ocat.pa4_f150 = new_cat[:, 1]
-    ocat.pa5_f150 = new_cat[:, 2]
-    ocat.pa5_f090 = new_cat[:, 3]
-    ocat.pa6_f150 = new_cat[:, 4]
-    ocat.pa6_f090 = new_cat[:, 5]
-    ocat.pa7_f040 = new_cat[:, 6]
-    ocat.pa7_f030 = new_cat[:, 7]
-    return ocat
-
-
-def crossmatch_freq_array(cats_in_order, N, radius):
-    """crossmatch btw 6 catalogs to get a matched catalog with inds from individual cats, the matched candidate should be detected by at least N arrays
-
-    Args:
-        cats_in_order: a list of 2D np.array of source, each with [[dec,ra]] in deg, order is pa4_f220,pa4_f150, pa5_f150,pa5f090,pa6f150,pa6f090. For non-existin cats, put an empty 2D array in the position in list
-        N: integer that candidate should be detected by at least N arrays
-        radius: radius to search for matches in arcmin
-    """
-    # arr_names = ["pa4", "pa4", "pa5", "pa5", "pa6", "pa6", "pa7", "pa7"]
-    keys = [
-        "pa4_f220",
-        "pa4_f150",
-        "pa5_f150",
-        "pa5_f090",
-        "pa6_f150",
-        "pa6_f090",
-        "pa7_f040",
-        "pa7_f030",
-    ]
-    rad_deg = radius / 60.0
-    skip = 0
-    cat_full = []
-    for i, cat in enumerate(cats_in_order):
-        if cat.shape[0] == 0:
-            skip += 1
-            continue
-        else:
-            patch1 = np.arange(0, cat.shape[0], dtype=int)[np.newaxis].T
-            patch2 = np.ones(shape=(cat.shape[0], 1)) * i  # i for index of the cats
-            cat = np.concatenate((cat, patch1), axis=1)
-            cat = np.concatenate((cat, patch2), axis=1)
-            if i - skip == 0:
-                cat_full = cat
-            else:
-                cat_full = np.concatenate((cat_full, cat), axis=0)
-    del cat
-    if len(cat_full) == 0:
-        new_cat = np.zeros((0, 8))
-    else:
-        tree = spatial.cKDTree(cat_full[:, 0:2])
-        idups_raw = tree.query_ball_tree(tree, rad_deg)
-        new_cat = []
-        idups_merge = []
-        for i, srcs in enumerate(idups_raw):
-            if i == 0:
-                idups_merge.append(srcs)
-            else:
-                overlap = 0
-                for srcs2 in idups_merge:
-                    if bool(set(srcs) & set(srcs2)):
-                        srcs_m = list(set(srcs + srcs2))
-                        idups_merge.remove(srcs2)
-                        idups_merge.append(srcs_m)
-                        overlap = 1
-                if overlap == 0:
-                    idups_merge.append(srcs)
-        idups_redu = []
-        for i, srcs in enumerate(idups_merge):
-            inds = cat_full[srcs, -1].astype(int)
-            keys_srcs = np.unique(np.array(keys)[inds])
-            if keys_srcs.shape[0] > N - 1 and srcs not in idups_redu:
-                row = np.ones(8) * float("nan")
-                for src in srcs:
-                    ind = int(cat_full[src, -1])
-                    row[ind] = cat_full[src, -2]
-                new_cat.append(row)
-                idups_redu.append(srcs)
-
-    new_cat = np.array(new_cat)
-    if len(new_cat) == 0:
-        new_cat = np.zeros((0, 8))
-    ocat_dtype = [
-        ("pa4_f220", "1f"),
-        ("pa4_f150", "1f"),
-        ("pa5_f150", "1f"),
-        ("pa5_f090", "1f"),
-        ("pa6_f150", "1f"),
-        ("pa6_f090", "1f"),
-        ("pa7_f040", "1f"),
-        ("pa7_f030", "1f"),
-    ]
-    ocat = np.zeros(len(new_cat), ocat_dtype).view(np.recarray)
-    ocat.pa4_f220 = new_cat[:, 0]
-    ocat.pa4_f150 = new_cat[:, 1]
-    ocat.pa5_f150 = new_cat[:, 2]
-    ocat.pa5_f090 = new_cat[:, 3]
-    ocat.pa6_f150 = new_cat[:, 4]
-    ocat.pa6_f090 = new_cat[:, 5]
-    ocat.pa7_f040 = new_cat[:, 6]
-    ocat.pa7_f030 = new_cat[:, 7]
-    return ocat
-
-
-def crossmatch_col(cats, cross_idx):
-    new_cats = []
-    keys = [
-        "pa4_f220",
-        "pa4_f150",
-        "pa5_f150",
-        "pa5_f090",
-        "pa6_f150",
-        "pa6_f090",
-        "pa7_f040",
-        "pa7_f030",
-    ]
-    for i, df in enumerate(cats):
-        indices = cross_idx[keys[i]]
-        indices = [x for x in indices if str(x) != "nan"]
-        df["crossmatch"] = False
-        df.loc[indices, "crossmatch"] = True
-        new_cats.append(df)
-    return new_cats
-
-
-def merge_cats(full_cats_in_order, N, radius, saveps=False, type="arr", ctime=None):
-    """merge btw 6 catalogs to get one single merged catalog, with on overall ra, dec, position error
-
-    Args:
-        full_cats_in_order: fits.open(file)[1].data from fits file, order is pa4_f220,pa4_f150, pa5_f150,pa5f090,pa6f150,pa6f090.
-        N: integer that candidate should be detected by at least N arrays
-        radius: radius to search for matches in arcmin
-        saveps: if True, output point sources
-        type: arr if crossmatch by array, freq if crossmatch by freq or arr, cond if crossmatch by array but use freq if only one array's data is available
-        ctime: ctime in the file name with string format, needed if type is cond
-    """
-    cats_pos = []
-    for i in range(8):
-        cat_data = full_cats_in_order[i]
-        cat_array = np.zeros((cat_data.shape[0], 2))
-        cat_array[:, 0] = cat_data.ra
-        cat_array[:, 1] = cat_data.dec
-        cats_pos.append(cat_array)
-    if type == "arr":
-        matched_inds = crossmatch_array(cats_pos, N, radius)
-    elif type == "freq":
-        matched_inds = crossmatch_freq_array(cats_pos, N, radius)
-    elif type == "cond":
-        matched_inds = crossmatch_array_new(cats_pos, N, radius, ctime)
-    else:
-        raise ValueError("type must be arr, freq, or cond")
-    new_cats = crossmatch_col(full_cats_in_order, matched_inds)
-    cat_merge = []
-    for i in range(matched_inds.shape[0]):
-        ra = 0
-        dec = 0
-        invar = 0
-        # count = 0
-        ps = 0
-        for j, [arr, freq] in enumerate(
-            [
-                ["pa4", "f220"],
-                ["pa4", "f150"],
-                ["pa5", "f150"],
-                ["pa5", "f090"],
-                ["pa6", "f150"],
-                ["pa6", "f090"],
-                ["pa7", "f040"],
-                ["pa7", "f030"],
-            ]
-        ):
-            key = "%s_%s" % (arr, freq)
-            if not math.isnan(matched_inds[key][i]):
-                inputs = None  # TODO: Fix
-                fwhm = inputs.get_fwhm_arcmin(arr, freq)
-                fwhm_deg = fwhm / 60.0
-                idx_src = int(matched_inds[key][i])
-                cat_src = full_cats_in_order[j]
-                ps += cat_src.source[idx_src]
-                ra_src = cat_src.ra[idx_src]
-                dec_src = cat_src.dec[idx_src]
-                snr_src = cat_src.snr[idx_src]
-                pos_err = fwhm_deg / snr_src
-                ra += ra_src / pos_err**2.0
-                dec += dec_src / pos_err**2.0
-                invar += 1 / pos_err**2.0
-        if invar == 0:
-            print("ZeroDivisionError: Candidate Removed")
-            continue
-        ra /= invar
-        dec /= invar
-        var_w = 0
-        if (ps > 0) and not saveps:
-            continue
-        for j, [arr, freq] in enumerate(
-            [
-                ["pa4", "f220"],
-                ["pa4", "f150"],
-                ["pa5", "f150"],
-                ["pa5", "f090"],
-                ["pa6", "f150"],
-                ["pa6", "f090"],
-                ["pa7", "f040"],
-                ["pa7", "f030"],
-            ]
-        ):
-            key = "%s_%s" % (arr, freq)
-            if not math.isnan(matched_inds[key][i]):
-                idx_src = int(matched_inds[key][i])
-                ra_src = full_cats_in_order[j].ra[idx_src]
-                dec_src = full_cats_in_order[j].dec[idx_src]
-                dist = pixell_utils.angdist(
-                    np.array([ra_src, dec_src]) * pixell_utils.degree,
-                    np.array([ra, dec]) * pixell_utils.degree,
-                )
-                inputs = None  # TODO: FIX
-                fwhm = inputs.get_fwhm_arcmin(arr, freq)
-                fwhm_deg = fwhm / 60.0
-                pos_err = fwhm_deg / snr_src
-                var_w += (dist**2.0) * (1 / pos_err**2.0)
-        var_w /= invar
-        pos_err_arcmin = var_w**0.5 * 60 / pixell_utils.degree
-        cat_merge.append(
-            np.array(
-                [
-                    ra,
-                    dec,
-                    pos_err_arcmin,
-                    matched_inds["pa4_f220"][i],
-                    matched_inds["pa4_f150"][i],
-                    matched_inds["pa5_f150"][i],
-                    matched_inds["pa5_f090"][i],
-                    matched_inds["pa6_f150"][i],
-                    matched_inds["pa6_f090"][i],
-                    matched_inds["pa7_f040"][i],
-                    matched_inds["pa7_f030"][i],
-                ]
-            )
-        )
-    if len(cat_merge) == 0:
-        cat_merge = np.zeros((0, 11))
-    else:
-        cat_merge = np.vstack(cat_merge)
-    ocat_dtype = [
-        ("ra", "1f"),
-        ("dec", "1f"),
-        ("pos_err_arcmin", "1f"),
-        ("pa4_f220", "1f"),
-        ("pa4_f150", "1f"),
-        ("pa5_f150", "1f"),
-        ("pa5_f090", "1f"),
-        ("pa6_f150", "1f"),
-        ("pa6_f090", "1f"),
-        ("pa7_f040", "1f"),
-        ("pa7_f030", "1f"),
-    ]
-    ocat = np.zeros(len(cat_merge), ocat_dtype).view(np.recarray)
-    ocat.ra = cat_merge[:, 0]
-    ocat.dec = cat_merge[:, 1]
-    ocat.pos_err_arcmin = cat_merge[:, 2]
-    ocat.pa4_f220 = cat_merge[:, 3]
-    ocat.pa4_f150 = cat_merge[:, 4]
-    ocat.pa5_f150 = cat_merge[:, 5]
-    ocat.pa5_f090 = cat_merge[:, 6]
-    ocat.pa6_f150 = cat_merge[:, 7]
-    ocat.pa6_f090 = cat_merge[:, 8]
-    ocat.pa7_f040 = cat_merge[:, 9]
-    ocat.pa7_f030 = cat_merge[:, 10]
-    return ocat, new_cats
-
-
-def candidate_output(
-    cross_ocat: pd.DataFrame,
-    cats: dict,
-    odir: str,
-    ctime: float,
-    verbosity=0,
-    overwrite=True,
-):
-    # for each candidate, write output
-    keys = list(cats.keys())
-    for i in range(len(cross_ocat)):
-        ocat = {
-            "idx": [],
-            "band": [],
-            "ra": [],
-            "dec": [],
-            "flux": [],
-            "dflux": [],
-            "snr": [],
-            "ctime": [],
-        }
-
-        # create output dir
-        can_odir = op.join(odir, f"{ctime}_{str(i)}")
-        if not op.exists(can_odir):
-            os.makedirs(can_odir)
-
-        # get ra, dec, flux, dflux, snr, ctime for each candidate
-        for k in keys:
-            # check if detection exists
-            if np.isnan(cross_ocat[k][i]):
-                continue
-
-            # get ra, dec, flux, dflux, snr, ctime
-            ocat["idx"].append(cross_ocat[k][i])
-            ocat["band"].append(k)
-            ocat["ra"].append(cats[k]["ra"][int(cross_ocat[k][i])])
-            ocat["dec"].append(cats[k]["dec"][int(cross_ocat[k][i])])
-            ocat["flux"].append(cats[k]["flux"][int(cross_ocat[k][i])])
-            ocat["dflux"].append(cats[k]["dflux"][int(cross_ocat[k][i])])
-            ocat["snr"].append(cats[k]["snr"][int(cross_ocat[k][i])])
-            ocat["ctime"].append(cats[k]["ctime"][int(cross_ocat[k][i])])
-
-        # Convert to table
-        output = pd.DataFrame(ocat)
-        output = Table.from_pandas(output)
-
-        # metadata
-        output.meta["ra_deg"] = cross_ocat["ra"][i]
-        output.meta["dec_deg"] = cross_ocat["dec"][i]
-        output.meta["err_am"] = cross_ocat["pos_err_arcmin"][i]
-        output.meta["maptime"] = ctime
-
-        # write output
-        output.write(
-            op.join(can_odir, f"{ctime}_{str(i)}_ocat.fits"), overwrite=overwrite
-        )
-
-        if verbosity > 0:
-            print(f"Wrote {op.join(can_odir, f'{ctime}_{str(i)}_ocat.fits')}")
-
-    return
-
-
-def blazar_crossmatch(pos, tol=0.5):
-    """
-    Mask blazar 3C 454.3 using icrs position from SIMBAD
-
-    Args:
-        pos:  np.array of positions [[dec, ra]] in deg
-        tol: tolerance in deg
-
-    Returns:
-        1 in blazar, 0 if not in blazar
-    """
-    blazar_pos = SkyCoord("22h53m57.7480438728s", "+16d08m53.561508864s", frame="icrs")
-    can_pos = SkyCoord(pos[:, 1], pos[:, 0], frame="icrs", unit="deg")
-    sep = blazar_pos.separation(can_pos).to(u.deg)
-    match = []
-    for s in sep:
-        if s.value < tol:
-            match.append(1)
-        else:
-            match.append(0)
-    return np.array(match).astype(bool)
diff --git a/sotrplib/sims/sim_maps.py b/sotrplib/sims/sim_maps.py
index 682a766b..e53af1ef 100644
--- a/sotrplib/sims/sim_maps.py
+++ b/sotrplib/sims/sim_maps.py
@@ -283,6 +283,7 @@ def inject_sources(
     observation_time: float | enmap.ndmap,
     freq: str = "f090",
     arr: str = None,
+    instrument: str = None,
     debug: bool = False,
     log: FilteringBoundLogger | None = None,
 ):
@@ -327,7 +328,7 @@ def inject_sources(
         raise ValueError("Input map must be a ProcessableMap or enmap.ndmap object")
 
     mapres = abs(wcs.wcs.cdelt[0]) * u.deg
-    fwhm = get_fwhm(freq, arr=arr)
+    fwhm = get_fwhm(freq, arr=arr, instrument=instrument)
     fwhm_pixels = (fwhm / mapres).value
 
     log.info("inject_sources.injecting_sources", n_sources=len(sources))
@@ -496,7 +497,9 @@ def inject_simulated_sources(
     catalog_sources += inject_random_sources(
         mapdata,
         sim_params,
-        fwhm_arcmin=get_fwhm(mapdata.frequency, arr=mapdata.array),
+        fwhm_arcmin=get_fwhm(
+            mapdata.frequency, arr=mapdata.array, instrument=mapdata.instrument
+        ),
         add_noise=False,
         log=log,
     )
diff --git a/sotrplib/sources/blind.py b/sotrplib/sources/blind.py
index 688ed2cc..9abc7f58 100644
--- a/sotrplib/sources/blind.py
+++ b/sotrplib/sources/blind.py
@@ -41,10 +41,14 @@ def __init__(
         log: FilteringBoundLogger | None = None,
         parameters: BlindSearchParameters | None = None,
         pixel_mask: enmap.ndmap | None = None,
+        thumbnail_half_width: AstroPydanticQuantity[u.deg] = AstroPydanticQuantity(
+            u.Quantity(0.1, "deg")
+        ),
     ):
         self.log = log or get_logger()
         self.parameters = parameters if parameters else BlindSearchParameters()
         self.pixel_mask = pixel_mask
+        self.thumbnail_half_width = thumbnail_half_width
         self.log.info("sigma_clip_blind_search.initialized", parameters=self.parameters)
 
     def search(
@@ -75,7 +79,10 @@ def search(
         )
         for source in extracted_sources:
             ## want to extract the thumbnail at the map location, but then apply ra,dec offsets
-            source.extract_thumbnail(input_map, reproject_thumb=True)
+            ## set thumbnail radius to 0.1 deg by default, but for large fwhm do 5*fwhm (i.e. for SAT it will be ~1.5 deg)
+            source.extract_thumbnail(
+                input_map, thumb_width=self.thumbnail_half_width, reproject_thumb=True
+            )
             if pointing_residuals:
                 source_pos = pointing_residuals.apply_offset_at_position(
                     SkyCoord(ra=source.ra, dec=source.dec)
diff --git a/sotrplib/sources/force.py b/sotrplib/sources/force.py
index 2bc9f1c2..66589502 100644
--- a/sotrplib/sources/force.py
+++ b/sotrplib/sources/force.py
@@ -146,7 +146,11 @@ def force(
         catalogs: list[SourceCatalog],
         pointing_residuals: MapPointingOffset | None = None,
     ) -> list[MeasuredSource]:
-        fwhm = get_fwhm(freq=input_map.frequency, arr=input_map.array)
+        fwhm = get_fwhm(
+            freq=input_map.frequency,
+            arr=input_map.array,
+            instrument=input_map.instrument,
+        )
         source_list = list(
             itertools.chain(*[c.forced_photometry_sources(input_map) for c in catalogs])
         )
@@ -253,7 +257,11 @@ def __init__(
     def force(
         self, input_map: ProcessableMap, catalogs: list[SourceCatalog]
     ) -> list[MeasuredSource]:
-        fwhm = get_fwhm(freq=input_map.frequency, arr=input_map.array)
+        fwhm = get_fwhm(
+            freq=input_map.frequency,
+            arr=input_map.array,
+            instrument=input_map.instrument,
+        )
         ## get all forced photometry sources from catalogs
         source_list = list(
             itertools.chain(*[c.forced_photometry_sources(input_map) for c in catalogs])
diff --git a/sotrplib/sources/forced_photometry.py b/sotrplib/sources/forced_photometry.py
index 3f2caedc..60e565b7 100644
--- a/sotrplib/sources/forced_photometry.py
+++ b/sotrplib/sources/forced_photometry.py
@@ -507,7 +507,11 @@ def scipy_2d_gaussian_fit(
             reprojected=reproject_thumb,
             fwhm_guess=fwhm
             if fwhm is not None
-            else get_fwhm(freq=input_map.frequency, arr=input_map.array),
+            else get_fwhm(
+                freq=input_map.frequency,
+                arr=input_map.array,
+                instrument=input_map.instrument,
+            ),
             force_center=forced_source.flux < flux_lim_fit_centroid
             if forced_source.flux is not None
             else False,
diff --git a/sotrplib/utils/utils.py b/sotrplib/utils/utils.py
index bb0aa2cf..f9302727 100644
--- a/sotrplib/utils/utils.py
+++ b/sotrplib/utils/utils.py
@@ -273,9 +273,11 @@ def get_map_groups(
     return map_groups, map_group_time_range, time_bins
 
 
-def get_frequency(freq: str, arr: str = None):
+def get_frequency(freq: str, arr: str = None, instrument: str | None = None):
+    ## instrument can be SOLAT or SOSAT, defaults to SOLAT
     # Array not used yet, but could have per array/freq band centers
     ## from lat white paper https://arxiv.org/pdf/2503.00636
+
     frequency = {
         "f030": 30 * u.GHz,
         "f040": 40 * u.GHz,
@@ -284,11 +286,19 @@ def get_frequency(freq: str, arr: str = None):
         "f220": 220 * u.GHz,
         "f280": 280 * u.GHz,
     }
+    if instrument is not None and instrument.upper() == "SOSAT":
+        frequency = {
+            "f090": 93 * u.GHz,
+            "f150": 145 * u.GHz,
+            "f220": 225 * u.GHz,
+            "f280": 280 * u.GHz,
+        }
 
     return frequency[freq]
 
 
-def get_fwhm(freq: str, arr: str = None):
+def get_fwhm(freq: str, arr: str = None, instrument: str | None = None):
+    ## instrument can be SOLAT or SOSAT , defaults to SOLAT
     # Array not used yet, but could have per array/freq fwhm
     ## from lat white paper https://arxiv.org/pdf/2503.00636
     fwhm = {
@@ -299,6 +309,13 @@ def get_fwhm(freq: str, arr: str = None):
         "f220": 1.0 * u.arcmin,
         "f280": 0.9 * u.arcmin,
     }
+    if instrument is not None and instrument.upper() == "SOSAT":
+        fwhm = {
+            "f090": 27.4 * u.arcmin,
+            "f150": 17.6 * u.arcmin,
+            "f220": 13.5 * u.arcmin,
+            "f280": 12.1 * u.arcmin,
+        }
 
     return fwhm[freq]
 
@@ -353,6 +370,7 @@ def get_cut_radius(
     map_res_arcmin: float,
     arr: str,
     freq: str,
+    instrument: str | None = None,
     fwhm_arcmin: float = None,
     matched_filtered: bool = False,
     source_amplitudes: list = [],
@@ -381,7 +399,7 @@ def get_cut_radius(
     """
 
     if not fwhm_arcmin:
-        fwhm_arcmin = get_fwhm(freq, arr).to_value(u.arcmin)
+        fwhm_arcmin = get_fwhm(freq, arr, instrument=instrument).to_value(u.arcmin)
 
     ## I guess this is needed for filtering wings?
     mf_factor = 2**0.5 if matched_filtered else 1.0