diff --git a/devel/roman/diagnose_bilinear_vs_direct.py b/devel/roman/diagnose_bilinear_vs_direct.py
index 68eaa4e1..e9fbcd75 100644
--- a/devel/roman/diagnose_bilinear_vs_direct.py
+++ b/devel/roman/diagnose_bilinear_vs_direct.py
@@ -96,7 +96,7 @@ def run_scan(model, scas, grid, stamp_size, scale, margin, aberrations):
 
                 prof = galsim.roman.getPSF(
                     int(sca),
-                    model.filter,
+                    model.filter_name,
                     SCA_pos=star.image_pos,
                     pupil_bin=piff.roman_psf.pupil_bin,
                     wcs=star.image.wcs,
@@ -104,12 +104,7 @@ def run_scan(model, scas, grid, stamp_size, scale, margin, aberrations):
                     wavelength=None if model.chromatic else model.bandpass.effective_wavelength,
                 )
                 direct_image = star.image.copy()
-                prof.drawImage(
-                    direct_image,
-                    method=model._method,
-                    center=star.image_pos,
-                    bandpass=model.bandpass if model.chromatic else None,
-                )
+                model._draw_profile_to_image(prof, direct_image, center=star.image_pos, star=star)
 
                 max_abs, rms, frac_peak, frac_l2 = point_metrics(
                     model_star.image.array, direct_image.array
@@ -301,12 +296,12 @@ def main():
 
     piff.roman_psf.pupil_bin = int(args.pupil_bin)
     model = piff.Roman(
-        filter="H158",
         chromatic=False,
         max_zernike=args.max_zernike,
         aberration_interp="constant",
         aberration_prior_sigma=1.0e6,
     )
+    model.set_bandpass(galsim.roman.getBandpasses()["H158"])
 
     scas = parse_scas(args.scas)
     aberrations = parse_aberrations(args.aberrations, model.param_len)
diff --git a/devel/roman/piff.yaml b/devel/roman/piff.yaml
index fb545ed4..39662aa8 100644
--- a/devel/roman/piff.yaml
+++ b/devel/roman/piff.yaml
@@ -1,3 +1,5 @@
+# Note: this run takes ~ 80 minutes to run.
+
 input:
 
     # This file is quite large, so not saved in the repo.
@@ -9,6 +11,10 @@ input:
 
     cat_file_name: 'stars_13906_test_v251105.parquet'
 
+    bandpass:
+        type: RomanBandpass
+        name: H158
+
     ra_col: ra
     dec_col: dec
     properties:
@@ -26,16 +32,15 @@ select:
     hsm_size_reject: 5    # number of inter-quartile ranges from median to reject size estimate.
     max_pixel_cut: 50000  # reject stars with a pixel value higher than this. (sort of)
     min_sep: 1.0          # reject stars with a neighbor closer than 1 arcsec away.
-    #nstars: 100
 
 output:
     dir: 'output'
-    file_name: 'ffov_13906_11.piff'
+    file_name: 'ffov_13906_17.piff'
 
     stats:
         - 
           type: StarImages
-          file_name: 'stars_11.png'
+          file_name: 'stars_17.png'
           nplot: 0  # all stars
 
 psf:
@@ -46,12 +51,11 @@ psf:
 
       -
         type: RomanOptics
-        filter: H158
         chromatic: False
         max_zernike: 12    # Fit Zernike coefficients 4-12 inclusive.
         aberration_prior_sigma: 3.e-3  # Gaussian prior on all Zernike coefficients in absolute units
         aberration_interp: constant    # linear is possible but doesn't seem to work better.
-        nominal_interp: bilinear       # five_points is possible but doesn't work better on
+        nominal_interp: bilinear       # five_point is possible but doesn't work better on
                                        # this sim. May be different on real data.
         nproc: 6
 
diff --git a/devel/roman/piff_chromatic.yaml b/devel/roman/piff_chromatic.yaml
new file mode 100644
index 00000000..17d900a4
--- /dev/null
+++ b/devel/roman/piff_chromatic.yaml
@@ -0,0 +1,88 @@
+# Note: this run takes ~ 380 minutes to run.
+
+input:
+
+    # This file is quite large, so not saved in the repo.
+    # It and the cat_file are available from here:
+    #   https://drive.google.com/drive/folders/1akvHjdKSTMppPTcNfXjEy9SVq_hFBDBc
+    image_file_name: 'ffov_13906_test_v251020.fits'
+    image_hdu: [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]  # hdu = SCA
+    badpix_zeros: True
+
+    cat_file_name: 'stars_13906_test_v251105.parquet'
+
+    # For now just use a single SED for all stars.  This is the SED for a K type star.
+    sed_file_name: 'seds/XSL_DR3_release/xsl_spectrum_X0066_merged_scl.fits'
+    sed_wave_type: 'nm'
+    sed_flux_type: 'erg cm**(-2) s**(-1) angstrom**(-1)'
+    sed_max_samples: 2  # Just keep a linear approximation to the sed across the bandpass.
+    bandpass:
+        type: RomanBandpass
+        name: H158
+
+    ra_col: ra
+    dec_col: dec
+    properties:
+        sca: '$@image_num + 1'
+    noise: BKGNDVAR
+    gain: EQVGAIN
+    sky: median
+
+    stamp_size: 25
+
+select:
+
+    min_snr: 50           # reject very faint stars
+    max_snr_weight: 500   # don't over-weight very bright stars
+    hsm_size_reject: 5    # number of inter-quartile ranges from median to reject size estimate.
+    max_pixel_cut: 50000  # reject stars with a pixel value higher than this. (sort of)
+    min_sep: 1.0          # reject stars with a neighbor closer than 1 arcsec away.
+
+output:
+    dir: 'output'
+    file_name: 'ffov_13906_18.piff'
+
+    stats:
+        - 
+          type: StarImages
+          file_name: 'stars_18.png'
+          nplot: 0  # all stars
+
+psf:
+
+    type: Sum
+
+    components:
+
+      -
+        type: RomanOptics
+        chromatic: True
+        max_zernike: 12    # Fit Zernike coefficients 4-12 inclusive.
+        aberration_prior_sigma: 3.e-3  # Gaussian prior on all Zernike coefficients.
+        aberration_interp: constant    # linear is possible but doesn't seem to work better.
+        nominal_interp: bilinear       # five_point is possible but doesn't work better on
+                                       # this sim. May be different on real data.
+        nproc: 6
+
+      - 
+        type: Simple
+        model:
+            type: PixelGrid
+            scale: 0.11
+            size: 5
+        interp:
+            type: BasisPolynomial
+            order: 2
+
+    outliers:
+      -
+        type: Chisq
+        nsigma: 10
+        max_remove: 0.03
+      -
+        type: Centroid
+        max_offset: 0.2 # arcsec
+
+
+  
+verbose: 2
diff --git a/piff/config.py b/piff/config.py
index ec827d69..5398ced6 100644
--- a/piff/config.py
+++ b/piff/config.py
@@ -159,11 +159,12 @@ def process(config, logger=None):
             config['select'][key] = config['input'].pop(key)
 
     # read in the input images
-    objects, wcs, pointing = Input.process(config['input'], logger=logger)
+    objects, wcs, pointing, bandpass = Input.process(config['input'], logger=logger)
     stars = Select.process(config.get('select',{}), objects, logger=logger)
 
     psf = PSF.process(config['psf'], logger=logger)
-    psf.fit(stars, wcs, pointing, logger=logger)
+    psf.set_context(wcs, pointing, bandpass)
+    psf.fit(stars, logger=logger)
 
     # Attach these for reference
     psf.initial_objects = objects
diff --git a/piff/convolvepsf.py b/piff/convolvepsf.py
index 04568554..3d6970b9 100644
--- a/piff/convolvepsf.py
+++ b/piff/convolvepsf.py
@@ -106,6 +106,13 @@ def interp_property_names(self):
             names.update(c.interp_property_names)
         return names
 
+    def set_context(self, wcs, pointing, bandpass):
+        super().set_context(wcs, pointing, bandpass)
+        if isinstance(self.components, list):
+            # If it's a list, then building components has been completed.
+            for comp in self.components:
+                comp.set_context(wcs, pointing, bandpass)
+
     @classmethod
     def parseKwargs(cls, config_psf, logger):
         """Parse the psf field of a configuration dict and return the kwargs to use for
diff --git a/piff/input.py b/piff/input.py
index d07afabb..6764560d 100644
--- a/piff/input.py
+++ b/piff/input.py
@@ -22,9 +22,11 @@
 import os
 import galsim
 import yaml
+import bisect
 
-from .util import run_multi
+from .util import make_flat, run_multi
 from .star import Star, StarData
+from .config import LoggerWrapper
 
 class Input(object):
     """The base class for handling inputs for building a Piff model.
@@ -40,6 +42,11 @@ class Input(object):
 
     nproc = 1  # Sub-classes can overwrite this as an instance attribute.
 
+    # These shouldn't ever be used, but give these defaults before they are set correctly.
+    wcs_list = []
+    pointing = None
+    bandpass = None
+
     @classmethod
     def process(cls, config_input, logger=None):
         """Parse the input field of the config dict.
@@ -47,11 +54,12 @@ def process(cls, config_input, logger=None):
         :param config_input:    The configuration dict.
         :param logger:          A logger object for logging debug info. [default: None]
 
-        :returns: stars, wcs, pointing
+        :returns: stars, wcs, pointing, bandpass
 
         stars is a list of Star instances with the initial data.
         wcs is a dict of WCS solutions indexed by chipnum.
         pointing is either a galsim.CelestialCoord or None.
+        bandpass is either a galsim.Bandpass or None.
         """
         # Get the class to use for handling the input data
         # Default type is 'Files'
@@ -78,7 +86,10 @@ def process(cls, config_input, logger=None):
         # Get the pointing (the coordinate center of the field of view)
         pointing = input_handler.getPointing(logger)
 
-        return stars, wcs, pointing
+        # Get the observing bandpass, if any
+        bandpass = input_handler.getBandpass(logger)
+
+        return stars, wcs, pointing, bandpass
 
     @classmethod
     def __init_subclass__(cls):
@@ -103,7 +114,6 @@ def makeStars(self, logger=None):
 
         :returns: a list of Star instances
         """
-        from .config import LoggerWrapper
         logger = LoggerWrapper(logger)
 
         if self.nimages == 1:
@@ -148,6 +158,15 @@ def getPointing(self, logger=None):
         """
         return self.pointing
 
+    def getBandpass(self, logger=None):
+        """Get the observing bandpass for the input images.
+
+        :param logger:      A logger object for logging debug info. [default: None]
+
+        :returns: a galsim.Bandpass or None.
+        """
+        return self.bandpass
+
 
 class InputFiles(Input):
     """An Input handler that just takes a list of image files and catalog files.
@@ -258,6 +277,18 @@ class InputFiles(Input):
                         SEDs from FITS tables. [default: ``WAVE``]
         :sed_flux_key:  FITS-table column name to use for flux when loading SEDs
                         from FITS tables. [default: ``FLUX``]
+        :sed_tol:       Relative tolerance for reducing the number of samples of input SEDs
+                        based on the calculated flux over the given bandpass. The SEDs are each
+                        reduced once at read time, and the reduced SED object is then shared
+                        by all matching stars. Values <= 0 disable this step. [default: 0]
+        :sed_max_samples:
+                        Maximum number of samples to keep in the effective SED. If set, use a
+                        greedy knot-selection algorithm to choose up to this many samples, then
+                        refit their values for the best piecewise-linear approximation of the
+                        original effective SED. Values <= 0 disable this step. [default: 0]
+        :bandpass:      GalSim bandpass config dict describing the observing bandpass for the
+                        input image(s). Required when using ``sed_col`` or ``sed_file_name``.
+                        [default: None]
 
     .. note::
 
@@ -387,7 +418,6 @@ class InputFiles(Input):
 
     def __init__(self, config, logger=None):
         import copy
-        from .config import LoggerWrapper
         logger = LoggerWrapper(logger)
 
         req = { 'image_file_name': str,
@@ -410,6 +440,8 @@ def __init__(self, config, logger=None):
                 'sed_flux_type' : str,
                 'sed_wave_key' : str,
                 'sed_flux_key' : str,
+                'sed_tol' : float,
+                'sed_max_samples' : int,
                 'sky_col' : str,
                 'gain_col' : str,
                 'flag_col' : str,
@@ -435,7 +467,7 @@ def __init__(self, config, logger=None):
                 'noise' : str,
                 'nstars' : int,
               }
-        ignore = [ 'nproc', 'nimages', 'ra', 'dec', 'wcs' ]
+        ignore = [ 'nproc', 'nimages', 'ra', 'dec', 'wcs', 'bandpass' ]
 
         # We're going to change the config dict a bit. Make a copy so we don't mess up the
         # user's original dict (in case they care).
@@ -675,6 +707,14 @@ def __init__(self, config, logger=None):
             sed_flux_type = params.get('sed_flux_type', None)
             sed_wave_key = params.get('sed_wave_key', 'WAVE')
             sed_flux_key = params.get('sed_flux_key', 'FLUX')
+            sed_tol = params.get('sed_tol', 0.0)
+            sed_max_samples = params.get('sed_max_samples', 0)
+            if 'bandpass' in config:
+                bandpass = galsim.config.BuildBandpass(config, 'bandpass', base, logger)[0]
+            else:
+                if sed_col is not None or sed_file_name is not None:
+                    raise ValueError("bandpass is required when using sed_col or sed_file_name")
+                bandpass = None
             sky_col = params.get('sky_col', None)
             gain_col = params.get('gain_col', None)
             gain = params.get('gain', None)
@@ -709,6 +749,9 @@ def __init__(self, config, logger=None):
                     'sed_flux_type': sed_flux_type,
                     'sed_wave_key': sed_wave_key,
                     'sed_flux_key': sed_flux_key,
+                    'sed_tol': sed_tol,
+                    'sed_max_samples': sed_max_samples,
+                    'bandpass': bandpass,
                     'sky_col' : sky_col,
                     'gain_col' : gain_col,
                     'sky' : sky,
@@ -719,6 +762,7 @@ def __init__(self, config, logger=None):
                     'stamp_size' : self.stamp_size})
 
         self.use_partial = config.get('use_partial', False)
+        self.bandpass = bandpass
 
         # Read all the wcs's, since we'll need this for the pointing, which in turn we'll
         # need for when we make the stars.
@@ -743,7 +787,6 @@ def load_images(self, stars, logger=None):
 
         :returns: a new list of Stars with the images information loaded.
         """
-        from .config import LoggerWrapper
         logger = LoggerWrapper(logger)
 
         images = {}
@@ -775,7 +818,6 @@ def getRawImageData(self, image_num, logger=None):
     def _getRawImageData(image_kwargs, cat_kwargs, wcs,
                          invert_weight, remove_signal_from_weight,
                          config=None, logger=None):
-        from .config import LoggerWrapper
         logger = LoggerWrapper(logger)
         image, weight = InputFiles.readImage(logger=logger, **image_kwargs)
 
@@ -1124,9 +1166,106 @@ def _parse_flux_type(flux_type):
             except Exception:
                 raise ValueError("Invalid sed_flux_type %r" % (flux_type,))
 
+    @staticmethod
+    def _trapezoid_weights(wave):
+        weights = np.empty_like(wave)
+        weights[0] = 0.5 * (wave[1] - wave[0])
+        weights[-1] = 0.5 * (wave[-1] - wave[-2])
+        weights[1:-1] = 0.5 * (wave[2:] - wave[:-2])
+        return weights
+
+    @staticmethod
+    def _fit_piecewise_linear_values(wave, flux, sample_wave):
+        # Given an original piecewise linear function defined by (wave, flux) pairs and a
+        # (smaller) set of target sample locations, sample_wave, this function solves for flux
+        # values at those sample locations which give a good approximation to the original.
+        #
+        # More specifically, it solves for the sample values whose linearly interpolated
+        # LookupTable is the best weighted least-squares approximation to the original SED.
+        # This is a standard linear least-squares projection onto the space of piecewise-linear
+        # functions spanned by the hat-function basis on sample_wave.
+
+        # Make the design matrix A which gives the new piecewise linear approximation at the
+        # location of each original flux value.
+        # If (w,f) are the original pairs, and (s,g) are the samples, then there is a row in A
+        # for each w, and a column for each s.
+        #    A @ g = f'
+        # gives the approximate flux values, which we want to be close to the original f.
+        n = len(wave)
+        m = len(sample_wave)
+        A = np.zeros((n, m), dtype=float)
+        idx = np.searchsorted(sample_wave, wave, side='right') - 1
+        idx = np.clip(idx, 0, m-2)
+        s0 = sample_wave[idx]
+        s1 = sample_wave[idx+1]
+        # Each dense sample depends only on the two neighboring retained samples.
+        # The coefficient is the barycentric weight of the linear interpolation on the interval.
+        # i.e. the fraction of the interval to each neighboring sample.
+        t = (wave - s0) / (s1 - s0)
+        A[np.arange(n), idx] = 1.0 - t
+        A[np.arange(n), idx+1] = t
+        A[wave == sample_wave[0], :] = 0.0    # These 4 should already be computed correctly,
+        A[wave == sample_wave[0], 0] = 1.0    # but this makes them exact without any possible
+        A[wave == sample_wave[-1], :] = 0.0   # numerical rounding errors.
+        A[wave == sample_wave[-1], -1] = 1.0
+        # Weight by the local wavelength interval so the fit roughly respects the integral
+        # over wavelength rather than just pointwise errors on a nonuniform grid.
+        sqrt_w = np.sqrt(InputFiles._trapezoid_weights(wave))
+        Aw = A * sqrt_w[:, np.newaxis]
+        fw = flux * sqrt_w
+        g = np.linalg.lstsq(Aw, fw, rcond=None)[0]
+
+        # g is now the vector that gives the best least squares approximation f'.
+        return g
+
+    @staticmethod
+    def _select_greedy_samples(wave, flux, candidate_wave, sed_max_samples):
+        # This function selects a subset of candidate_wave locations that would be most
+        # effective for making a piecewise-linear approximation to the SED defined by
+        # (wave, flux).
+
+        if len(candidate_wave) <= sed_max_samples:
+            return candidate_wave
+
+        # Greedy knot placement for piecewise-linear approximation:
+        # keep the endpoints, repeatedly refit the current piecewise-linear model, then add the
+        # candidate wavelength with the largest weighted residual.  This is the same general idea
+        # as orthogonal matching pursuit / greedy basis selection, specialized to a linear spline
+        # basis with fixed candidate knot locations.
+
+        # Start with just the two end points.  We always want these.
+        selected = [0, len(candidate_wave)-1]
+
+        # Use the original dense SED definition for the target values of our fit.
+        # This is piecewise interpolated values using all the candidate waves.
+        target = np.array([np.interp(w, wave, flux) for w in candidate_wave], dtype=float)
+
+        # Get the weights that would be used in the approximation using all candidates.
+        # These will stay fixed throughout.  The selection process will use the combination
+        # of the residual at each point and the weight.
+        cand_weights = InputFiles._trapezoid_weights(candidate_wave)
+
+        # Add candidate waves to our selection list, refitting each time.
+        while len(selected) < sed_max_samples:
+            sample_wave = candidate_wave[selected]
+            sample_flux = InputFiles._fit_piecewise_linear_values(wave, flux, sample_wave)
+            approx = np.interp(candidate_wave, sample_wave, sample_flux)
+
+            # Use weighted residuals so the chosen samples target the integral-relevant
+            # approximation error, not just the largest pointwise deviation on the candidate grid.
+            resid = np.abs(target - approx) * np.sqrt(cand_weights)
+
+            # Add the item with the highest weighted residual
+            resid[selected] = -1.0  # Don't reselect something we already have.
+            bisect.insort(selected, np.argmax(resid))
+
+        return candidate_wave[selected]
+
     @staticmethod
     def _read_sed_file(sed_file_name, sed_wave_type, sed_flux_type,
-                       sed_wave_key, sed_flux_key):
+                       sed_wave_key, sed_flux_key, sed_tol, sed_max_samples, bandpass,
+                       logger):
+        logger = LoggerWrapper(logger)
         cache_key = (sed_file_name, sed_wave_type, sed_flux_type, sed_wave_key, sed_flux_key)
         if cache_key in InputFiles._sed_cache:
             return InputFiles._sed_cache[cache_key]
@@ -1166,6 +1305,46 @@ def _read_sed_file(sed_file_name, sed_wave_type, sed_flux_type,
 
             table = galsim.LookupTable(wave, flux, interpolant='linear')
             sed = galsim.SED(table, wave_type=wave_type, flux_type=flux_type)
+        logger.info("Read sed file %s",sed_file_name)
+
+        # We will only ever use the sed in combination with the bandpass, so we make
+        # that combination here and then possibly thin the result.
+        sed = sed * bandpass
+
+        # This means that the bandpass for drawing needs to be a flat bandpass with the same
+        # limits as the original, which we call flat_bandpass.
+        flat_bandpass = make_flat(bandpass)
+
+        wave = np.array(sed.wave_list, dtype=float)
+        if len(wave) == 0:
+            # If the sed doesn't have a wave_list, we can't thin it.  (This is very rare!)
+            sed = sed.withFlux(1.0, flat_bandpass)
+            InputFiles._sed_cache[cache_key] = sed
+            return sed
+
+        if sed_tol > 0 or sed_max_samples > 0:
+            flux = np.array([sed(w) for w in wave], dtype=float)
+            sample_wave = wave
+            if sed_tol > 0:
+                # sed_tol means use GalSim's thin function.
+                sample_wave = np.array(sed.thin(rel_err=sed_tol).wave_list, dtype=float)
+            if sed_max_samples > 0:
+                # If requested, reduce the number of samples down to the given maximum.
+                if sed_max_samples == 1:
+                    raise ValueError("sed_max_samples must be >= 2")
+                sample_wave = InputFiles._select_greedy_samples(
+                    wave, flux, sample_wave, sed_max_samples
+                )
+
+            # Either version of thinning ends with adjusting the flux values to optimize the
+            # piecewise-linear approximation (the LookupTable).
+            sample_flux = InputFiles._fit_piecewise_linear_values(wave, flux, sample_wave)
+            table = galsim.LookupTable(sample_wave, sample_flux, interpolant='linear')
+            sed = galsim.SED(table, wave_type='nm', flux_type='fphotons')
+            logger.info("  Thinned SED to use %d wave samples", len(sample_wave))
+
+        # Normalize to unit flux.
+        sed = sed.withFlux(1.0, flat_bandpass)
 
         InputFiles._sed_cache[cache_key] = sed
         return sed
@@ -1187,6 +1366,7 @@ def readStarCatalog(cat_file_name, cat_hdu, x_col, y_col,
                         ra_col, dec_col, ra_units, dec_units, image,
                         flag_col, skip_flag, use_flag, property_cols, sed_col, sed_wave_type,
                         sed_file_name, sed_flux_type, sed_wave_key, sed_flux_key,
+                        sed_tol, sed_max_samples, bandpass,
                         properties, image_num, sky_col, gain_col, sky, gain, satur,
                         trust_pos, nstars, stamp_size, config, logger):
         """Read in the star catalogs and return lists of positions for each star in each image.
@@ -1214,6 +1394,9 @@ def readStarCatalog(cat_file_name, cat_hdu, x_col, y_col,
         :param sed_flux_type:   Flux type for SED files.
         :param sed_wave_key:    FITS-table wavelength column name for SED files.
         :param sed_flux_key:    FITS-table flux column name for SED files.
+        :param sed_tol:         Relative tolerance for reducing loaded SEDs.
+        :param sed_max_samples: Maximum number of retained SED samples.
+        :param bandpass:        The galsim.Bandpass for the observation.
         :param sky_col:         A column with sky (background) levels.
         :param gain_col:        A column with gain values.
         :param sky:             Either a float value for the sky to use for all objects or a str
@@ -1369,26 +1552,37 @@ def safe_to_image(wcs, ra, dec):
                     sed_flux_type=sed_flux_type,
                     sed_wave_key=sed_wave_key,
                     sed_flux_key=sed_flux_key,
+                    sed_tol=sed_tol,
+                    sed_max_samples=sed_max_samples,
+                    bandpass=bandpass,
+                    logger=logger,
                 ))
-            extra_props['sed'] = np.array(sed_values, dtype=object)
+            extra_props['sed_eff'] = np.array(sed_values, dtype=object)
             extra_props['sed_file_name'] = np.array(sed_file_name_values, dtype=object)
             extra_props['sed_wave_type'] = np.array([sed_wave_type] * len(cat), dtype=object)
             extra_props['sed_flux_type'] = np.array([sed_flux_type] * len(cat), dtype=object)
             extra_props['sed_wave_key'] = np.array([sed_wave_key] * len(cat), dtype=object)
             extra_props['sed_flux_key'] = np.array([sed_flux_key] * len(cat), dtype=object)
+            extra_props['sed_tol'] = np.full(len(cat), sed_tol, dtype=float)
+            extra_props['sed_max_samples'] = np.full(len(cat), sed_max_samples, dtype=int)
         elif sed_file_name is not None:
             sed_file_name = InputFiles._resolve_sed_file_name(sed_file_name, cat_file_name)
             sed = InputFiles._read_sed_file(
                 sed_file_name, sed_wave_type=sed_wave_type,
                 sed_flux_type=sed_flux_type, sed_wave_key=sed_wave_key,
-                sed_flux_key=sed_flux_key
+                sed_flux_key=sed_flux_key, sed_tol=sed_tol,
+                sed_max_samples=sed_max_samples,
+                bandpass=bandpass,
+                logger=logger,
             )
-            extra_props['sed'] = np.array([sed] * len(cat), dtype=object)
+            extra_props['sed_eff'] = np.array([sed] * len(cat), dtype=object)
             extra_props['sed_file_name'] = np.array([sed_file_name] * len(cat), dtype=object)
             extra_props['sed_wave_type'] = np.array([sed_wave_type] * len(cat), dtype=object)
             extra_props['sed_flux_type'] = np.array([sed_flux_type] * len(cat), dtype=object)
             extra_props['sed_wave_key'] = np.array([sed_wave_key] * len(cat), dtype=object)
             extra_props['sed_flux_key'] = np.array([sed_flux_key] * len(cat), dtype=object)
+            extra_props['sed_tol'] = np.full(len(cat), sed_tol, dtype=float)
+            extra_props['sed_max_samples'] = np.full(len(cat), sed_max_samples, dtype=int)
 
         # If we used a flag column, keep it as a property.
         if flag_col is not None:
@@ -1504,7 +1698,6 @@ def setPointing(self, ra, dec, logger=None):
            individual WCS functions. [Not implemented currently.]
         """
         import fitsio
-        from .config import LoggerWrapper
         logger = LoggerWrapper(logger)
 
         if (ra is None) != (dec is None):
diff --git a/piff/optical_model.py b/piff/optical_model.py
index 2e745061..efc50687 100644
--- a/piff/optical_model.py
+++ b/piff/optical_model.py
@@ -165,7 +165,7 @@ class Optical(Model):
     """
     _type_name = 'Optical'
     _method = 'auto'
-    _centered = False
+    _centered = True
     _model_can_be_offset = False
 
     def __init__(self, template=None, gsparams=None, atmo_type='VonKarman', logger=None, **kwargs):
diff --git a/piff/psf.py b/piff/psf.py
index 6dd00cef..668da935 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -20,6 +20,7 @@
 import galsim
 
 from .star import Star, StarData
+from .util import make_flat
 
 class PSF(object):
     """The base class for describing a PSF model across a field of view.
@@ -41,6 +42,10 @@ class PSF(object):
     # Normally overridden by subclasses.  But this gives a reasonable default for
     # tests that bypass places where it would be potentially set in normal operations.
     degenerate_points = False
+    wcs = None
+    pointing = None
+    bandpass = None
+    flat_bandpass = None
 
     @classmethod
     def process(cls, config_psf, logger=None):
@@ -54,14 +59,21 @@ def process(cls, config_psf, logger=None):
         several components.
 
         This function merely creates a "blank" PSF object.  It does not actually do any
-        part of the solution yet.  Typically this will be followed by fit:
+        part of the solution yet.  Typically this will be followed by set_context and fit:
 
+            >>> stars, wcs, pointing, bandpass = piff.Input.process(config['input'])
             >>> psf = piff.PSF.process(config['psf'])
-            >>> stars, wcs, pointing = piff.Input.process(config['input'])
-            >>> psf.fit(stars, wcs, pointing)
+            >>> psf.set_context(wcs, pointing, bandpass)
+            >>> psf.fit(stars)
 
         at which point, the ``psf`` instance would have a solution to the PSF model.
 
+        .. note::
+
+            Shared runtime context such as wcs, pointing, and bandpass should be attached
+            with set_context before calling fit.  The old pattern of providing wcs and
+            pointing to fit is still supported, but deprecated.
+
         :param config_psf:  A dict specifying what type of PSF to build along with the
                             appropriate kwargs for building it.
         :param logger:      A logger object for logging debug info. [default: None]
@@ -95,6 +107,25 @@ def process(cls, config_psf, logger=None):
 
         return psf
 
+    def set_context(self, wcs, pointing=None, bandpass=None):
+        """Attach shared runtime context to this PSF instance.
+
+        It is permissible for pointing to be None if the WCS is not Celestial.
+        Also, bandpass may be None if the psf to be fit is not chromatic.
+        If pointing or bandpass is None, the existing value is left unchanged.
+
+        :param wcs:         A dict of WCS solutions indexed by chipnum.
+        :param pointing:    A galsim.CelestialCoord object giving the telescope pointing.
+                            [default: None]
+        :param bandpass:    A galsim.Bandpass giving the observing bandpass. [default: None]
+        """
+        self.wcs = wcs
+        if pointing is not None:
+            self.pointing = pointing
+        if bandpass is not None:
+            self.bandpass = bandpass
+            self.flat_bandpass = make_flat(bandpass)
+
     def set_num(self, num):
         """If there are multiple components involved in the fit, set the number to use
         for this model.
@@ -261,22 +292,27 @@ def reflux(self, star, logger=None):
 
         if self.fit_center and not star.data.properties.get('trust_pos',False):
             psf_prof, method = self._getRawProfile(model_star)
+            if isinstance(psf_prof, galsim.ChromaticObject):
+                psf_prof = psf_prof * star.data.properties['sed_eff']
+                draw_args = (self.flat_bandpass,)
+            else:
+                draw_args = ()
 
             # Use finite different to approximate d(model)/duc, d(model)/dvc
             duv = 1.e-5
             temp = star.image.copy()
             center = star.fit.center
             du_prof = psf_prof.shift(center[0] + duv, center[1]) * new_flux
-            du_prof.drawImage(temp, method=method, center=star.image_pos)
+            du_prof.drawImage(*draw_args, image=temp, method=method, center=star.image_pos)
             dmduc = (temp.array.ravel() - model) / duv
             dv_prof = psf_prof.shift(center[0], center[1] + duv) * new_flux
-            dv_prof.drawImage(temp, method=method, center=star.image_pos)
+            dv_prof.drawImage(*draw_args, image=temp, method=method, center=star.image_pos)
             dmdvc = (temp.array.ravel() - model) / duv
 
             # Also dmdflux
             dflux = 1.e-5 * max(abs(new_flux), 1.e-5)  # Guard against division by 0
             df_prof = psf_prof.shift(center[0], center[1]) * (new_flux + dflux)
-            df_prof.drawImage(temp, method=method, center=star.image_pos)
+            df_prof.drawImage(*draw_args, image=temp, method=method, center=star.image_pos)
             dmdf = (temp.array.ravel() - model) / dflux
 
             # Now construct the design matrix for this minimization
@@ -301,7 +337,7 @@ def reflux(self, star, logger=None):
             # Extract the values we want.
             df, duc, dvc = x
 
-            if self.include_model_centroid and psf_prof.centroid != galsim.PositionD(0,0):
+            if self.include_model_centroid:
                 # In addition to shifting to the best fit center location, also shift
                 # by the centroid of the model itself, so the next next pass through the
                 # fit will be closer to centered.  In practice, this converges pretty quickly.
@@ -373,13 +409,14 @@ def remove_outliers(self, stars, iteration, logger):
             nremoved = 0
         return stars, nremoved
 
-    def fit(self, stars, wcs, pointing, logger=None, convert_funcs=None, draw_method=None):
+    def fit(self, stars, wcs=None, pointing=None, logger=None,
+            convert_funcs=None, draw_method=None):
         """Fit interpolated PSF model to star data using standard sequence of operations.
 
         :param stars:           A list of Star instances.
-        :param wcs:             A dict of WCS solutions indexed by chipnum.
+        :param wcs:             A dict of WCS solutions indexed by chipnum. [deprecated]
         :param pointing:        A galsim.CelestialCoord object giving the telescope pointing.
-                                [Note: pointing should be None if the WCS is not a CelestialWCS]
+                                [deprecated]
         :param logger:          A logger object for logging debug info. [default: None]
         :param convert_funcs:   An optional list of function to apply to the profiles being fit
                                 before drawing it onto the image.  This is used by composite PSFs
@@ -391,8 +428,9 @@ def fit(self, stars, wcs, pointing, logger=None, convert_funcs=None, draw_method
         from .config import LoggerWrapper
         logger = LoggerWrapper(logger)
 
-        self.wcs = wcs
-        self.pointing = pointing
+        if self.wcs is None:
+            logger.error("WARNING: wcs should now be set with set_context before calling fit.")
+            self.set_context(wcs, pointing, self.bandpass)
 
         # Initialize stars as needed by the PSF modeling class.
         stars = self.initialize_flux_center(stars, logger=logger)
@@ -766,9 +804,12 @@ def _write(self, writer, name, logger):
             if hasattr(self, 'stars'):
                 Star.write(self.stars, w, 'stars')
                 logger.verbose("Wrote the PSF stars to name %s", w.get_full_name('stars'))
-            if hasattr(self, 'wcs'):
+            if self.wcs is not None:
                 w.write_wcs_map('wcs', self.wcs, self.pointing)
                 logger.verbose("Wrote the PSF WCS to name %s", w.get_full_name('wcs'))
+            if self.bandpass is not None:
+                w.write_bandpass('bandpass', self.bandpass)
+                logger.verbose("Wrote the PSF bandpass to name %s", w.get_full_name('bandpass'))
             self._finish_write(w, logger=logger)
 
     @classmethod
@@ -822,15 +863,15 @@ def _read(cls, reader, name, logger):
 
         with reader.nested(name) as r:
             # Read the stars, wcs, pointing values
-            stars = Star.read(r, 'stars')
+            wcs, pointing = r.read_wcs_map('wcs', logger=logger)
+            bandpass = r.read_bandpass('bandpass')
+            stars = Star.read(r, 'stars', bandpass=bandpass)
             if stars is not None:
                 logger.debug("stars = %s", stars)
                 psf.stars = stars
-            wcs, pointing = r.read_wcs_map('wcs', logger=logger)
             if wcs is not None:
-                logger.debug("wcs = %s, pointing = %s",wcs,pointing)
-                psf.wcs = wcs
-                psf.pointing = pointing
+                logger.debug("wcs = %s, pointing = %s, bandpass = %s", wcs, pointing, bandpass)
+                psf.set_context(wcs, pointing, bandpass)
 
             # Just in case the class needs to do something else at the end.
             psf._finish_read(r, logger)
diff --git a/piff/readers.py b/piff/readers.py
index 8d9a2b27..68b5bada 100644
--- a/piff/readers.py
+++ b/piff/readers.py
@@ -21,6 +21,21 @@
 import fitsio
 import galsim
 import copy
+import numpy as np
+
+def _read_chunked_objs(data, prefix, encoding='bytes'):
+    import base64
+    import pickle
+    nchunks_key = f'{prefix}_nchunks'
+    if nchunks_key in data.dtype.names:
+        nchunks = data[nchunks_key][0]
+    else:
+        # Prior to version 1.7, there was only 1 chunked item, which is now wcs_nchunks.
+        nchunks = data['nchunks'][0]
+    keys = [f'{prefix}_str_{i:04d}' for i in range(nchunks)]
+    serialized_list = [b''.join(s) for s in zip(*[data[key] for key in keys])]
+    serialized_list = [base64.b64decode(s) for s in serialized_list]
+    return [pickle.loads(s, encoding=encoding) for s in serialized_list]
 
 
 class FitsReader:
@@ -108,40 +123,21 @@ def read_wcs_map(self, name, logger):
             return None, None
 
         assert 'chipnums' in self._fits[extname].get_colnames()
-        assert 'nchunks' in self._fits[extname].get_colnames()
+        assert ('wcs_nchunks' in self._fits[extname].get_colnames() or
+                'nchunks' in self._fits[extname].get_colnames())
 
         data = self._fits[extname].read()
 
         chipnums = data['chipnums']
-        nchunks = data['nchunks']
-        nchunks = nchunks[0]  # These are all equal, so just take first one.
-
-        wcs_keys = [ 'wcs_str_%04d'%i for i in range(nchunks) ]
-        wcs_str = [ data[key] for key in wcs_keys ] # Get all wcs_str columns
-        try:
-            wcs_str = [ b''.join(s) for s in zip(*wcs_str) ]  # Rejoint into single string each
-        except TypeError:  # pragma: no cover
-            # fitsio 1.0 returns strings
-            wcs_str = [ ''.join(s) for s in zip(*wcs_str) ]  # Rejoint into single string each
-
-        wcs_str = [ base64.b64decode(s) for s in wcs_str ] # Convert back from b64 encoding
-        # Convert back into wcs objects
-        try:
-            wcs_list = [ pickle.loads(s, encoding='bytes') for s in wcs_str ]
-        except Exception:
-            # If the file was written by py2, the bytes encoding might raise here,
-            # or it might not until we try to use it.
-            wcs_list = [ pickle.loads(s, encoding='latin1') for s in wcs_str ]
-
-        wcs = dict(zip(chipnums, wcs_list))
-
         try:
+            wcs_list = _read_chunked_objs(data, 'wcs')
+            wcs = dict(zip(chipnums, wcs_list))
             # If this doesn't work, then the file was probably written by py2, not py3
             repr(wcs)
         except Exception:
             logger.verbose('Failed to decode wcs with bytes encoding.')
             logger.verbose('Retry with encoding="latin1" in case file written with python 2.')
-            wcs_list = [ pickle.loads(s, encoding='latin1') for s in wcs_str ]
+            wcs_list = _read_chunked_objs(data, 'wcs', encoding='latin1')
             wcs = dict(zip(chipnums, wcs_list))
             repr(wcs)
 
@@ -169,6 +165,13 @@ def read_wcs_map(self, name, logger):
 
         return wcs, pointing
 
+    def read_bandpass(self, name):
+        extname = self.get_full_name(name)
+        if extname not in self._fits:
+            return None
+        data = self._fits[extname].read()
+        return _read_chunked_objs(data, 'bandpass')[0]
+
     @contextmanager
     def nested(self, name):
         """Return a context manager that yields a new `FitsReader` that reads
diff --git a/piff/roman/roman_psf.py b/piff/roman/roman_psf.py
index 834dc18f..256179ad 100644
--- a/piff/roman/roman_psf.py
+++ b/piff/roman/roman_psf.py
@@ -28,6 +28,7 @@
 from ..psf import PSF
 from ..star import Star
 from ..config import LoggerWrapper
+from ..util import make_flat
 from ..util import run_multi
 
 # Global control of GalSim Roman pupil-plane resolution in getPSF calls.
@@ -145,13 +146,13 @@ class RomanOpticalModel(Model):
 
     _type_name = 'Roman'
     _method = 'auto'
-    _centered = False
+    _centered = True
     _model_can_be_offset = False
 
     def __init__(
         self,
-        filter,
         chromatic=True,
+        filter_name=None,
         max_zernike=22,
         aberration_interp='constant',
         nominal_interp='bilinear',
@@ -160,12 +161,14 @@ def __init__(
         logger=None,
     ):
         self.logger = logger
-        self.filter = filter
         self.chromatic = chromatic
         self.max_zernike = int(max_zernike)
         self.aberration_interp = str(aberration_interp)
         self.nominal_interp = str(nominal_interp)
         self.nproc = int(nproc)
+        self.bandpass = None
+        self.flat_bandpass = None
+        self.filter_name = filter_name
         self.set_num(None)
 
         if self.max_zernike < 4 or self.max_zernike > 22:
@@ -174,21 +177,14 @@ def __init__(
             raise ValueError("aberration_interp must be one of 'global', 'constant', 'linear'")
         if self.nominal_interp not in ('bilinear', 'five_point'):
             raise ValueError("nominal_interp must be one of 'bilinear', 'five_point'")
-
-        # Notation: filter is the string name of the filter.
-        #           bandpass is the galsim.Bandpass object with the transmission function.
-        bandpasses = galsim.roman.getBandpasses()
-        if self.filter not in bandpasses:
-            raise ValueError("Roman filter %r is not a valid GalSim Roman bandpass" % self.filter)
-        self.bandpass = bandpasses[self.filter]
         # Save the original for serialization.
         self.orig_prior_sigma = self._parse_prior_sigma(aberration_prior_sigma)
         self.prior_sigma = self._expand_prior_sigma(self.orig_prior_sigma)
         self.prior_invsigsq = 1.0/self.prior_sigma**2 if self.prior_sigma is not None else None
 
         self.kwargs = {
-            'filter': self.filter,
             'chromatic': self.chromatic,
+            'filter_name': self.filter_name,
             'max_zernike': self.max_zernike,
             'aberration_interp': self.aberration_interp,
             'nominal_interp': self.nominal_interp,
@@ -223,6 +219,32 @@ def initialize_iteration(self):
     def clear_cache(self):
         self._corner_cache = {}
 
+    def set_bandpass(self, bandpass, filter_name=None):
+        if self.filter_name is None:
+            if filter_name is not None:
+                self.filter_name = filter_name
+            elif getattr(bandpass, 'name', None) is not None:
+                self.filter_name = bandpass.name
+            else:
+                raise ValueError(
+                    "RomanOptics requires filter_name when bandpass has no name attribute"
+                )
+        self.bandpass = bandpass
+        self.flat_bandpass = make_flat(bandpass)
+
+    def _require_bandpass(self):
+        if self.bandpass is None:
+            raise ValueError("RomanOptics requires bandpass to be set before use")
+        return self.bandpass
+
+    def _get_sed_eff(self, star):
+        sed_eff = star.data.properties.get('sed_eff')
+        if sed_eff is None:
+            raise ValueError(
+                "RomanOptics with chromatic=True requires each star to have an 'sed_eff' property"
+            )
+        return sed_eff
+
     def _get_roman_five_point_data(self, sca):
         if sca in self._roman_five_point_data:
             return self._roman_five_point_data[sca]
@@ -285,13 +307,12 @@ def _make_extra_aberrations(self, params):
         extra_aberrations[4:] = params
         return extra_aberrations
 
-    def _draw_profile_to_image(self, prof, image, center):
-        prof.drawImage(
-            image,
-            method=self._method,
-            center=center,
-            bandpass=self.bandpass if self.chromatic else None
-        )
+    def _draw_profile_to_image(self, prof, image, center, star):
+        if self.chromatic:
+            prof = prof * self._get_sed_eff(star)
+            prof.drawImage(self.flat_bandpass, image=image, method=self._method, center=center)
+        else:
+            prof.drawImage(image=image, method=self._method, center=center)
 
     def _solve_params(self, aw, bw, params):
         ata = aw.T.dot(aw)
@@ -449,10 +470,7 @@ def _fit_sca_group(self, sca, stars, convert_funcs, logger=None, draw_method=Non
             image = star.image.array
             weight = star.weight.array
             base_image = self._draw_model_image_from_samples(
-                star,
-                sample_base_profiles,
-                convert_func=convert_func,
-            )
+                    star, sample_base_profiles, convert_func)
             resid = image - base_image
             sw = np.sqrt(weight.ravel())
             bw = resid.ravel() * sw
@@ -468,11 +486,7 @@ def _fit_sca_group(self, sca, stars, convert_funcs, logger=None, draw_method=Non
             scale = 1.0 / step
             # Use these adjusted sample profiles for all stars on the SCA.
             for star, convert_func, _, _, _, jac, base_image in group_data:
-                im1 = self._draw_model_image_from_samples(
-                    star,
-                    sample_plus_profiles,
-                    convert_func=convert_func,
-                )
+                im1 = self._draw_model_image_from_samples(star, sample_plus_profiles, convert_func)
                 jac[:, i] = ((im1 - base_image) * scale).ravel()
 
         # Build one SCA-level normal equation from all stars.
@@ -493,8 +507,7 @@ def _fit_sca_group(self, sca, stars, convert_funcs, logger=None, draw_method=Non
         out = []
         for star, convert_func, weight in solved:
             model_image = self._draw_model_image_from_samples(
-                star, sample_sca_profiles, convert_func=convert_func
-            )
+                    star, sample_sca_profiles, convert_func)
             chisq = np.sum(weight * (star.image.array - model_image) ** 2)
             dof = np.count_nonzero(weight) - 3
             out.append(
@@ -511,7 +524,7 @@ def draw(self, star, copy_image=True):
         params = star.fit.get_params(self._num)
         prof = self.getProfile(params, star=star).shift(star.fit.center) * star.fit.flux
         image = star.image.copy() if copy_image else star.image
-        self._draw_profile_to_image(prof, image, center=star.image_pos)
+        self._draw_profile_to_image(prof, image, star.image_pos, star)
         return Star(star.data.withNew(image=image), star.fit)
 
     def getProfile(self, params=None, star=None, cache=True):
@@ -528,13 +541,13 @@ def getProfile(self, params=None, star=None, cache=True):
         sample_profiles = self._get_sample_profiles(star, params, cache=cache)
         return self._interpolate_samples(star, sample_profiles)
 
-    def _draw_model_image_from_samples(self, star, sample_profiles, convert_func=None):
+    def _draw_model_image_from_samples(self, star, sample_profiles, convert_func):
         prof = self._interpolate_samples(star, sample_profiles)
         prof = prof.shift(star.fit.center) * star.fit.flux
         if convert_func is not None:
             prof = convert_func(prof)
         image = star.image.copy()
-        self._draw_profile_to_image(prof, image, star.image_pos)
+        self._draw_profile_to_image(prof, image, star.image_pos, star)
         return image.array
 
     def _get_sample_profiles(self, star, params, cache=True, sca=None):
@@ -549,7 +562,8 @@ def _get_sample_profiles(self, star, params, cache=True, sca=None):
             if same_wcs and np.array_equal(cached_params, params):
                 return cached_profiles
 
-        wavelength = None if self.chromatic else self.bandpass.effective_wavelength
+        bandpass = self._require_bandpass()
+        wavelength = None if self.chromatic else bandpass.effective_wavelength
         if self.nominal_interp == 'five_point':
             points = self._get_roman_five_point_data(sca)['points']
         else:
@@ -573,7 +587,7 @@ def _get_sample_profiles(self, star, params, cache=True, sca=None):
             profiles.append(
                 galsim.roman.getPSF(
                     sca,
-                    self.filter,
+                    self.filter_name,
                     SCA_pos=pt,
                     pupil_bin=pupil_bin,
                     wcs=wcs,
@@ -627,6 +641,11 @@ class RomanOpticsPSF(PSF):
     def __init__(self, model, interp, outliers=None, chisq_thresh=0.1, min_iter=2, max_iter=30):
         self.model = model
         self.interp = interp
+        if isinstance(model, Model):
+            self.bandpass = model.bandpass
+        else:
+            # During read, model is temporarily a placeholder integer until _finish_read.
+            self.bandpass = None
         self.outliers = outliers
         self.chisq_thresh = chisq_thresh
         self.min_iter = min_iter
@@ -653,9 +672,14 @@ def set_num(self, num):
         # Only call set_num if they are actually built.
         if isinstance(self.model, Model):
             self.model.set_num(num)
-        if isinstance(self.interp, Interp):
             self.interp.set_num(num)
 
+    def set_context(self, wcs, pointing=None, bandpass=None):
+        super().set_context(wcs, pointing, bandpass)
+        # During read, model is temporarily a placeholder integer until _finish_read.
+        if isinstance(self.model, Model):
+            self.model.set_bandpass(self.bandpass)
+
     @property
     def interp_property_names(self):
         return self.interp.property_names
@@ -782,5 +806,8 @@ def _finish_read(self, reader, logger):
         for key in chisq_dict:
             setattr(self, key, chisq_dict[key])
         self.model = Model.read(reader, 'model')
+        if self.bandpass is not None:
+            # Now that model is built, we can update its bandpass to the correct value.
+            self.model.set_bandpass(self.bandpass)
         self.interp = Interp.read(reader, 'interp')
         self.outliers = Outliers.read_all(reader, 'outliers')
diff --git a/piff/singlechip.py b/piff/singlechip.py
index ce1d4f57..990e961f 100644
--- a/piff/singlechip.py
+++ b/piff/singlechip.py
@@ -24,18 +24,20 @@
 from .util import make_dtype, adjust_value, run_multi
 
 # Used by SingleChipPSF.fit
-def single_chip_run(chipnum, single_psf, stars, wcs, pointing, convert_funcs, draw_method, logger):
+def single_chip_run(chipnum, single_psf, stars, wcs, convert_funcs, draw_method, logger):
     # Make a copy of single_psf for each chip
     psf_chip = copy.deepcopy(single_psf)
 
     # Break the list of stars up into a list for each chip
     stars_chip = [ s for s in stars if s['chipnum'] == chipnum ]
+
+    # Set the wcs correctly for this chip
     wcs_chip = { chipnum : wcs[chipnum] }
+    psf_chip.set_context(wcs_chip)
 
     # Run the psf_chip fit function using this stars and wcs (and the same pointing)
     logger.info("Building solution for chip %s with %d stars", chipnum, len(stars_chip))
-    psf_chip.fit(stars_chip, wcs_chip, pointing, logger=logger, convert_funcs=convert_funcs,
-                 draw_method=draw_method)
+    psf_chip.fit(stars_chip, logger=logger, convert_funcs=convert_funcs, draw_method=draw_method)
 
     return psf_chip
 
@@ -69,6 +71,12 @@ def set_num(self, num):
         if isinstance(self.single_psf, PSF):
             self.single_psf.set_num(num)
 
+    def set_context(self, wcs, pointing=None, bandpass=None):
+        super().set_context(wcs, pointing, bandpass)
+        # wcs will be different for each chip.  This will be set in single_chip_run.
+        if isinstance(self.single_psf, PSF):
+            self.single_psf.set_context(None, pointing, bandpass)
+
     @property
     def interp_property_names(self):
         return self.single_psf.interp_property_names
@@ -91,17 +99,18 @@ def parseKwargs(cls, config_psf, logger):
         config_psf['type'] = config_psf.pop('single_type', 'Simple')
 
         # Now the regular PSF process function can process the dict.
-        single_psf = PSF.process(config_psf, logger)
+        single_psf = PSF.process(config_psf, logger=logger)
 
         return { 'single_psf' : single_psf, 'nproc' : nproc }
 
-    def fit(self, stars, wcs, pointing, logger=None, convert_funcs=None, draw_method=None):
+    def fit(self, stars, wcs=None, pointing=None, logger=None, convert_funcs=None,
+            draw_method=None):
         """Fit interpolated PSF model to star data using standard sequence of operations.
 
         :param stars:           A list of Star instances.
-        :param wcs:             A dict of WCS solutions indexed by chipnum.
+        :param wcs:             A dict of WCS solutions indexed by chipnum. [deprecated]
         :param pointing:        A galsim.CelestialCoord object giving the telescope pointing.
-                                [Note: pointing should be None if the WCS is not a CelestialWCS]
+                                [deprecated]
         :param logger:          A logger object for logging debug info. [default: None]
         :param convert_funcs:   An optional list of function to apply to the profiles being fit
                                 before drawing it onto the image.  This is used by composite PSFs
@@ -112,12 +121,13 @@ def fit(self, stars, wcs, pointing, logger=None, convert_funcs=None, draw_method
         """
         from .config import LoggerWrapper
         logger = LoggerWrapper(logger)
-        self.wcs = wcs
-        self.pointing = pointing
+        if self.wcs is None:
+            logger.error("WARNING: wcs should now be set with set_context before calling fit.")
+            self.set_context(wcs, pointing)
         self.psf_by_chip = {}
 
-        chipnums = list(wcs.keys())
-        args = [(chipnum, self.single_psf, stars, wcs, pointing, convert_funcs, draw_method)
+        chipnums = list(self.wcs.keys())
+        args = [(chipnum, self.single_psf, stars, self.wcs, convert_funcs, draw_method)
                 for chipnum in chipnums]
 
         output = run_multi(single_chip_run, self.nproc, raise_except=False,
diff --git a/piff/star.py b/piff/star.py
index dca4052b..99d9703a 100644
--- a/piff/star.py
+++ b/piff/star.py
@@ -34,7 +34,7 @@ class Star(object):
     Piff functions such as:
 
         stars = input_handler.makeStars(logger)
-        stars, wcs = piff.Input.process(config['input'])
+        stars, wcs, pointing, bandpass = piff.Input.process(config['input'])
         target_star = piff.Star.makeTarget(x=x, y=y, wcs=wcs, ...)
 
     However, a star can be constructed directly from a StarData instance and a StarFit instance.
@@ -347,7 +347,7 @@ def write(cls, stars, writer, name):
         # Start with the data properties
         prop_keys = list(stars[0].data.properties)
         # SED objects are large and not directly serializable; persist SED metadata only.
-        prop_keys = [key for key in prop_keys if key != 'sed']
+        prop_keys = [key for key in prop_keys if key != 'sed_eff']
         # Do the position ones first
         for key in [ 'x', 'y', 'u', 'v' ]:
             dtypes.append( (key, float) )
@@ -438,11 +438,12 @@ def read_coords_params(cls, fits, extname):
         return coords, params
 
     @classmethod
-    def read(cls, reader, name):
+    def read(cls, reader, name, bandpass=None):
         """Read a list of stars via an open reader object.
 
         :param reader:      A reader object that encapsulates the serialization format.
-        :param name:         Name associated with the stars in the serialized output.
+        :param name:        Name associated with the stars in the serialized output.
+        :param bandpass:    galsim.Bandpass to use when reconstructing stored SEDs.
         :returns: a list of Star instances, or None if there aren't any
         """
         metadata = {}
@@ -523,7 +524,7 @@ def read(cls, reader, name):
                 prop_list[i]['is_reserve'] = True
             if 'flag_psf' in data.dtype.names and row['flag_psf']:
                 prop_list[i]['is_flagged'] = True
-        cls._restore_seds(prop_list)
+        cls._restore_seds(prop_list, bandpass)
 
         wcs_list = [ galsim.JacobianWCS(*jac) for jac in zip(dudx,dudy,dvdx,dvdy) ]
         pos_list = [ galsim.PositionD(*pos) for pos in zip(x_list,y_list) ]
@@ -541,23 +542,31 @@ def read(cls, reader, name):
         return stars
 
     @classmethod
-    def _restore_seds(cls, prop_list):
+    def _restore_seds(cls, prop_list, bandpass):
         from .input import InputFiles
 
         for prop in prop_list:
-            if 'sed' in prop or 'sed_file_name' not in prop:
+            if 'sed_eff' in prop or 'sed_file_name' not in prop:
                 continue
+            if bandpass is None:
+                raise ValueError("bandpass is required to restore stars with serialized SEDs")
             sed_file_name = prop['sed_file_name']
             sed_wave_key = prop.get('sed_wave_key', 'WAVE')
             sed_flux_key = prop.get('sed_flux_key', 'FLUX')
             sed_wave_type = prop.get('sed_wave_type')
             sed_flux_type = prop.get('sed_flux_type')
-            prop['sed'] = InputFiles._read_sed_file(
+            sed_tol = prop.get('sed_tol', 0.0)
+            sed_max_samples = prop.get('sed_max_samples', 0)
+            prop['sed_eff'] = InputFiles._read_sed_file(
                 sed_file_name,
                 sed_wave_type=sed_wave_type,
                 sed_flux_type=sed_flux_type,
                 sed_wave_key=sed_wave_key,
                 sed_flux_key=sed_flux_key,
+                sed_tol=sed_tol,
+                sed_max_samples=sed_max_samples,
+                bandpass=bandpass,
+                logger=None
             )
 
     @staticmethod
diff --git a/piff/sumpsf.py b/piff/sumpsf.py
index 5c210e1e..6a5ccbe4 100644
--- a/piff/sumpsf.py
+++ b/piff/sumpsf.py
@@ -106,6 +106,12 @@ def interp_property_names(self):
             names.update(c.interp_property_names)
         return names
 
+    def set_context(self, wcs, pointing, bandpass):
+        super().set_context(wcs, pointing, bandpass)
+        if isinstance(self.components, list):
+            for comp in self.components:
+                comp.set_context(wcs, pointing, bandpass)
+
     @classmethod
     def parseKwargs(cls, config_psf, logger):
         """Parse the psf field of a configuration dict and return the kwargs to use for
@@ -298,7 +304,7 @@ def _getRawProfile(self, star):
             method = methods[0]
 
         # Add them up.
-        return galsim.Sum(profiles), method
+        return galsim.Add(profiles), method
 
     def _finish_write(self, writer, logger):
         """Finish the writing process with any class-specific steps.
diff --git a/piff/util.py b/piff/util.py
index 3d970fca..b1581b00 100644
--- a/piff/util.py
+++ b/piff/util.py
@@ -33,6 +33,14 @@ def ensure_dir(target):
     if d != '' and  not os.path.exists(d):
         os.makedirs(d)
 
+
+def make_flat(bandpass):
+    """Return a flat bandpass with the same limits as ``bandpass``."""
+    blue = bandpass.blue_limit
+    red = bandpass.red_limit
+    tab = galsim.LookupTable([blue, red], [1, 1], interpolant='linear')
+    return galsim.Bandpass(tab, 'nm', blue_limit=blue, red_limit=red)
+
 def make_dtype(key, value):
     """A helper function that makes a dtype appropriate for a given value
 
diff --git a/piff/writers.py b/piff/writers.py
index fa527ce6..7941b92b 100644
--- a/piff/writers.py
+++ b/piff/writers.py
@@ -25,6 +25,33 @@
 from .util import make_dtype, adjust_value
 
 
+def _make_chunked_cols(objs, prefix):
+    # Some serialized payloads are too long for a single FITS table column
+    # (max width 28799). Split them into 2**14-sized chunks, which is safely
+    # below that limit. This is needed for some GalSim WCS objects (notably
+    # Pixmappy) and also for serialized Bandpass objects.
+
+    import base64
+    import pickle
+
+    serialized_list = [base64.b64encode(pickle.dumps(obj)) for obj in objs]
+    nrows = len(serialized_list)
+    max_len = max(len(s) for s in serialized_list)
+    chunk_size = 2**14
+    nchunks = max_len // chunk_size + 1
+    cols = [[nchunks] * nrows]
+    dtypes = [(f"{prefix}_nchunks", int)]
+
+    chunk_size = (max_len + nchunks - 1) // nchunks
+    chunks = [
+        [s[i : i + chunk_size] for i in range(0, max_len, chunk_size)]
+        for s in serialized_list
+    ]
+    cols.extend(zip(*chunks))
+    dtypes.extend((f"{prefix}_str_{i:04d}", bytes, chunk_size) for i in range(nchunks))
+    return cols, dtypes
+
+
 class FitsWriter:
     """A writer object that writes to multiple FITS HDUs.
 
@@ -96,40 +123,20 @@ def write_table(self, name, array, metadata=None):
         self._fits.write_table(array, extname=self.get_full_name(name), header=header)
 
     def write_wcs_map(self, name, wcs_map, pointing):
-        """Write a regular a map of WCS objects and an optoinal pointing coord.
+        """Write a regular a map of WCS objects and an optional pointing coord.
 
         :param name:      Name used to save this WCS map.
         :param wcs_map:   A `dict` mapping `int` chipnum to `galsim.BaseWCS`.
         :param pointing:  A `galsim.CelestialCoord`, or `None`.
-        :param logger:    A logger object for logging debug info.
         """
-        import base64
-        import pickle
         # Start with the chipnums
         chipnums = list(wcs_map.keys())
         cols = [chipnums]
         dtypes = [("chipnums", int)]
 
-        # GalSim WCS objects can be serialized via pickle
-        wcs_str = [base64.b64encode(pickle.dumps(w)) for w in wcs_map.values()]
-        max_len = np.max([len(s) for s in wcs_str])
-        # Some GalSim WCS serializations are rather long.  In particular, the Pixmappy one
-        # is longer than the maximum length allowed for a column in a fits table (28799).
-        # So split it into chunks of size 2**14 (mildly less than this maximum).
-        chunk_size = 2**14
-        nchunks = max_len // chunk_size + 1
-        cols.append([nchunks] * len(chipnums))
-        dtypes.append(("nchunks", int))
-
-        # Update to size of chunk we actually need.
-        chunk_size = (max_len + nchunks - 1) // nchunks
-
-        chunks = [
-            [s[i : i + chunk_size] for i in range(0, max_len, chunk_size)]
-            for s in wcs_str
-        ]
-        cols.extend(zip(*chunks))
-        dtypes.extend(("wcs_str_%04d" % i, bytes, chunk_size) for i in range(nchunks))
+        wcs_cols, wcs_dtypes = _make_chunked_cols(list(wcs_map.values()), "wcs")
+        cols.extend(wcs_cols)
+        dtypes.extend(wcs_dtypes)
 
         if pointing is not None:
             # Currently, there is only one pointing for all the chips, but write it out
@@ -142,6 +149,13 @@ def write_wcs_map(self, name, wcs_map, pointing):
         data = np.array(list(zip(*cols)), dtype=dtypes)
         self.write_table(name, data)
 
+    def write_bandpass(self, name, bandpass):
+        """Write a serialized bandpass object in its own extension.
+        """
+        cols, dtypes = _make_chunked_cols([bandpass], 'bandpass')
+        data = np.array(list(zip(*cols)), dtype=dtypes)
+        self.write_table(name, data)
+
     @contextmanager
     def nested(self, name):
         """Return a context manager that yields a new `FitsWriter` that nests
diff --git a/tests/test_input.py b/tests/test_input.py
index 968ba94e..93eb4894 100644
--- a/tests/test_input.py
+++ b/tests/test_input.py
@@ -16,6 +16,7 @@
 import os
 import shutil
 import galsim
+import galsim.roman
 import numpy as np
 import fitsio
 import piff
@@ -23,6 +24,7 @@
 import pyarrow
 import pyarrow.parquet
 
+from piff.util import make_flat
 from piff_test_helper import get_script_name, timer, CaptureLog
 
 @pytest.fixture(scope="module", autouse=True)
@@ -2117,7 +2119,7 @@ def test_sky():
 
 @timer
 def test_sed_col():
-    """Test per-star SED loading from a catalog filename column (`input.sed_col`).
+    """Test per-star SED loading from a catalog filename column (input.sed_col).
     """
     if __name__ == '__main__':
         logger = piff.config.setup_logger(verbose=2)
@@ -2134,25 +2136,31 @@ def test_sed_col():
         'sed_col': 'sed_file',
         'sed_wave_type': 'nm',
         'sed_flux_type': 'erg cm**(-2) s**(-1) angstrom**(-1)',
+        'bandpass': {'type': 'RomanBandpass', 'name': 'H158'},
     }
     input = piff.InputFiles(config, logger=logger)
     _, _, image_pos, extra_props = input.getRawImageData(0)
+    bandpass = galsim.roman.getBandpasses()['H158']
+    flat_bandpass = make_flat(bandpass)
     assert len(image_pos) == 100
-    assert len(extra_props['sed']) == 100
+    assert len(extra_props['sed_eff']) == 100
 
     stars = input.makeStars(logger=logger)
     assert len(stars) > 0
     for star in stars[:8]:
-        assert isinstance(star['sed'], galsim.SED)
-        assert np.isfinite(float(star['sed'](1000.0)))
+        assert isinstance(star['sed_eff'], galsim.SED)
+        np.testing.assert_allclose(
+            star['sed_eff'].calculateFlux(flat_bandpass), 1.0, rtol=1.0e-12, atol=0.0
+        )
+        assert star['sed_tol'] == 0.0
 
     # Repeated file names should reuse the same cached SED object.
-    assert extra_props['sed'][0] is extra_props['sed'][4]
-    assert extra_props['sed'][1] is extra_props['sed'][5]
-    assert extra_props['sed'][0] is not extra_props['sed'][1]
+    assert extra_props['sed_eff'][0] is extra_props['sed_eff'][4]
+    assert extra_props['sed_eff'][1] is extra_props['sed_eff'][5]
+    assert extra_props['sed_eff'][0] is not extra_props['sed_eff'][1]
 
     # Distinct template spectra should have different flux at a representative wavelength.
-    sed_vals = [float(extra_props['sed'][i](1000.0)) for i in range(4)]
+    sed_vals = [float(extra_props['sed_eff'][i](1500.0)) for i in range(4)]
     assert len(np.unique(np.round(sed_vals, 12))) > 1
 
     # Invalid sed_col should raise.
@@ -2163,18 +2171,18 @@ def test_sed_col():
     # With a sed_file_name dict, sed_col entries act as labels into the mapping.
     config3 = dict(config, sed_col='sed_label', sed_file_name='test_sed_map.yaml')
     _, _, _, extra_props2 = piff.InputFiles(config3, logger=logger).getRawImageData(0)
-    assert extra_props2['sed'][0] is extra_props2['sed'][4]
-    assert extra_props2['sed'][1] is extra_props2['sed'][5]
-    assert extra_props2['sed'][0] is not extra_props2['sed'][1]
+    assert extra_props2['sed_eff'][0] is extra_props2['sed_eff'][4]
+    assert extra_props2['sed_eff'][1] is extra_props2['sed_eff'][5]
+    assert extra_props2['sed_eff'][0] is not extra_props2['sed_eff'][1]
     np.testing.assert_allclose(
-        float(extra_props2['sed'][0](1000.0)),
-        float(extra_props['sed'][0](1000.0)),
+        float(extra_props2['sed_eff'][0](1500.0)),
+        float(extra_props['sed_eff'][0](1500.0)),
         rtol=1.0e-12,
         atol=0.0,
     )
     np.testing.assert_allclose(
-        float(extra_props2['sed'][1](1000.0)),
-        float(extra_props['sed'][1](1000.0)),
+        float(extra_props2['sed_eff'][1](1500.0)),
+        float(extra_props['sed_eff'][1](1500.0)),
         rtol=1.0e-12,
         atol=0.0,
     )
@@ -2198,7 +2206,7 @@ def test_sed_col():
 
 @timer
 def test_single_sed():
-    """Test single-file SED mode (`sed_file_name`)
+    """Test single-file SED mode (sed_file_name)
     """
     if __name__ == '__main__':
         logger = piff.config.setup_logger(verbose=2)
@@ -2214,20 +2222,37 @@ def test_single_sed():
         'sed_file_name': 'vega.txt',
         'sed_wave_type': 'Angstrom',
         'sed_flux_type': 'flambda',
+        'bandpass': {
+            'type': 'Eval',
+            'str': "galsim.Bandpass(lambda wave: 1.0, 'Angstrom', "
+                   "blue_limit=1000.0, red_limit=2000.0)",
+        },
     }
     input = piff.InputFiles(config, logger=logger)
     _, _, image_pos, extra_props = input.getRawImageData(0)
     assert len(image_pos) == 100
-    assert len(extra_props['sed']) == 100
+    assert len(extra_props['sed_eff']) == 100
 
+    bandpass = galsim.Bandpass(lambda wave: 1.0, 'Angstrom', blue_limit=1000, red_limit=2000)
     ref_sed = galsim.SED('vega.txt', wave_type='Angstrom', flux_type='flambda')
+    ref_sed = ref_sed.withFlux(1.0, bandpass)
     np.testing.assert_allclose(
-        float(extra_props['sed'][0](100.0)),
+        float(extra_props['sed_eff'][0](100.0)),
         float(ref_sed(100.0)),
         rtol=1.0e-12,
         atol=0.0,
     )
 
+    # An analytic SED expression combined with an analytic bandpass can still have no wave_list
+    # after forming sed * bandpass, which should take the len(wave) == 0 branch in _read_sed_file.
+    config['sed_file_name'] = '1./wave'
+    input = piff.InputFiles(config, logger=logger)
+    _, _, image_pos, extra_props = input.getRawImageData(0)
+    assert len(image_pos) == 100
+    assert len(extra_props['sed_eff']) == 100
+    assert extra_props['sed_eff'][0] is extra_props['sed_eff'][1]
+    assert len(extra_props['sed_eff'][0].wave_list) == 0
+
     # A single XSL SED file can be applied to all stars.
     config = {
         'dir': 'input',
@@ -2238,14 +2263,17 @@ def test_single_sed():
         'sed_file_name': 'xsl_spectrum_X0203_merged.fits',
         'sed_wave_type': 'nm',
         'sed_flux_type': 'erg cm**(-2) s**(-1) angstrom**(-1)',
+        'bandpass': {'type': 'RomanBandpass', 'name': 'H158'},
     }
+    piff.InputFiles._sed_cache = {}
     input = piff.InputFiles(config, logger=logger)
     _, _, _, props = input.getRawImageData(0)
-    assert len(props['sed']) == 100
-    assert props['sed'][0] is props['sed'][1]
+    assert len(props['sed_eff']) == 100
+    assert props['sed_eff'][0] is props['sed_eff'][1]
+    assert props['sed_tol'][0] == 0.0
     np.testing.assert_allclose(
-        float(props['sed'][0](1000.0)),
-        float(props['sed'][37](1000.0)),
+        float(props['sed_eff'][0](1500.0)),
+        float(props['sed_eff'][37](1500.0)),
         rtol=0.0,
         atol=0.0,
     )
@@ -2322,6 +2350,128 @@ def test_single_sed():
         input.getRawImageData(0)
     assert "sed_flux_key" in str(e.value)
 
+    config6 = dict(config)
+    del config6['bandpass']
+    with pytest.raises(ValueError) as e:
+        piff.InputFiles(config6, logger=logger)
+    assert "bandpass is required" in str(e.value)
+
+
+@timer
+def test_sed_thin():
+    """Test optional input-time SED reduction (sed_tol, sed_max_samples)."""
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=2)
+    else:
+        logger = None
+
+    config = {
+        'dir': 'input',
+        'image_file_name': 'test_input_image_00.fits',
+        'cat_file_name': 'test_input_cat_00.fits',
+        'x_col': 'x',
+        'y_col': 'y',
+        'sed_file_name': 'xsl_spectrum_X0203_merged.fits',
+        'sed_wave_type': 'nm',
+        'sed_flux_type': 'erg cm**(-2) s**(-1) angstrom**(-1)',
+        'bandpass': {'type': 'RomanBandpass', 'name': 'H158'},
+    }
+
+    bandpass = galsim.roman.getBandpasses()['H158']
+    flat_bandpass = make_flat(bandpass)
+
+    def weighted_line_fit(wave, flux):
+        weights = np.empty_like(wave)
+        weights[0] = 0.5 * (wave[1] - wave[0])
+        weights[-1] = 0.5 * (wave[-1] - wave[-2])
+        weights[1:-1] = 0.5 * (wave[2:] - wave[:-2])
+        X = np.column_stack([np.ones_like(wave), wave])
+        Xw = X * np.sqrt(weights)[:, np.newaxis]
+        yw = flux * np.sqrt(weights)
+        return np.linalg.lstsq(Xw, yw, rcond=None)[0]
+
+    # First default is no thinning.
+    piff.InputFiles._sed_cache = {}
+    _, _, _, props = piff.InputFiles(config, logger=logger).getRawImageData(0)
+    sed = props['sed_eff'][0]
+    assert props['sed_eff'][0] is props['sed_eff'][1]
+    assert props['sed_tol'][0] == 0.0
+    assert props['sed_max_samples'][0] == 0
+    np.testing.assert_allclose(
+        sed.calculateFlux(flat_bandpass),
+        1.0,
+        rtol=1.0e-12,
+        atol=0.0,
+    )
+
+    # Now test with a non-trivial thinning value.
+    piff.InputFiles._sed_cache = {}
+    thin_config = dict(config, sed_tol=0.1)
+    _, _, _, thin_props = piff.InputFiles(thin_config, logger=logger).getRawImageData(0)
+    thin_sed = thin_props['sed_eff'][0]
+    assert thin_props['sed_eff'][0] is thin_props['sed_eff'][1]
+    assert thin_props['sed_tol'][0] == 0.1
+    assert thin_props['sed_max_samples'][0] == 0
+    assert thin_sed is not sed
+    assert len(thin_sed.wave_list) < len(sed.wave_list)
+    np.testing.assert_allclose(
+        thin_sed.calculateFlux(flat_bandpass), 1.0, rtol=1.0e-12, atol=0.0
+    )
+    assert thin_sed.blue_limit == bandpass.blue_limit
+    assert thin_sed.red_limit == bandpass.red_limit
+
+    # Hard cap the number of samples and refit their values for the best piecewise-linear
+    # approximation of the original effective SED.
+    piff.InputFiles._sed_cache = {}
+    sample_config = dict(config, sed_max_samples=2)
+    _, _, _, sample_props = piff.InputFiles(sample_config, logger=logger).getRawImageData(0)
+    sample_sed = sample_props['sed_eff'][0]
+    assert sample_props['sed_eff'][0] is sample_props['sed_eff'][1]
+    assert sample_props['sed_tol'][0] == 0.0
+    assert sample_props['sed_max_samples'][0] == 2
+    assert len(sample_sed.wave_list) == 2
+    np.testing.assert_allclose(
+        sample_sed.calculateFlux(flat_bandpass), 1.0, rtol=1.0e-12, atol=0.0
+    )
+    wave = np.array(sed.wave_list, dtype=float)
+    flux = np.array([sed(w) for w in wave], dtype=float)
+    sample_wave = np.array(sample_sed.wave_list, dtype=float)
+    sample_flux = np.array([sample_sed(w) for w in sample_wave], dtype=float)
+    fit0, fit1 = weighted_line_fit(wave, flux)
+    sample0, sample1 = weighted_line_fit(sample_wave, sample_flux)
+    np.testing.assert_allclose(sample1, fit1, rtol=1.0e-3, atol=0.0)
+
+    # The two options can be used together: first use sed_tol to choose candidate samples,
+    # then trim to a maximum count and refit their values.
+    piff.InputFiles._sed_cache = {}
+    combo_config = dict(config, sed_tol=0.1, sed_max_samples=4)
+    _, _, _, combo_props = piff.InputFiles(combo_config, logger=logger).getRawImageData(0)
+    combo_sed = combo_props['sed_eff'][0]
+    assert combo_props['sed_tol'][0] == 0.1
+    assert combo_props['sed_max_samples'][0] == 4
+    assert len(combo_sed.wave_list) == 4
+    np.testing.assert_allclose(
+        combo_sed.calculateFlux(flat_bandpass), 1.0, rtol=1.0e-12, atol=0.0
+    )
+
+    # A large enough sed_max_samples should be a no-op after the sed_tol thinning stage.
+    piff.InputFiles._sed_cache = {}
+    no_op_config = dict(config, sed_tol=0.1, sed_max_samples=100)
+    _, _, _, no_op_props = piff.InputFiles(no_op_config, logger=logger).getRawImageData(0)
+    no_op_sed = no_op_props['sed_eff'][0]
+    assert no_op_props['sed_tol'][0] == 0.1
+    assert no_op_props['sed_max_samples'][0] == 100
+    assert len(no_op_sed.wave_list) == len(thin_sed.wave_list)
+    np.testing.assert_array_equal(
+        [no_op_sed(w) for w in no_op_sed.wave_list],
+        [thin_sed(w) for w in thin_sed.wave_list],
+    )
+
+    # sed_max_samples must be >= 2 if provided.
+    piff.InputFiles._sed_cache = {}
+    bad_config = dict(config, sed_max_samples=1)
+    with np.testing.assert_raises(ValueError):
+        piff.InputFiles(bad_config, logger=logger).getRawImageData(0)
 
 @timer
 def test_sed_star_io():
@@ -2341,15 +2491,21 @@ def test_sed_star_io():
         'sed_col': 'sed_file',
         'sed_wave_type': 'nm',
         'sed_flux_type': 'erg cm**(-2) s**(-1) angstrom**(-1)',
+        'sed_tol': 0.1,
+        'sed_max_samples': 6,
+        'bandpass': {'type': 'RomanBandpass', 'name': 'H158'},
     }
+    piff.InputFiles._sed_cache = {}
     stars = piff.InputFiles(config, logger=logger).makeStars(logger=logger)
     assert len(stars) > 10
-    assert 'sed' in stars[0].data.properties
+    assert 'sed_eff' in stars[0].data.properties
     assert 'sed_file_name' in stars[0].data.properties
     assert 'sed_wave_type' in stars[0].data.properties
     assert 'sed_flux_type' in stars[0].data.properties
     assert 'sed_wave_key' in stars[0].data.properties
     assert 'sed_flux_key' in stars[0].data.properties
+    assert 'sed_tol' in stars[0].data.properties
+    assert 'sed_max_samples' in stars[0].data.properties
 
     file_name = os.path.join('output', 'sed_star_io.fits')
     with piff.writers.FitsWriter.open(file_name) as w:
@@ -2363,26 +2519,131 @@ def test_sed_star_io():
     assert 'sed_flux_type' in cols
     assert 'sed_wave_key' in cols
     assert 'sed_flux_key' in cols
+    assert 'sed_tol' in cols
+    assert 'sed_max_samples' in cols
 
+    bandpass = galsim.roman.getBandpasses()['H158']
     with piff.readers.FitsReader.open(file_name) as r:
-        stars2 = piff.Star.read(r, 'stars')
-
+        with pytest.raises(ValueError) as e:
+            piff.Star.read(r, 'stars')
+    assert "bandpass is required" in str(e.value)
+    with piff.readers.FitsReader.open(file_name) as r:
+        stars2 = piff.Star.read(r, 'stars', bandpass=bandpass)
     for s1, s2 in zip(stars, stars2):
-        assert isinstance(s2['sed'], galsim.SED)
+        assert isinstance(s2['sed_eff'], galsim.SED)
         assert s1['sed_file_name'] == s2['sed_file_name']
         assert s1['sed_wave_type'] == s2['sed_wave_type']
         assert s1['sed_flux_type'] == s2['sed_flux_type']
         assert s1['sed_wave_key'] == s2['sed_wave_key']
         assert s1['sed_flux_key'] == s2['sed_flux_key']
+        assert s1['sed_tol'] == s2['sed_tol'] == 0.1
+        assert s1['sed_max_samples'] == s2['sed_max_samples'] == 6
         np.testing.assert_allclose(
-            float(s2['sed'](1000.0)),
-            float(s1['sed'](1000.0)),
+            float(s2['sed_eff'](1500.0)),
+            float(s1['sed_eff'](1500.0)),
             rtol=1.0e-12,
             atol=0.0,
         )
     # Same file names should still share one cached SED after deserialization.
-    assert stars2[0]['sed'] is stars2[4]['sed']
-    assert stars2[1]['sed'] is stars2[5]['sed']
+    assert stars2[0]['sed_eff'] is stars2[4]['sed_eff']
+    assert stars2[1]['sed_eff'] is stars2[5]['sed_eff']
+
+
+@timer
+def test_wcs_bandpass_io():
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=2)
+    else:
+        logger = None
+
+    file_name = os.path.join('output', 'wcs_bandpass_io.fits')
+    wcs = {1: galsim.PixelScale(0.2)}
+    pointing = galsim.CelestialCoord(6 * galsim.hours, -30 * galsim.degrees)
+    bandpass = galsim.roman.getBandpasses()['H158']
+
+    with piff.writers.FitsWriter.open(file_name) as w:
+        w.write_wcs_map('wcs', wcs, pointing)
+        w.write_bandpass('bandpass', bandpass)
+
+    with piff.readers.FitsReader.open(file_name) as r:
+        wcs2, pointing2 = r.read_wcs_map('wcs', logger=logger)
+        bandpass2 = r.read_bandpass('bandpass')
+
+    assert list(wcs2.keys()) == [1]
+    assert repr(wcs2[1]) == repr(wcs[1])
+    assert pointing2 == pointing
+    assert bandpass2.blue_limit == bandpass.blue_limit
+    assert bandpass2.red_limit == bandpass.red_limit
+    np.testing.assert_allclose(bandpass2(1500.0), bandpass(1500.0), rtol=0.0, atol=0.0)
+
+
+@timer
+def test_input_bandpass():
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=2)
+    else:
+        logger = None
+
+    config = {
+        'input': {
+            'dir': 'input',
+            'image_file_name': 'test_input_image_00.fits',
+            'cat_file_name': 'test_input_cat_00.fits',
+            'flag_col': 'flag',
+            'use_flag': 1,
+            'skip_flag': 4,
+            'stamp_size': 15,
+            'bandpass': {'type': 'RomanBandpass', 'name': 'H158'},
+        },
+        'psf': {
+            'model': {'type': 'Gaussian', 'fastfit': True, 'include_pixel': False},
+            'interp': {'type': 'Mean'},
+            'min_iter': 1,
+            'max_iter': 1,
+        },
+    }
+
+    objects, wcs, pointing, bandpass = piff.Input.process(config['input'], logger=logger)
+    assert len(objects) > 0
+    assert wcs is not None
+    assert pointing is None
+    assert bandpass is not None
+
+    # PSF.process creates a PSF object without the wcs, pointing, bandpass context yet.
+    psf = piff.PSF.process(config['psf'], logger=logger)
+    assert psf.wcs is None
+    assert psf.pointing is None
+    assert psf.bandpass is None
+
+    # set_context assigns them as attributes of the psf.
+    psf.set_context(wcs, pointing, bandpass)
+    assert psf.wcs == wcs
+    assert psf.pointing == pointing
+    assert psf.bandpass == bandpass
+
+    psf2 = piff.process(config, logger=logger)
+    assert psf2.wcs == wcs
+    assert psf2.pointing == pointing
+    assert psf2.bandpass == bandpass
+
+    # The old pattern was to provide wcs and pointing to fit rather than process.
+    # It is still supported, but should emit a deprecation warning and be equivalent.
+    with CaptureLog() as cl:
+        psf_old = piff.PSF.process(config['psf'], logger=cl.logger)
+        psf_old.fit(objects, wcs, pointing, logger=cl.logger)
+    assert "wcs should now be set with set_context" in cl.output
+    assert psf_old.bandpass is None  # Old pattern didn't include bandpass
+    test_im_old = psf_old.draw(x=123.4, y=234.5, chipnum=0)
+    test_im_new = psf2.draw(x=123.4, y=234.5, chipnum=0)
+    np.testing.assert_allclose(test_im_old.array, test_im_new.array, atol=1.0e-12)
+
+    # Check roundtrip through file
+    file_name = os.path.join('output', 'input_bandpass.piff')
+    psf2.write(file_name, logger=logger)
+    psf3 = piff.read(file_name, logger=logger)
+    assert psf3.wcs == wcs
+    assert psf3.pointing == pointing
+    assert psf3.bandpass == bandpass
 
 
 if __name__ == '__main__':
@@ -2399,3 +2660,5 @@ def test_sed_star_io():
     test_sed_col()
     test_single_sed()
     test_sed_star_io()
+    test_wcs_bandpass_io()
+    test_input_bandpass()
diff --git a/tests/test_pixel.py b/tests/test_pixel.py
index 816a2868..67f53eac 100644
--- a/tests/test_pixel.py
+++ b/tests/test_pixel.py
@@ -985,9 +985,9 @@ def test_single_image():
         # Process the star data
         model = piff.PixelGrid(0.2, 16)
         interp = piff.BasisPolynomial(order=order)
-        pointing = None     # wcs is not Celestial here, so pointing needs to be None.
         psf = piff.SimplePSF(model, interp)
-        psf.fit(orig_stars, {0:image1.wcs}, pointing, logger=logger)
+        psf.set_context(wcs={0:image1.wcs})
+        psf.fit(orig_stars, logger=logger)
 
         # Check that the interpolation is what it should be
         print('target.flux = ',target_star.fit.flux)
@@ -1186,7 +1186,7 @@ def test_des_image():
 
         # Largely copied from Gary's fit_des.py, but using the Piff input_handler to
         # read the input files.
-        stars, wcs, pointing = piff.Input.process(config['input'], logger=logger)
+        stars, wcs, pointing, _ = piff.Input.process(config['input'], logger=logger)
         if nstars is not None:
             stars = stars[:nstars]
 
@@ -1200,7 +1200,8 @@ def test_des_image():
 
         # Make a psf
         psf = piff.SimplePSF(model, interp)
-        psf.fit(stars, wcs, pointing, logger=logger)
+        psf.set_context(wcs, pointing)
+        psf.fit(stars, logger=logger)
 
         # The difference between the images of the fitted stars and the originals should be
         # consistent with noise.  Keep track of how many don't meet that goal.
@@ -1257,7 +1258,7 @@ def test_des_image():
     print('start piffify')
     piff.piffify(config)
     print('read stars')
-    stars, wcs, pointing = piff.Input.process(config['input'])
+    stars, wcs, pointing, _ = piff.Input.process(config['input'])
     print('read psf')
     psf = piff.read(psf_file)
     nreserve = len([s for s in psf.stars if s.is_reserve])
@@ -1294,7 +1295,7 @@ def test_des_image():
         p = subprocess.Popen( [piffify_exe, 'pixel_des.yaml'] )
         p.communicate()
         print('read stars')
-        stars, wcs, pointing = piff.Input.process(config['input'])
+        stars, wcs, pointing, _ = piff.Input.process(config['input'])
         print('read psf')
         psf = piff.read(psf_file)
         stars = [psf.model.initialize(s) for s in stars]
@@ -1417,7 +1418,7 @@ def test_des2():
     piff.piffify(config)
     t1 = time.time()
     print('Time for direct ATA solution = ',t1-t0)
-    stars, wcs, pointing = piff.Input.process(config['input'])
+    stars, wcs, pointing, _ = piff.Input.process(config['input'])
     psf = piff.read(psf_file)
     stars = psf.initialize_flux_center(stars)
 
@@ -1897,9 +1898,9 @@ def test_convergence_centering_failed():
         )
         stars.append(piff.Star(data, None))
 
-
     piffResult = piff.PSF.process(piffConfig)
-    piffResult.fit(stars, wcs, pointing)
+    piffResult.set_context(wcs)
+    piffResult.fit(stars)
 
     assert piffResult.niter == 6, 'Maximum number of iterations in this example must be 6.'
     assert np.shape(piffResult.stars[28].fit.A) == (0, 625), 'Centroid failed for this star, expected shape of A is (0,625)'
@@ -1935,23 +1936,24 @@ def test_too_small():
         }
     }
 
-    piffResult = piff.PSF.process(piffConfig)
     # Run on a single CCD, and in image coords rather than sky coords.
     wcs = {0: galsim.PixelScale(1.0)}
-    pointing = None
+    piffResult = piff.PSF.process(piffConfig)
+    piffResult.set_context(wcs)
 
     # The original behavior was to complete successfully, but the model had crazy values
     # off the edge of the constrained region.
     with CaptureLog(2) as cl:
         with np.testing.assert_raises(RuntimeError):
-            piffResult.fit(stars, wcs, pointing, logger=cl.logger)
+            piffResult.fit(stars, logger=cl.logger)
     assert "Image size (21,21) is too small to constrain this PixelGrid." in cl.output
     assert "Removed 12 stars in initialize" in cl.output
 
     # This configuration works.
     piffConfig['model']['size'] = 21
     piffResult = piff.PSF.process(piffConfig)
-    piffResult.fit(stars, wcs, pointing)
+    piffResult.set_context(wcs)
+    piffResult.fit(stars)
     im = piffResult.draw(10, 10)
     assert im(10,10) > 0
 
diff --git a/tests/test_roman.py b/tests/test_roman.py
index e45651bb..c0c9bff9 100644
--- a/tests/test_roman.py
+++ b/tests/test_roman.py
@@ -23,11 +23,19 @@
 import piff
 import galsim
 import pytest
+import fitsio
 
+from piff.util import make_flat
 from piff.roman import RomanOpticsPSF, RomanOpticalModel, RomanSCAInterp
 from piff_test_helper import timer
 
 
+def make_roman_model(**kwargs):
+    model = RomanOpticalModel(**kwargs)
+    model.set_bandpass(galsim.roman.getBandpasses()['H158'])
+    return model
+
+
 @contextmanager
 def fast_pupil_bin(value=16):
     """Temporarily raise Roman `pupil_bin` to speed up the getPSF calls in unit tests.
@@ -48,14 +56,16 @@ def test_roman_optics():
     """Check RomanOptics basic construction, drawing, and SCA/chipnum property handling.
     """
     with fast_pupil_bin():
+        bandpass = galsim.roman.getBandpasses()['H158']
+
         # Check basic construction.
         psf = piff.PSF.process(
-            {'type': 'RomanOptics', 'filter': 'H158', 'chromatic': False, 'max_zernike': 6}
+            {'type': 'RomanOptics', 'chromatic': False, 'max_zernike': 6}
         )
         assert isinstance(psf, RomanOpticsPSF)
         logger = piff.config.setup_logger()
         assert psf.interp_property_names == ('sca',)
-        assert psf.fit_center is False
+        assert psf.fit_center is True
         assert psf.include_model_centroid is False
         assert psf.model.aberration_interp == 'constant'
         assert psf.model.nominal_interp == 'bilinear'
@@ -65,7 +75,6 @@ def test_roman_optics():
         psf_global = piff.PSF.process(
             {
                 'type': 'RomanOptics',
-                'filter': 'H158',
                 'chromatic': False,
                 'max_zernike': 6,
                 'aberration_interp': 'global',
@@ -85,8 +94,7 @@ def test_roman_optics():
 
         stars, nremoved = psf.initialize_params([star], logger=logger)
         assert nremoved == 0
-        psf.wcs = {5: galsim.PixelScale(0.11)}
-        psf.pointing = None
+        psf.set_context({5: galsim.PixelScale(0.11)}, None, bandpass)
         psf.interp.set_sca_solution({5: np.zeros(psf.model.param_len)})
         prof, method = psf.get_profile(x=123.4, y=456.7, sca=5)
         assert prof is not None
@@ -95,6 +103,42 @@ def test_roman_optics():
         assert image.array.shape == (48, 48)  # Default stamp_size=48 in draw function.
         assert np.isclose(image.array.sum(), 1.0, rtol=0.05)
 
+        chromatic_psf = piff.PSF.process(
+            {
+                'type': 'RomanOptics',
+                'chromatic': True,
+                'max_zernike': 6,
+            }
+        )
+        chromatic_psf.set_context({5: galsim.PixelScale(0.11)}, None, bandpass)
+        sed = galsim.SED(lambda w: 1.0, wave_type='nm', flux_type='fphotons')
+        chromatic_star = piff.Star.makeTarget(
+            x=123.4,
+            y=456.7,
+            stamp_size=25,
+            scale=0.11,
+            properties={'sca': 5, 'sed_eff': sed},
+        ).withFlux(1.0, (0.0, 0.0))
+        chromatic_stars, _ = chromatic_psf.initialize_params([chromatic_star], logger=logger)
+        chromatic_psf.interp.set_sca_solution({5: np.zeros(chromatic_psf.model.param_len)})
+        chromatic_model_star = chromatic_psf.drawStar(chromatic_stars[0])
+        assert np.isfinite(chromatic_model_star.image.array).all()
+        assert chromatic_model_star.image.array.sum() > 0.0
+        chromatic_psf.drawStar(chromatic_stars[0])
+        assert chromatic_stars[0].data.properties['sed_eff'] is sed
+
+        missing_sed_star = piff.Star.makeTarget(
+            x=123.4,
+            y=456.7,
+            stamp_size=25,
+            scale=0.11,
+            properties={'sca': 5},
+        ).withFlux(1.0, (0.0, 0.0))
+        missing_sed_stars, _ = chromatic_psf.initialize_params([missing_sed_star], logger=logger)
+        with pytest.raises(ValueError) as err:
+            chromatic_psf.drawStar(missing_sed_stars[0])
+        assert "requires each star to have an 'sed_eff' property" in str(err.value)
+
         # Helper function to use as a side-effect of getPSF to record what sca argument is used.
         used_sca = [] # Use a list so we can modify in place.
         real_get_psf = galsim.roman.getPSF  # Save this here to avoid recursion when patched.
@@ -159,7 +203,6 @@ def traced_get_psf(*args, **kwargs):
         psf1 = piff.PSF.process(
             {
                 'type': 'RomanOptics',
-                'filter': 'H158',
                 'chromatic': False,
                 'max_zernike': 6,
                 'outliers': [
@@ -168,6 +211,7 @@ def traced_get_psf(*args, **kwargs):
                 ],
             }
         )
+        psf1.set_context(None, None, bandpass)
         assert psf1.outliers is not None
         assert isinstance(psf1.outliers, list)
         fn = os.path.join('output', 'roman_outliers_write_test.piff')
@@ -178,27 +222,71 @@ def traced_get_psf(*args, **kwargs):
         assert psf2.outliers[0]._type_name == 'Chisq'
         assert psf2.outliers[1]._type_name == 'Centroid'
 
+        # Round trip WCS + bandpass through the read path where set_context is called
+        # before _finish_read has replaced the placeholder model.
+        psf1 = piff.PSF.process(
+            {
+                'type': 'RomanOptics',
+                'chromatic': True,
+                'max_zernike': 6,
+            }
+        )
+        wcs = {5: galsim.PixelScale(0.11)}
+        psf1.set_context(wcs, None, bandpass)
+        fn = os.path.join('output', 'roman_context_write_test.piff')
+        psf1.write(fn)
+        psf2 = piff.read(fn)
+        assert psf2.wcs == wcs
+        assert psf2.bandpass == bandpass
+        assert psf2.model.bandpass == bandpass
+        assert psf2.model.filter_name == 'H158'
+
     # max_zernike must be >= 4
     with pytest.raises(ValueError) as err:
-        RomanOpticalModel(filter='H158', chromatic=False, max_zernike=3)
+        RomanOpticalModel(chromatic=False, max_zernike=3)
     assert "range 4..22" in str(err.value)
 
     # max_zernike must be <= 22
     with pytest.raises(ValueError) as err:
-        RomanOpticalModel(filter='H158', chromatic=False, max_zernike=23)
+        RomanOpticalModel(chromatic=False, max_zernike=23)
     assert "range 4..22" in str(err.value)
 
-    # Check invalid filter string
+    # A custom bandpass can be used as long as the Roman filter_name is given explicitly.
+    model = RomanOpticalModel(chromatic=False, max_zernike=6)
+    custom_bandpass = galsim.Bandpass(lambda wave: 1.0, 'nm', blue_limit=1000, red_limit=2000)
+    with pytest.raises(ValueError) as err:
+        model.getProfile(params=np.zeros(model.param_len), star=star)
+    assert "requires bandpass to be set before use" in str(err.value)
+    with pytest.raises(ValueError) as err:
+        model.set_bandpass(custom_bandpass)
+    assert "requires filter_name" in str(err.value)
+    model.set_bandpass(custom_bandpass, filter_name='H158')
+    assert model.bandpass == custom_bandpass
+    assert model.filter_name == 'H158'
+
+    # It's ok, but unnecessary to give the filter_name explicitly with a roman bandpass.
+    model.set_bandpass(galsim.roman.getBandpasses()['H158'], filter_name='H158')
+    assert model.filter_name == 'H158'      # The explicitly set filter_name.
+    assert model.bandpass.name == 'H158'    # The one that galsim.roman sets in the bandpass obj.
+
+    # Invalid filter_name fails when GalSim is actually asked to build the Roman PSF.
+    model.filter_name = None  # Lets set_bandpass actually work.
+    model.set_bandpass(
+        galsim.Bandpass(lambda wave: 1.0, 'nm', blue_limit=1000, red_limit=2000),
+        filter_name='NotAFilter',
+    )
     with pytest.raises(ValueError) as err:
-        RomanOpticalModel(filter='NotAFilter', chromatic=False, max_zernike=6)
-    assert "not a valid GalSim Roman bandpass" in str(err.value)
+        model.getProfile(params=np.zeros(model.param_len), star=star)
+    assert "valid Roman bandpass" in str(err.value)
 
+    # Invalid aberration_interp
     with pytest.raises(ValueError) as err:
-        RomanOpticalModel(filter='H158', chromatic=False, max_zernike=6, aberration_interp='bad')
+        RomanOpticalModel(chromatic=False, max_zernike=6, aberration_interp='bad')
     assert "must be one of" in str(err.value)
 
+    # Invalid nominal_interp
     with pytest.raises(ValueError) as err:
-        RomanOpticalModel(filter='H158', chromatic=False, max_zernike=6, nominal_interp='bad')
+        RomanOpticalModel(chromatic=False, max_zernike=6, nominal_interp='bad')
     assert "must be one of" in str(err.value)
 
 @timer
@@ -206,9 +294,11 @@ def test_corner_cache():
     """Verify corner-profile caching reuses one 4-corner set for same SCA and params.
     """
     with fast_pupil_bin():
+        bandpass = galsim.roman.getBandpasses()['H158']
         psf = piff.PSF.process(
-            {'type': 'RomanOptics', 'filter': 'H158', 'chromatic': False, 'max_zernike': 6}
+            {'type': 'RomanOptics', 'chromatic': False, 'max_zernike': 6}
         )
+        psf.set_context(None, None, bandpass)
         logger = piff.config.setup_logger()
         stars = [
             piff.Star.makeTarget(
@@ -242,12 +332,12 @@ def test_corner_cache():
         psf5 = piff.PSF.process(
             {
                 'type': 'RomanOptics',
-                'filter': 'H158',
                 'chromatic': False,
                 'max_zernike': 6,
                 'nominal_interp': 'five_point',
             }
         )
+        psf5.set_context(None, None, bandpass)
         stars5, _ = psf5.initialize_params(stars, logger=logger)
         psf5.drawStar(stars5[0])
         params5 = stars5[0].fit.params
@@ -260,12 +350,7 @@ def test_corner_cache():
 def test_five_point_weights():
     """Check five-point interpolation weights for corner/center and quadratic basis terms.
     """
-    model = RomanOpticalModel(
-        filter='H158',
-        chromatic=False,
-        max_zernike=6,
-        nominal_interp='five_point',
-    )
+    model = make_roman_model(chromatic=False, max_zernike=6, nominal_interp='five_point')
     sca = 7
     size = model.sca_size
 
@@ -318,12 +403,7 @@ def test_fit():
     """Check local convergence of single-star Roman fits in constant mode.
     """
     with fast_pupil_bin():
-        model = RomanOpticalModel(
-            filter='H158',
-            chromatic=False,
-            max_zernike=6,
-            aberration_prior_sigma=1.0e6,
-        )
+        model = make_roman_model(chromatic=False, max_zernike=6, aberration_prior_sigma=1.0e6)
         star = piff.Star.makeTarget(
             x=64.0,
             y=64.0,
@@ -338,7 +418,7 @@ def test_fit():
         # realistic while still providing measurable signal in this unit test.
         truth_params = np.array([0.004, -0.003, 0.005])
         prof = model.getProfile(truth_params, star=star)
-        model._draw_profile_to_image(prof * star.fit.flux, star.image, star.image_pos)
+        model._draw_profile_to_image(prof * star.fit.flux, star.image, star.image_pos, star)
         fitted = model.fit(star)
         print('fit = ',fitted.fit.params)
         # 1 pass isn't great, but after 2 passes, the agreement is sub percent.
@@ -384,18 +464,10 @@ def test_aberration_prior():
     """Validate the use of priors on the aberration values.
     """
     # Check that fitted aberrations stay closer to 0 when strong prior is applied.
-    weak_prior = RomanOpticalModel(
-        filter='H158',
-        chromatic=False,
-        max_zernike=6,
-        aberration_prior_sigma=1.0e6,
-    )
-    strong_prior = RomanOpticalModel(
-        filter='H158',
-        chromatic=False,
-        max_zernike=6,
-        aberration_prior_sigma=[0.02],  # List with 1 element treated as scalar.
-    )
+    weak_prior = make_roman_model(chromatic=False, max_zernike=6, aberration_prior_sigma=1.0e6)
+    strong_prior = make_roman_model(
+        chromatic=False, max_zernike=6, aberration_prior_sigma=[0.02]
+    )  # List with 1 element treated as scalar.
     aw = np.zeros((10, weak_prior.param_len))
     aw[:, 0] = 12.3
     aw[:, 1] = 1.7
@@ -407,12 +479,7 @@ def test_aberration_prior():
     assert abs(new_strong[0]) < abs(new_weak[0])
 
     # None means no prior, which is equivalent to infinite prior.
-    no_prior = RomanOpticalModel(
-        filter='H158',
-        chromatic=False,
-        max_zernike=6,
-        aberration_prior_sigma=None,
-    )
+    no_prior = make_roman_model(chromatic=False, max_zernike=6, aberration_prior_sigma=None)
     ata = np.eye(no_prior.param_len)
     atb = np.ones(no_prior.param_len)
     p = np.zeros(no_prior.param_len)
@@ -424,12 +491,7 @@ def test_aberration_prior():
     np.testing.assert_allclose(atb, atb0)
 
     # Compare the solve to what you get with an actually infinite prior.
-    inf_prior = RomanOpticalModel(
-        filter='H158',
-        chromatic=False,
-        max_zernike=6,
-        aberration_prior_sigma=np.inf
-    )
+    inf_prior = make_roman_model(chromatic=False, max_zernike=6, aberration_prior_sigma=np.inf)
     no_prior_result, _ = no_prior._solve_params(aw, bw, p0)
     inf_prior_result, _ = inf_prior._solve_params(aw, bw, p0)
     print('no prior',no_prior_result)
@@ -438,26 +500,26 @@ def test_aberration_prior():
 
     # prior length must match number of zernikes used
     with pytest.raises(ValueError) as err:
-        RomanOpticalModel(filter='H158', chromatic=False, max_zernike=6,
-                          aberration_prior_sigma=[1.0, 2.0])
+        make_roman_model(chromatic=False, max_zernike=6, aberration_prior_sigma=[1.0, 2.0])
     assert "scalar or length 3" in str(err.value)
 
     with pytest.raises(ValueError) as err:
-        RomanOpticalModel(filter='H158', chromatic=False, max_zernike=6,
-                          aberration_prior_sigma=[])
+        make_roman_model(chromatic=False, max_zernike=6, aberration_prior_sigma=[])
     assert "scalar or length 3" in str(err.value)
 
     # Same for linear mode (in particular the allowed length is 3 here, not 9,
     # which is the full param_len).
     with pytest.raises(ValueError) as err:
-        piff.roman.RomanOpticalModel(filter='H158', chromatic=False, max_zernike=6,
-                                     aberration_interp='linear',
-                                     aberration_prior_sigma=[1.0] * 9)
+        make_roman_model(
+            chromatic=False,
+            max_zernike=6,
+            aberration_interp='linear',
+            aberration_prior_sigma=[1.0] * 12,
+        )
     assert "scalar or length 3" in str(err.value)
 
     # In linear mode, scalar and one-point vectors tile to the (a, b, c) coefficient blocks.
-    linear_scalar = piff.roman.RomanOpticalModel(
-        filter='H158',
+    linear_scalar = make_roman_model(
         chromatic=False,
         max_zernike=6,
         aberration_interp='linear',
@@ -465,8 +527,7 @@ def test_aberration_prior():
     )
     np.testing.assert_allclose(linear_scalar.prior_sigma, np.full(9, 0.05))
 
-    linear_vec = piff.roman.RomanOpticalModel(
-        filter='H158',
+    linear_vec = make_roman_model(
         chromatic=False,
         max_zernike=6,
         aberration_interp='linear',
@@ -476,13 +537,11 @@ def test_aberration_prior():
 
     # priors cannot be <= 0
     with pytest.raises(ValueError) as err:
-        RomanOpticalModel(filter='H158', chromatic=False, max_zernike=6,
-                          aberration_prior_sigma=[1.0, 0.0, 1.0])
+        make_roman_model(chromatic=False, max_zernike=6, aberration_prior_sigma=[1.0, 0.0, 1.0])
     assert "must all be > 0" in str(err.value)
 
     with pytest.raises(ValueError) as err:
-        RomanOpticalModel(filter='H158', chromatic=False, max_zernike=6,
-                          aberration_prior_sigma=[1.0, -1.0, 1.0])
+        make_roman_model(chromatic=False, max_zernike=6, aberration_prior_sigma=[1.0, -1.0, 1.0])
     assert "must all be > 0" in str(err.value)
 
 
@@ -493,7 +552,6 @@ def test_linear_prior_io():
     psf = piff.PSF.process(
         {
             'type': 'RomanOptics',
-            'filter': 'H158',
             'chromatic': False,
             'max_zernike': 6,
             'aberration_interp': 'linear',
@@ -525,7 +583,6 @@ def test_linear_prior_io():
     psf_scalar = piff.PSF.process(
         {
             'type': 'RomanOptics',
-            'filter': 'H158',
             'chromatic': False,
             'max_zernike': 6,
             'aberration_interp': 'linear',
@@ -557,12 +614,7 @@ def test_fit_many():
     """Check accuracy of fitting multiple stars using fit_many.
     """
     with fast_pupil_bin():
-        model = RomanOpticalModel(
-            filter='H158',
-            chromatic=False,
-            max_zernike=6,
-            aberration_prior_sigma=1.0e6,
-        )
+        model = make_roman_model(chromatic=False, max_zernike=6, aberration_prior_sigma=1.0e6)
         stars = [
             piff.Star.makeTarget(
                 x=64.2,
@@ -594,7 +646,7 @@ def test_fit_many():
         ]
         for star in stars:
             prof = model.getProfile(truth_params, star=star)
-            model._draw_profile_to_image(prof * star.fit.flux, star.image, star.image_pos)
+            model._draw_profile_to_image(prof * star.fit.flux, star.image, star.image_pos, star)
 
         # 1 pass isn't great, but after 2 passes, the agreement is sub percent.
         # And 3 is within 0.1% agreement.
@@ -616,12 +668,8 @@ def test_fit_many_nproc():
     """Check `fit_many` multiprocessing path and accuracy with `nproc > 1`.
     """
     with fast_pupil_bin():
-        model = piff.roman.RomanOpticalModel(
-            filter='H158',
-            chromatic=False,
-            max_zernike=6,
-            aberration_prior_sigma=1.0e6,
-            nproc=2,
+        model = make_roman_model(
+            chromatic=False, max_zernike=6, aberration_prior_sigma=1.0e6, nproc=2
         )
         stars = [
             piff.Star.makeTarget(
@@ -653,7 +701,7 @@ def test_fit_many_nproc():
         ]
         for star in stars:
             prof = model.getProfile(truth_params, star=star)
-            model._draw_profile_to_image(prof * star.fit.flux, star.image, star.image_pos)
+            model._draw_profile_to_image(prof * star.fit.flux, star.image, star.image_pos, star)
 
         for _ in range(3):
             stars = model.fit_many(stars)
@@ -664,6 +712,110 @@ def test_fit_many_nproc():
         assert list(model._corner_cache.keys()) == [5]
 
 
+@timer
+def test_chromatic():
+    """Check realistic chromatic RomanOptics fitting with different stellar SEDs.
+    """
+    with fast_pupil_bin():
+        bandpass = galsim.roman.getBandpasses()['H158']
+        flat_bandpass = make_flat(bandpass)
+        psf = piff.PSF.process(
+            {
+                'type': 'RomanOptics',
+                'chromatic': True,
+                'max_zernike': 6,
+                'max_iter': 1,  # Not used until later when we run fit() directly.
+            }
+        )
+        psf.set_context(galsim.PixelScale(0.11), None, bandpass)
+        logger = piff.config.setup_logger()
+
+        def read_xsl_sed(file_name):
+            from astropy import units as u
+            data = fitsio.read(os.path.join(os.path.dirname(__file__), 'input', file_name))
+            table = galsim.LookupTable(data['WAVE'], data['FLUX'], interpolant='linear')
+            return galsim.SED(
+                table,
+                wave_type=u.nm,
+                flux_type=u.Unit('erg cm**(-2) s**(-1) angstrom**(-1)'),
+            )
+
+        # The Piff input process will make effective sed's (sed*bp) and thin them.
+        # Do that here for the purpose of this test.
+        sed_blue = read_xsl_sed('xsl_spectrum_X0529_merged.fits')
+        sed_red = read_xsl_sed('xsl_spectrum_X0066_merged_scl.fits')
+        sed_eff_blue = (sed_blue * bandpass).thin(rel_err=0.1).withFlux(1.0, flat_bandpass)
+        sed_eff_red = (sed_red * bandpass).thin(rel_err=0.1).withFlux(1.0, flat_bandpass)
+        print('sed_eff_blue wave_list = ',sed_eff_blue.wave_list)
+        print('sed_eff_red wave_list = ',sed_eff_red.wave_list)
+        assert len(sed_eff_blue.wave_list) <= 10
+        assert len(sed_eff_red.wave_list) <= 10
+
+        flux = 1.e5
+        truth_params = np.array([0.004, -0.003, 0.005])
+
+        # Make the stars with images using a less-thinned sed and bandpass.
+        stars = [
+            piff.Star.makeTarget(x=64.2, y=64.1, stamp_size=25, scale=0.11,
+                                 properties={'sca': 5, 'sed_eff': sed_eff_blue}).withFlux(flux),
+            piff.Star.makeTarget(x=171.8, y=162.7, stamp_size=25, scale=0.11,
+                                 properties={'sca': 5, 'sed_eff': sed_eff_red}).withFlux(flux),
+        ]
+        bp = bandpass.thin(rel_err=0.05)
+        for star, sed in zip(stars, [sed_blue, sed_red]):
+            prof = galsim.roman.getPSF(
+                5, 'H158',
+                SCA_pos=star.image_pos,
+                pupil_bin=piff.roman.roman_psf.pupil_bin,
+                wcs=star.image.wcs,
+                extra_aberrations=psf.model._make_extra_aberrations(truth_params),
+                wavelength=None,
+            )
+            prof = prof * sed.thin(rel_err=0.05).withFlux(flux, bp)
+            prof.drawImage(bp, image=star.image, center=star.image_pos)
+
+        # Do the fit
+        stars, nremoved = psf.initialize_params(stars, logger=logger)
+        assert nremoved == 0
+
+        for _ in range(2):
+            stars, nremoved = psf.single_iteration(
+                stars,
+                logger=logger,
+                convert_funcs=None,
+                draw_method=None,
+            )
+            assert nremoved == 0
+
+        p0 = stars[0].fit.params
+        p1 = stars[1].fit.params
+        print('chromatic p0 = ', p0)
+        print('chromatic p1 = ', p1)
+        print('       truth = ', truth_params)
+        np.testing.assert_allclose(p0, p1, atol=2.0e-6, rtol=0.0)
+        np.testing.assert_allclose(psf.interp.sca_mean[5], p0, atol=2.0e-6, rtol=0.0)
+
+        for star in stars:
+            model_star = psf.drawStar(star)
+            diff = model_star.image.array - star.image.array
+            frac_l2 = np.sqrt(np.sum(diff * diff) / np.sum(star.image.array ** 2))
+            frac_peak = np.max(np.abs(diff)) / np.max(np.abs(star.image.array))
+            print('  frac_l2 = ', frac_l2)
+            print('  frac_peak = ', frac_peak)
+            assert frac_l2 < 0.01
+            assert frac_peak < 0.01
+
+        # Exercise the chromatic reflux branch through a real fit() call.
+        psf.fit(stars, logger=logger)
+        assert psf.niter == 1
+        assert psf.nremoved == 0
+        for star in psf.stars:
+            print(star.fit.flux, star.fit.center)
+            assert flux/2 < star.fit.flux < flux*2  # < factor of 2 change from reflux.
+            assert 0 < np.abs(star.fit.center[0]) < 0.05  # centroid shift should be small.
+            assert 0 < np.abs(star.fit.center[1]) < 0.05
+
+
 @timer
 def test_bilinear_vs_five_point():
     """Exercise practical fit behavior and quantify bilinear vs five-point differences.
@@ -709,15 +861,13 @@ def make_star(x, y):
         # Note: GalSim internally does bilinear interpolation of the nominal aberrations.
         # This is not the same bilinear that we do in Piff -- we do either bilinear or
         # 5-point interpolation of the profiles.  These are not equivalent.
-        model_bilinear = RomanOpticalModel(
-            filter=filt,
+        model_bilinear = make_roman_model(
             chromatic=False,
             max_zernike=6,
             nominal_interp='bilinear',
             aberration_prior_sigma=1.0e6,
         )
-        model_five = RomanOpticalModel(
-            filter=filt,
+        model_five = make_roman_model(
             chromatic=False,
             max_zernike=6,
             nominal_interp='five_point',
@@ -791,8 +941,7 @@ def test_fit_linear():
             0.0, 0.0, 0.0,
         ])
         for nominal_interp in ['bilinear', 'five_point']:
-            model = piff.roman.RomanOpticalModel(
-                filter='H158',
+            model = make_roman_model(
                 chromatic=False,
                 max_zernike=6,
                 aberration_interp='linear',
@@ -812,7 +961,7 @@ def test_fit_linear():
             stars = [model.initialize(star) for star in stars]
             for star in stars:
                 prof = model.getProfile(truth_params, star=star)
-                model._draw_profile_to_image(prof * star.fit.flux, star.image, star.image_pos)
+                model._draw_profile_to_image(prof * star.fit.flux, star.image, star.image_pos, star)
 
             for _ in range(4):
                 stars = model.fit_many(stars)
@@ -842,8 +991,7 @@ def test_fit_linear_gradient():
     """Check linear-mode recovery when aberrations vary as a + b x + c y across the SCA.
     """
     with fast_pupil_bin():
-        model = piff.roman.RomanOpticalModel(
-            filter='H158',
+        model = make_roman_model(
             chromatic=False,
             max_zernike=6,
             aberration_interp='linear',
@@ -948,8 +1096,7 @@ def test_fit_linear_nproc():
     """Check linear-mode fit_many behavior with multiprocessing enabled.
     """
     with fast_pupil_bin():
-        model = piff.roman.RomanOpticalModel(
-            filter='H158',
+        model = make_roman_model(
             chromatic=False,
             max_zernike=6,
             aberration_interp='linear',
@@ -982,7 +1129,7 @@ def test_fit_linear_nproc():
         ])
         for star in stars:
             prof = model.getProfile(truth_params, star=star)
-            model._draw_profile_to_image(prof * star.fit.flux, star.image, star.image_pos)
+            model._draw_profile_to_image(prof * star.fit.flux, star.image, star.image_pos, star)
 
         for _ in range(4):
             stars = model.fit_many(stars)
@@ -991,21 +1138,21 @@ def test_fit_linear_nproc():
             np.testing.assert_allclose(star.fit.params[3:], truth_params[3:], atol=5.e-9, rtol=0.0)
 
 
-
 @timer
 def test_optics_convert_funcs():
     """Check aberration recovery when fitting with a nontrivial convert_func (profile shear).
     """
     with fast_pupil_bin():
+        bandpass = galsim.roman.getBandpasses()['H158']
         psf = piff.PSF.process(
             {
                 'type': 'RomanOptics',
-                'filter': 'H158',
                 'chromatic': False,
                 'max_zernike': 6,
                 'aberration_prior_sigma': 1.0e6,
             }
         )
+        psf.set_context(None, None, bandpass)
         logger = piff.config.setup_logger()
         stars = [
             piff.Star.makeTarget(
@@ -1034,7 +1181,7 @@ def apply_shear(prof):
         for star in stars:
             # True profile is sheared version of optical PSF.
             prof = psf.model.getProfile(truth_params, star=star).shear(g1=0.01, g2=-0.005)
-            psf.model._draw_profile_to_image(prof, star.image, star.image_pos)
+            psf.model._draw_profile_to_image(prof, star.image, star.image_pos, star)
 
         for _ in range(3):
             stars, nremoved = psf.single_iteration(
@@ -1060,7 +1207,7 @@ def test_sca_interp():
     with fast_pupil_bin():
         # Start with some basic exercises of interp machinery.
         psf = piff.PSF.process(
-            {'type': 'RomanOptics', 'filter': 'H158', 'chromatic': False, 'max_zernike': 6}
+            {'type': 'RomanOptics', 'chromatic': False, 'max_zernike': 6}
         )
         assert type(psf.interp).__name__ == 'RomanSCAInterp'
 
@@ -1119,7 +1266,6 @@ def test_sca_interp():
         psf_global = piff.PSF.process(
             {
                 'type': 'RomanOptics',
-                'filter': 'H158',
                 'chromatic': False,
                 'max_zernike': 6,
                 'aberration_interp': 'global',
@@ -1141,7 +1287,6 @@ def test_sca_interp():
         psf = piff.PSF.process(
             {
                 'type': 'RomanOptics',
-                'filter': 'H158',
             }
         )
         psf.write(fn)
@@ -1159,6 +1304,7 @@ def test_sca_interp():
     test_aberration_prior()
     test_fit_many()
     test_fit_many_nproc()
+    test_chromatic()
     test_bilinear_vs_five_point()
     test_fit_linear()
     test_fit_linear_gradient()
diff --git a/tests/test_simple.py b/tests/test_simple.py
index c325981a..7b9707f6 100644
--- a/tests/test_simple.py
+++ b/tests/test_simple.py
@@ -269,7 +269,7 @@ def test_single_image():
         },
         'output' : { 'file_name' : psf_file },
     }
-    orig_stars, wcs, pointing = piff.Input.process(config['input'], logger)
+    orig_stars, wcs, pointing, _ = piff.Input.process(config['input'], logger)
 
     # Use a SimplePSF to process the stars data this time.
     interp = piff.Mean()
@@ -873,7 +873,7 @@ def test_load_images():
                'cat_file_name': cat_file,
                'sky': 10
              }
-    orig_stars, wcs, pointing = piff.Input.process(config, logger)
+    orig_stars, wcs, pointing, _ = piff.Input.process(config, logger)
 
     # Fit these with a simple Mean, Gaussian
     model = piff.Gaussian()
diff --git a/tests/test_sizemag.py b/tests/test_sizemag.py
index 4146fb79..e6fb4472 100644
--- a/tests/test_sizemag.py
+++ b/tests/test_sizemag.py
@@ -64,7 +64,7 @@ def test_smallbright():
         'type': 'SmallBright'
     }
 
-    objects, _, _ = piff.Input.process(config['input'], logger=logger)
+    objects, _, _, _ = piff.Input.process(config['input'], logger=logger)
     stars = piff.Select.process(config['select'], objects, logger=logger)
 
     # This does a pretty decent job actually.  Finds 88 stars.
@@ -230,7 +230,7 @@ def test_sizemag():
         'type': 'SizeMag'
     }
 
-    objects, _, _ = piff.Input.process(config['input'], logger=logger)
+    objects, _, _, _ = piff.Input.process(config['input'], logger=logger)
     stars = piff.Select.process(config['select'], objects, logger=logger)
 
     # This finds more stars than the simple SmallBright selector found.
@@ -328,7 +328,7 @@ def test_sizemag():
             # Use all the input objects
         }
     }
-    objects, _, _ = piff.Input.process(config['input'], logger=logger)
+    objects, _, _, _ = piff.Input.process(config['input'], logger=logger)
     stars = piff.Select.process(config['select'], objects, logger=logger)
 
     # Check pathologically few input objects.
diff --git a/tests/test_stats.py b/tests/test_stats.py
index ca48015b..e020e41d 100644
--- a/tests/test_stats.py
+++ b/tests/test_stats.py
@@ -376,7 +376,7 @@ def test_rhostats_config():
     max_sep = 100
     bin_size = 0.1
     psf = piff.read(psf_file)
-    orig_stars, wcs, pointing = piff.Input.process(config['input'], logger)
+    orig_stars, wcs, pointing, _ = piff.Input.process(config['input'], logger)
     stats = piff.RhoStats(min_sep=min_sep, max_sep=max_sep, bin_size=bin_size)
     with np.testing.assert_raises(RuntimeError):
         stats.write('dummy')  # Cannot write before compute
@@ -471,7 +471,7 @@ def test_shapestats_config():
     # Test ShapeHistStats directly
     psf = piff.read(psf_file)
     shapeStats = piff.ShapeHistStats(nbins=5)  # default is sqrt(nstars)
-    orig_stars, wcs, pointing = piff.Input.process(config['input'], logger)
+    orig_stars, wcs, pointing, _ = piff.Input.process(config['input'], logger)
     with np.testing.assert_raises(RuntimeError):
         shapeStats.write()  # Cannot write before compute
     shapeStats.compute(psf, orig_stars)
@@ -554,7 +554,7 @@ def test_starstats_config():
     # check default nplot
     psf = piff.read(psf_file)
     starStats = piff.StarStats(include_ave=False)
-    orig_stars, wcs, pointing = piff.Input.process(config['input'], logger=logger)
+    orig_stars, wcs, pointing, _ = piff.Input.process(config['input'], logger=logger)
     orig_stars = piff.Select.process(config['select'], orig_stars, logger=logger)
     with np.testing.assert_raises(RuntimeError):
         starStats.write()  # Cannot write before compute
@@ -742,7 +742,7 @@ def test_hsmcatalog():
 
     # Use class directly, rather than through config.
     psf = piff.PSF.read(psf_file)
-    stars, _, _ = piff.Input.process(config['input'])
+    stars, _, _, _ = piff.Input.process(config['input'])
     stars = piff.Select.process(config['select'], stars)
     hsmcat = piff.stats.HSMCatalogStats()
     with np.testing.assert_raises(RuntimeError):
diff --git a/tests/test_wcs.py b/tests/test_wcs.py
index 2c7beee5..4c10e669 100644
--- a/tests/test_wcs.py
+++ b/tests/test_wcs.py
@@ -372,6 +372,21 @@ def test_single():
             print('max diff = ',np.max(np.abs(im.array - star.data.image.array)))
             np.testing.assert_almost_equal(im.array, star.data.image.array)
 
+    # Make sure old pattern of giving wcs, pointing in fit() still works.
+    objects, wcs_dict, pointing, bandpass = piff.Input.process(config['input'], logger=logger)
+    assert bandpass is None
+    with CaptureLog(level=2) as cl:
+        psf_old = piff.PSF.process(config['psf'], logger=cl.logger)
+        psf_old.fit(objects, wcs_dict, pointing, logger=cl.logger)
+    assert "wcs should now be set with set_context" in cl.output
+
+    for chipnum, data in [(1, data1), (2, data2)]:
+        x = data['x'][0]
+        y = data['y'][0]
+        im = psf.draw(x, y, chipnum=chipnum)
+        im_old = psf_old.draw(x, y, chipnum=chipnum)
+        np.testing.assert_almost_equal(im_old.array, im.array)
+
     # SingleChip doesn't use code paths that hit these functions, but they are officially
     # required methods for PSF classes, so just check them directly.
     assert psf.fit_center is True