diff --git a/piff/__init__.py b/piff/__init__.py
index f8c49de1..6da1e72c 100644
--- a/piff/__init__.py
+++ b/piff/__init__.py
@@ -113,6 +113,7 @@
 
 # Composite PSFs
 from .sumpsf import SumPSF
+from .convolvepsf import ConvolvePSF
 
 # Leave these in their own namespaces
 from . import util
diff --git a/piff/convolvepsf.py b/piff/convolvepsf.py
new file mode 100644
index 00000000..5fe8a9cb
--- /dev/null
+++ b/piff/convolvepsf.py
@@ -0,0 +1,314 @@
+# Copyright (c) 2016 by Mike Jarvis and the other collaborators on GitHub at
+# https://github.com/rmjarvis/Piff  All rights reserved.
+#
+# Piff is free software: Redistribution and use in source and binary forms
+# with or without modification, are permitted provided that the following
+# conditions are met:
+#
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the disclaimer given in the accompanying LICENSE
+#    file.
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the disclaimer given in the documentation
+#    and/or other materials provided with the distribution.
+
+"""
+.. module:: psf
+"""
+
+import numpy as np
+import galsim
+
+from .psf import PSF
+from .util import write_kwargs, read_kwargs
+from .star import Star, StarFit
+from .outliers import Outliers
+
+class ConvolvePSF(PSF):
+    """A PSF class that is the Convolution of two or more other PSF types.
+
+    A ConvolvePSF is built from an ordered list of other PSFs.
+
+    When fitting the convolution, the default pattern is that all components after the first one
+    are initialized as (approximately) a delta function, and the first component is fit as usual,
+    but just using a single iteration of the fit. Then the residuals of this model are fit using
+    the second component, and so on. Once all components are fit, outliers may be rejected, and
+    then the process is iterated.
+
+    This pattern can be tweaked somewhat using the initialization options available to
+    PSF models. If a component should be initialized to something other than a delta-function.
+    then one should explicitly set it.
+
+    :param components:  A list of PSF instances defining the components to be convolved.
+    :param outliers:    Optionally, an Outliers instance used to remove outliers.
+                        [default: None]
+    :param chisq_thresh: Change in reduced chisq at which iteration will terminate.
+                        [default: 0.1]
+    :param max_iter:    Maximum number of iterations to try. [default: 30]
+    """
+    _type_name = 'Convolve'
+
+    def __init__(self, components, outliers=None, chisq_thresh=0.1, max_iter=30):
+        self.components = components
+        self.outliers = outliers
+        self.chisq_thresh = chisq_thresh
+        self.max_iter = max_iter
+        self.kwargs = {
+            # If components is a list, mark the number of components here for I/O purposes.
+            # But if it's just a number, leave it alone.
+            'components': len(components) if isinstance(components, list) else components,
+            'outliers': 0,
+            'chisq_thresh': self.chisq_thresh,
+            'max_iter': self.max_iter,
+        }
+        self.chisq = 0.
+        self.last_delta_chisq = 0.
+        self.dof = 0
+        self.nremoved = 0
+        self.niter = 0
+        self.set_num(None)
+
+    def set_num(self, num):
+        """If there are multiple components involved in the fit, set the number to use
+        for this model.
+        """
+        if isinstance(self.components, list):
+            # Then building components has been completed.
+
+            # This array keeps track of which num to use for each component.
+            self._nums = np.empty(len(self.components), dtype=int)
+            self._num_components = 0  # It might not be len(self.components) if some of them are
+                                      # in turn composite types. Figure it out below.
+
+            k1 = 0 if num is None else num
+            for k, comp in enumerate(self.components):
+                self._nums[k] = k1
+                comp.set_num(k1)
+                k1 += comp.num_components
+                self._num_components += comp.num_components
+            self._num = self._nums[0]
+        else:
+            # else components are not yet built. This will be called again when that's done.
+            self._num = None
+
+    @property
+    def num_components(self):
+        return self._num_components
+
+    @property
+    def interp_property_names(self):
+        names = set()
+        for c in self.components:
+            names.update(c.interp_property_names)
+        return names
+
+    @classmethod
+    def parseKwargs(cls, config_psf, logger):
+        """Parse the psf field of a configuration dict and return the kwargs to use for
+        initializing an instance of the class.
+
+        :param config_psf:      The psf field of the configuration dict, config['psf']
+        :param logger:          A logger object for logging debug info. [default: None]
+
+        :returns: a kwargs dict to pass to the initializer
+        """
+        from .outliers import Outliers
+
+        kwargs = {}
+        kwargs.update(config_psf)
+        kwargs.pop('type',None)
+
+        if 'components' not in kwargs:
+            raise ValueError("components field is required in psf field for type=Convolve")
+
+        # make components
+        components = kwargs.pop('components')
+        kwargs['components'] = []
+        for comp in components:
+            kwargs['components'].append(PSF.process(comp, logger=logger))
+
+        if 'outliers' in kwargs:
+            outliers = Outliers.process(kwargs.pop('outliers'), logger=logger)
+            kwargs['outliers'] = outliers
+
+        return kwargs
+
+    def setup_params(self, stars):
+        """Make sure the stars have the right shape params object
+        """
+        new_stars = []
+        for star in stars:
+            if star.fit.params is None:
+                fit = star.fit.withNew(params=[None] * self.num_components,
+                                       params_var=[None] * self.num_components)
+                star = Star(star.data, fit)
+            else:
+                assert len(star.fit.params) > self._nums[-1]
+            new_stars.append(star)
+        return new_stars
+
+    def initialize_params(self, stars, logger, default_init=None):
+        nremoved = 0
+
+        logger.debug("Initializing components of ConvolvePSF")
+
+        # First make sure the params are the right shape.
+        stars = self.setup_params(stars)
+
+        # Now initialize all the components
+        for comp in self.components:
+            stars, nremoved1 = comp.initialize_params(stars, logger, default_init=default_init)
+            nremoved += nremoved1
+            # After the first one, set default_init to 'delta'
+            default_init='delta'
+
+        # If any components are degenerate, mark the convolution as degenerate.
+        self.degenerate_points = any([comp.degenerate_points for comp in self.components])
+
+        return stars, nremoved
+
+    def single_iteration(self, stars, logger, convert_funcs, draw_method):
+        nremoved = 0  # For this iteration
+
+        # Fit each component in order
+        for k, comp in enumerate(self.components):
+            logger.info("Starting work on component %d (%s)", k, comp._type_name)
+
+            # Update the convert_funcs to add a convolution by the other components.
+            new_convert_funcs = []
+            for k, star in enumerate(stars):
+                others, method = self._getRawProfile(star, skip=comp)
+
+                if others is None:
+                    cf = convert_funcs[k] if convert_funcs is not None else None
+                elif convert_funcs is None:
+                    cf = lambda prof: galsim.Convolve(prof, others)
+                else:
+                    cf = lambda prof: galsim.Convolve(convert_funcs[k](prof), others)
+                new_convert_funcs.append(cf)
+
+            stars, nremoved1 = comp.single_iteration(stars, logger, new_convert_funcs, method)
+            nremoved += nremoved1
+
+            # Update the current models for later components
+            stars = comp.interpolateStarList(stars)
+
+        return stars, nremoved
+
+    @property
+    def fit_center(self):
+        """Whether to fit the center of the star in reflux.
+
+        This is generally set in the model specifications.
+        If all component models includes a shift, then this is False.
+        Otherwise it is True.
+        """
+        return any([comp.fit_center for comp in self.components])
+
+    @property
+    def include_model_centroid(self):
+        """Whether a model that we want to center can have a non-zero centroid during iterations.
+        """
+        return any([comp.include_model_centroid for comp in self.components])
+
+    def interpolateStarList(self, stars):
+        """Update the stars to have the current interpolated fit parameters according to the
+        current PSF model.
+
+        :param stars:       List of Star instances to update.
+
+        :returns:           List of Star instances with their fit parameters updated.
+        """
+        stars = self.setup_params(stars)
+        for comp in self.components:
+            stars = comp.interpolateStarList(stars)
+        return stars
+
+    def interpolateStar(self, star):
+        """Update the star to have the current interpolated fit parameters according to the
+        current PSF model.
+
+        :param star:        Star instance to update.
+
+        :returns:           Star instance with its fit parameters updated.
+        """
+        star, = self.setup_params([star])
+        for comp in self.components:
+            star = comp.interpolateStar(star)
+        return star
+
+    def _drawStar(self, star):
+        params = star.fit.get_params(self._num)
+        prof, method = self._getRawProfile(star)
+        prof = prof.shift(star.fit.center) * star.fit.flux
+        image = prof.drawImage(image=star.image.copy(), method=method, center=star.image_pos)
+        return Star(star.data.withNew(image=image), star.fit)
+
+    def _getRawProfile(self, star, skip=None):
+        # Get each component profile
+        profiles = []
+        methods = []
+        for comp in self.components:
+            prof, method = comp._getRawProfile(star)
+            methods.append(method)
+            if comp is not skip:
+                profiles.append(prof)
+
+        # If any components already include the pixel, then keep no_pixel for the convolution.
+        assert all([m == 'no_pixel' or m == 'auto' for m in methods])
+        if any([m == 'no_pixel' for m in methods]):
+            method = 'no_pixel'
+        else:
+            method = 'auto'
+
+        # Convolve them.
+        if len(profiles) == 0:
+            return None, method
+        elif len(profiles) == 1:
+            return profiles[0], method
+        else:
+            return galsim.Convolve(profiles), method
+
+    def _finish_write(self, fits, extname, logger):
+        """Finish the writing process with any class-specific steps.
+
+        :param fits:        An open fitsio.FITS object
+        :param extname:     The base name of the extension to write to.
+        :param logger:      A logger object for logging debug info.
+        """
+        logger = galsim.config.LoggerWrapper(logger)
+        chisq_dict = {
+            'chisq' : self.chisq,
+            'last_delta_chisq' : self.last_delta_chisq,
+            'dof' : self.dof,
+            'nremoved' : self.nremoved,
+        }
+        write_kwargs(fits, extname + '_chisq', chisq_dict)
+        logger.debug("Wrote the chisq info to extension %s",extname + '_chisq')
+        for k, comp in enumerate(self.components):
+            comp._write(fits, extname + '_' + str(k), logger=logger)
+        if self.outliers:
+            self.outliers.write(fits, extname + '_outliers')
+            logger.debug("Wrote the PSF outliers to extension %s",extname + '_outliers')
+
+    def _finish_read(self, fits, extname, logger):
+        """Finish the reading process with any class-specific steps.
+
+        :param fits:        An open fitsio.FITS object
+        :param extname:     The base name of the extension to write to.
+        :param logger:      A logger object for logging debug info.
+        """
+        chisq_dict = read_kwargs(fits, extname + '_chisq')
+        for key in chisq_dict:
+            setattr(self, key, chisq_dict[key])
+
+        ncomponents = self.components
+        self.components = []
+        for k in range(ncomponents):
+            self.components.append(PSF._read(fits, extname + '_' + str(k), logger=logger))
+        if extname + '_outliers' in fits:
+            self.outliers = Outliers.read(fits, extname + '_outliers')
+        else:
+            self.outliers = None
+        # Set up all the num's properly now that everything is constructed.
+        self.set_num(None)
diff --git a/piff/gsobject_model.py b/piff/gsobject_model.py
index 3224ec74..ae5b2b2e 100644
--- a/piff/gsobject_model.py
+++ b/piff/gsobject_model.py
@@ -166,7 +166,7 @@ def getProfile(self, params):
 
         return prof
 
-    def _resid(self, params, star, convert_func):
+    def _resid(self, params, star, convert_func, draw_method):
         """Residual function to use with least_squares.
 
         Essentially `chi` from `chisq`, but not summed over pixels yet.
@@ -176,6 +176,8 @@ def _resid(self, params, star, convert_func):
         :param convert_func:    An optional function to apply to the profile being fit before
                                 drawing it onto the image.  This is used by composite PSFs to
                                 isolate the effect of just this model component.
+        :param draw_method:     The method to use with the GalSim drawImage command to determine
+                                the residuals. [if draw_method is None, use self._method]
 
         :returns: `chi` as a flattened numpy array.
         """
@@ -203,9 +205,10 @@ def _resid(self, params, star, convert_func):
 
         if convert_func is not None:
             prof = convert_func(prof)
-
+        if draw_method is None:
+            draw_method = self._method
         try:
-            prof.drawImage(model_image, method=self._method, center=image_pos)
+            prof.drawImage(model_image, method=draw_method, center=image_pos)
         except galsim.GalSimError:
             # Declare this a *bad* residual if GalSim ran into some kind of error.
             # (Most typically a GalSimFFTSizeError from getting to a place with bad stepk/maxk.)
@@ -241,7 +244,7 @@ def _get_params(self, star):
 
         return np.array([flux, du, dv, scale, g1, g2])
 
-    def least_squares_fit(self, star, logger=None, convert_func=None):
+    def least_squares_fit(self, star, logger=None, convert_func=None, draw_method=None):
         """Fit parameters of the given star using least-squares minimization.
 
         :param star:            A Star to fit.
@@ -249,6 +252,8 @@ def least_squares_fit(self, star, logger=None, convert_func=None):
         :param convert_func:    An optional function to apply to the profile being fit before
                                 drawing it onto the image.  This is used by composite PSFs to
                                 isolate the effect of just this model component. [default: None]
+        :param draw_method:     The method to use with galsim.DrawImage to determine the residuals.
+                                [default: 'auto' if include_pixel = True, else 'no_pixel']
 
         :returns: (flux, dx, dy, scale, g1, g2, flag)
         """
@@ -262,7 +267,8 @@ def least_squares_fit(self, star, logger=None, convert_func=None):
         if params[3]**2 < star.data.local_wcs.pixelArea():
             params[3] = star.data.local_wcs.pixelArea()**0.5
 
-        results = scipy.optimize.least_squares(self._resid, params, args=(star,convert_func),
+        results = scipy.optimize.least_squares(self._resid, params,
+                                               args=(star,convert_func,draw_method),
                                                **self._scipy_kwargs)
         #logger.debug(results)
         logger.debug("results = %s",results.x)
@@ -274,7 +280,7 @@ def least_squares_fit(self, star, logger=None, convert_func=None):
 
         return flux, du, dv, scale, g1, g2, var
 
-    def fit(self, star, fastfit=None, logger=None, convert_func=None):
+    def fit(self, star, fastfit=None, logger=None, convert_func=None, draw_method=None):
         """Fit the image either using HSM or least-squares minimization.
 
         If ``fastfit`` is True, then the galsim.hsm module will be used to estimate the
@@ -289,7 +295,11 @@ def fit(self, star, fastfit=None, logger=None, convert_func=None):
         :param logger:          A logger object for logging debug info. [default: None]
         :param convert_func:    An optional function to apply to the profile being fit before
                                 drawing it onto the image.  This is used by composite PSFs to
-                                isolate the effect of just this model component. [default: None]
+                                isolate the effect of just this model component. [default: None;
+                                if provided, fastfit is ignored and taken to be False.]
+        :param draw_method:     The method to use with the GalSim drawImage command to determine
+                                the residuals. [default: None, which uses 'auto' if
+                                include_pixel = True, else 'no_pixel']
 
         :returns: a new Star with the fitted parameters in star.fit
         """
@@ -305,8 +315,8 @@ def fit(self, star, fastfit=None, logger=None, convert_func=None):
         if fastfit:
             flux, du, dv, scale, g1, g2, var = self.moment_fit(star, logger=logger)
         else:
-            flux, du, dv, scale, g1, g2, var = self.least_squares_fit(star, logger=logger,
-                                                                      convert_func=convert_func)
+            flux, du, dv, scale, g1, g2, var = self.least_squares_fit(
+                star, logger=logger, convert_func=convert_func, draw_method=draw_method)
 
         # Make a StarFit object with these parameters
         params = [scale, g1, g2]
@@ -330,9 +340,10 @@ def fit(self, star, fastfit=None, logger=None, convert_func=None):
 
         if convert_func is not None:
             prof = convert_func(prof)
-
+        if draw_method is None:
+            draw_method = self._method
         try:
-            prof.drawImage(model_image, method=self._method, center=star.image_pos)
+            prof.drawImage(model_image, method=draw_method, center=star.image_pos)
         except galsim.GalSimError:
             # If we get the exception here, then set a large chisq, so it gets outlier rejected.
             chisq = 1.e300
diff --git a/piff/model.py b/piff/model.py
index 31944916..393a4057 100644
--- a/piff/model.py
+++ b/piff/model.py
@@ -134,7 +134,7 @@ def fit(self, star, convert_func=None):
         """
         raise NotImplementedError("Derived classes must define the fit function")
 
-    def draw(self, star, copy_image=True, center=None):
+    def draw(self, star, copy_image=True):
         """Draw the model on the given image.
 
         :param star:        A Star instance with the fitted parameters to use for drawing and a
@@ -142,8 +142,6 @@ def draw(self, star, copy_image=True, center=None):
         :param copy_image:  If False, will use the same image object.
                             If True, will copy the image and then overwrite it.
                             [default: True]
-        :param center:      An optional tuple (x,y) location for where to center the drawn profile
-                            in the image. [default: None, which draws at the star's location.]
 
         :returns: a new Star instance with the data field having an image of the drawn model.
         """
@@ -153,11 +151,7 @@ def draw(self, star, copy_image=True, center=None):
             image = star.image.copy()
         else:
             image = star.image
-        if center is None:
-            center = star.image_pos
-        else:
-            center = galsim.PositionD(*center)
-        prof.drawImage(image, method=self._method, center=center)
+        prof.drawImage(image, method=self._method, center=star.image_pos)
         return Star(star.data.withNew(image=image), star.fit)
 
     def write(self, fits, extname):
diff --git a/piff/optical_model.py b/piff/optical_model.py
index ac3a98c0..8ef0e9c7 100644
--- a/piff/optical_model.py
+++ b/piff/optical_model.py
@@ -121,6 +121,8 @@ class Optical(Model):
     :param oversampling     Numerical factor by which to oversample the FFT [default: 1]
     :param pupil_plane_im:  The name of a file containing the pupil plane image to use instead
                             of creating one from obscuration, struts, etc. [default: None]
+    :param base_aberrations: A list of aberrations to which any fitted aberrations are added.
+                            [default: None]
 
     The next two parameters are specific to our treatment of the Optical PSF. If given,
     they are used to create an additional galsim.PhaseScreen describing the mirror's own
@@ -161,7 +163,7 @@ class Optical(Model):
     """
     _type_name = 'Optical'
     _method = 'auto'
-    _centered = True
+    _centered = False
     _model_can_be_offset = False
 
     def __init__(self, template=None, gsparams=None, atmo_type='VonKarman', logger=None, **kwargs):
@@ -217,7 +219,7 @@ def __init__(self, template=None, gsparams=None, atmo_type='VonKarman', logger=N
                 self.gsparams = galsim.GSParams(**gsparams_kwargs)
 
         # Check that no unexpected parameters were passed in:
-        other_keys = ('sigma', 'mirror_figure_im', 'mirror_figure_scale')
+        other_keys = ('sigma', 'mirror_figure_im', 'mirror_figure_scale', 'base_aberrations')
         extra_kwargs = [k for k in kwargs if k not in opt_keys + gsparams_keys + other_keys]
         if len(extra_kwargs) > 0:
             raise TypeError('__init__() got an unexpected keyword argument %r'%extra_kwargs[0])
@@ -271,9 +273,13 @@ def __init__(self, template=None, gsparams=None, atmo_type='VonKarman', logger=N
         self.idx_g1 = self.idx_L0 + 1
         self.idx_g2 = self.idx_g1 + 1
         self.param_len = self.idx_g2 + 1
+        self.base_aberrations = np.zeros(self.nZ+1)
+        base_aberrations = kwargs.get('base_aberrations', None)
+        if base_aberrations is not None:
+            self.base_aberrations[0:len(base_aberrations)] += base_aberrations
 
     # fit is currently a DUMMY routine
-    def fit(self, star):
+    def fit(self, star, logger=None, convert_func=None, draw_method=None):
         """Warning: This method just updates the fit with the chisq and dof!
 
         :param star:    A Star instance
@@ -294,6 +300,21 @@ def fit(self, star):
         return Star(star.data,
                     star.fit.newParams(params, params_var=var, num=self._num, chisq=chisq, dof=dof))
 
+    def initialize(self, star, logger=None, default_init=None):
+        """Initialize the given star's fit parameters.
+
+        :param star:            The Star to initialize.
+        :param logger:          A logger object for logging debug info. [default: None]
+        :param default_init:    Ignored. [default: None]
+
+        :returns: a new initialized Star.
+        """
+        params = self.kwargs_to_params()
+        params_var = np.zeros_like(params)
+        fit = star.fit.newParams(params, params_var=params_var, num=self._num)
+        star = Star(star.data, fit)
+        return star
+
     @lru_cache(maxsize=8)
     def getOptics(self, zernike_coeff):
         """Get the strictly Optics PSF component.  Use two phase screens, one from
@@ -407,6 +428,9 @@ def getProfile(self,params=None,zernike_coeff=None,r0=None,L0=None,g1=0.,g2=0.):
         # otherwise parameters are in getProfile argument list
         if params is not None :
             *zernike_coeff,r0,L0,g1,g2 = params
+        aberrations = self.base_aberrations.copy()
+        if zernike_coeff is not None:
+            aberrations[:len(zernike_coeff)] += zernike_coeff
 
         # list of PSF components
         prof = []
@@ -419,13 +443,15 @@ def getProfile(self,params=None,zernike_coeff=None,r0=None,L0=None,g1=0.,g2=0.):
             prof.append(gaussian)
 
         # atmospheric kernel
-        atmopsf = self.getAtmosphere(r0, L0, g1, g2)
-        prof.append(atmopsf)
+        if self.atmo_type not in ['None', None]:
+            atmopsf = self.getAtmosphere(r0, L0, g1, g2)
+            prof.append(atmopsf)
 
         # optics
-        prof.append(self.getOptics(tuple(zernike_coeff)))
-
-        # convolve
-        prof = galsim.Convolve(prof,gsparams=self.gsparams)
+        optics = self.getOptics(tuple(aberrations))
 
-        return prof
+        if len(prof) == 0:
+            return optics
+        else:
+            prof.append(optics)
+            return galsim.Convolve(prof, gsparams=self.gsparams)
diff --git a/piff/pixelgrid.py b/piff/pixelgrid.py
index 3d7a525c..d275caf6 100644
--- a/piff/pixelgrid.py
+++ b/piff/pixelgrid.py
@@ -158,7 +158,7 @@ def initialize(self, star, logger=None, default_init=None):
         fit = star.fit.newParams(params=params, num=self._num)
         return Star(star.data, fit)
 
-    def fit(self, star, logger=None, convert_func=None):
+    def fit(self, star, logger=None, convert_func=None, draw_method=None):
         """Fit the Model to the star's data to yield iterative improvement on its PSF parameters
         and uncertainties.
 
@@ -167,9 +167,14 @@ def fit(self, star, logger=None, convert_func=None):
         :param convert_func:    An optional function to apply to the profile being fit before
                                 drawing it onto the image.  This is used by composite PSFs to
                                 isolate the effect of just this model component. [default: None]
+        :param draw_method:     The method to use with the GalSim drawImage command to determine
+                                the residuals. [PixelGrid always uses 'no_pixel'; this parameter
+                                is only present for API compatibility.  It must be either None
+                                or 'no_pixel'.]
 
         :returns: a new Star instance with updated fit information
         """
+        assert draw_method in (None, 'no_pixel')
         logger = galsim.config.LoggerWrapper(logger)
         # Get chisq Taylor expansion for linearized model
         star1 = self.chisq(star, logger=logger, convert_func=convert_func)
@@ -216,7 +221,7 @@ def fit(self, star, logger=None, convert_func=None):
         self.normalize(star)
         return star
 
-    def chisq(self, star, logger=None, convert_func=None):
+    def chisq(self, star, logger=None, convert_func=None, draw_method=None):
         """Calculate dependence of chi^2 = -2 log L(D|p) on PSF parameters for single star.
         as a quadratic form chi^2 = dp^T AT A dp - 2 bT A dp + chisq,
         where dp is the *shift* from current parameter values.  Returned Star
@@ -228,9 +233,14 @@ def chisq(self, star, logger=None, convert_func=None):
         :param convert_func:    An optional function to apply to the profile being fit before
                                 drawing it onto the image.  This is used by composite PSFs to
                                 isolate the effect of just this model component. [default: None]
+        :param draw_method:     The method to use with the GalSim drawImage command to determine
+                                the residuals. [PixelGrid always uses 'no_pixel'; this parameter
+                                is only present for API compatibility.  It must be either None
+                                or 'no_pixel'.]
 
         :returns: a new Star instance with updated fit parameters. (esp. A,b)
         """
+        assert draw_method in (None, 'no_pixel')
         logger = galsim.config.LoggerWrapper(logger)
         logger.debug('Start chisq function')
         logger.debug('initial params = %s',star.fit.get_params(self._num))
@@ -253,6 +263,7 @@ def chisq(self, star, logger=None, convert_func=None):
 
         if convert_func is not None:
             prof = convert_func(prof)
+            logger.debug(f'converted prof = {prof}')
 
         # Adjust the image now so that we can draw it in pixel coordinates, rather than letting
         # galsim.drawImage do the adjustment each time for each component, which would be
@@ -263,10 +274,16 @@ def chisq(self, star, logger=None, convert_func=None):
         image.wcs = galsim.PixelScale(1.0)
 
         # Draw the profile.
-        # Be careful here.  If prof was converted to something that isn't analytic in
-        # real space, this will need to be _drawFFT.
         prof = star.data.local_wcs.profileToImage(prof, offset=offset)
-        prof._drawReal(image)
+        if (convert_func is None or
+                (prof.is_analytic_x and not isinstance(prof, galsim.Convolution))):
+            draw_real = True
+            prof._drawReal(image)
+        else:
+            # If the convert func turns this into a Convolution, or something else that isn't
+            # analytic in real space, then use FFT drawing.
+            draw_real = False
+            prof.drawFFT(image)
         logger.debug('drawn flux = %s',image.array.sum())
         model = image.array.ravel() * star.fit.flux
 
@@ -298,27 +315,26 @@ def chisq(self, star, logger=None, convert_func=None):
             # In this case we need the basis_profile to have the right scale (rather than
             # incorporate it into the jacobian) so that convert_func will have the right size.
             basis_profile = basis_profile.dilate(self.scale)
+            logger.debug(f'basis_profile = {basis_profile}')
             # Find the net shift from the star.fit.center and the offset.
             jac = star.data.local_wcs.jacobian().inverse().getMatrix()
+            jac_det = abs(jac[0,0]*jac[1,1] - jac[0,1]*jac[1,0])
             dx = jac[0,0]*u0 + jac[0,1]*v0 + offset.x
             dy = jac[1,0]*u0 + jac[1,1]*v0 + offset.y
             du = du.ravel() * self.scale
             dv = dv.ravel() * self.scale
             for k, duk, dvk in zip(range(self._nparams), du,dv):
-                # This implementation removes most of the overhead, and is still relatively
-                # straightforward.
-                # The one wrinkle is that if the net profile is no longer analytic in real space,
-                # we'll need a call to _drawFFT rather than _drawReal.  That would be a lot slower
-                # I think, so we may want to limit convolutions in conjunction with PixelGrid.
-                # (That's less of an issue for models with few paramters.)
-
-                basis_profile_k = basis_profile._shift(duk,dvk)
-                basis_profile_k = convert_func(basis_profile_k)
-                # TODO: If the convert_func has transformed the profile into something that is
-                # not analytic_x, then we need to do the _drawFFT version here.
-                basis_profile_k._drawReal(image, jac, (dx,dy), 1.)
+                prof = basis_profile._shift(duk,dvk)
+                prof = convert_func(prof)
+                if draw_real:
+                    prof._drawReal(image, jac, (dx,dy), 1.)
+                    # Equivalent to:
+                    #prof = galsim._Transform(prof, jac, (dx,dy), 1./jac_det)
+                    #prof.drawReal(image)
+                else:
+                    prof = galsim._Transform(prof, jac, (dx,dy), 1./jac_det)
+                    prof.drawFFT(image)
                 A[:,k] = image.array.ravel()[mask]
-
         else:
             if 0:
                 # When we don't have the convert_func step, this calculation can be sped up
diff --git a/piff/psf.py b/piff/psf.py
index 6c81b149..fbb4ccd6 100644
--- a/piff/psf.py
+++ b/piff/psf.py
@@ -177,7 +177,7 @@ def initialize_params(self, stars, logger=None, default_init=None):
         # behavior is just to return the input stars list.
         return stars, 0
 
-    def single_iteration(self, stars, logger, convert_func):
+    def single_iteration(self, stars, logger, convert_funcs, draw_method):
         """Perform a single iteration of the solver.
 
         Note that some object might fail at some point in the fitting, so some objects can be
@@ -186,8 +186,11 @@ def single_iteration(self, stars, logger, convert_func):
 
         :param stars:           The list of stars to use for constraining the PSF.
         :param logger:          A logger object for logging progress.
-        :param convert_func:    An optional function to apply to the profile being fit before
-                                drawing it onto the image.
+        :param convert_funcs:   An optional list of function to apply to the profiles being fit
+                                before drawing it onto the image.  If not None, it should be the
+                                same length as stars.
+        :param draw_method:     The method to use with the GalSim drawImage command. If None,
+                                use the default method for the PSF model being fit.
 
         :returns: an updated list of all_stars, nremoved
         """
@@ -201,7 +204,13 @@ def fit_center(self):
         If all component models includes a shift, then this is False.
         Otherwise it is True.
         """
-        raise NotImplementedError("Derived classes must define the single_iteration function")
+        raise NotImplementedError("Derived classes must define the fit_center property")
+
+    @property
+    def include_model_centroid(self):
+        """Whether a model that we want to center can have a non-zero centroid during iterations.
+        """
+        raise NotImplementedError("Derived classes must define the include_model_centroid property")
 
     def reflux(self, star, logger=None):
         """Fit the PSF to the star's data, varying only the flux (and
@@ -284,12 +293,12 @@ def reflux(self, star, logger=None):
             AtA = Atw.dot(At.T)
             Atb = Atw.dot(resid)
             x = np.linalg.solve(AtA, Atb)
-            logger.debug('    centroid shift = %s,%s', x[0], x[1])
+            logger.debug('    centroid shift = %s,%s', x[1], x[2])
 
             # Extract the values we want.
             df, duc, dvc = x
 
-            if psf_prof.centroid != galsim.PositionD(0,0):
+            if self.include_model_centroid and psf_prof.centroid != galsim.PositionD(0,0):
                 # 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.
@@ -357,7 +366,7 @@ def remove_outliers(self, stars, iteration, logger):
             nremoved = 0
         return stars, nremoved
 
-    def fit(self, stars, wcs, pointing, logger=None, convert_func=None):
+    def fit(self, stars, wcs, pointing, 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.
@@ -365,9 +374,12 @@ def fit(self, stars, wcs, pointing, logger=None, convert_func=None):
         :param pointing:        A galsim.CelestialCoord object giving the telescope pointing.
                                 [Note: pointing should be None if the WCS is not a CelestialWCS]
         :param logger:          A logger object for logging debug info. [default: None]
-        :param convert_func:    An optional function to apply to the profile being fit before
-                                drawing it onto the image.  This is used by composite PSFs to
-                                isolate the effect of just this model component. [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
+                                to isolate the effect of just this model component.  If provided,
+                                it should be the same length as stars. [default: None]
+        :param draw_method:     The method to use with the GalSim drawImage command. If not given,
+                                use the default method for the PSF model being fit. [default: None]
         """
         logger = galsim.config.LoggerWrapper(logger)
 
@@ -389,7 +401,7 @@ def fit(self, stars, wcs, pointing, logger=None, convert_func=None):
 
             # Run a single iteration of the fitter.
             # Different PSF types do different things here.
-            stars, iter_nremoved = self.single_iteration(stars, logger, convert_func)
+            stars, iter_nremoved = self.single_iteration(stars, logger, convert_funcs, draw_method)
 
             # Update estimated poisson noise
             signals = self.drawStarList(stars)
@@ -659,7 +671,7 @@ def drawStarList(self, stars, copy_image=True):
             stars = self.interpolateStarList(stars)
         return [self._drawStar(star) for star in stars]
 
-    def drawStar(self, star, copy_image=True, center=None):
+    def drawStar(self, star, copy_image=True):
         """Generate PSF image for a given star.
 
         .. note::
@@ -677,8 +689,6 @@ def drawStar(self, star, copy_image=True, center=None):
         :param copy_image:  If False, will use the same image object.
                             If True, will copy the image and then overwrite it.
                             [default: True]
-        :param center:      An optional tuple (x,y) location for where to center the drawn profile
-                            in the image. [default: None, which draws at the star's location.]
 
         :returns:           Star instance with its image filled with rendered PSF
         """
@@ -690,9 +700,9 @@ def drawStar(self, star, copy_image=True, center=None):
         if star.fit is None or star.fit.get_params(self._num) is None:
             star = self.interpolateStar(star)
         # Render the image
-        return self._drawStar(star, center=center)
+        return self._drawStar(star)
 
-    def _drawStar(self, star, center=None):
+    def _drawStar(self, star):
         # Derived classes may choose to override any of the above functions
         # But they have to at least override this one and interpolateStar to implement
         # their actual PSF model.
diff --git a/piff/simplepsf.py b/piff/simplepsf.py
index 87ced184..995b1f36 100644
--- a/piff/simplepsf.py
+++ b/piff/simplepsf.py
@@ -157,7 +157,7 @@ def initialize_params(self, stars, logger=None, default_init=None):
 
         return stars, nremoved
 
-    def single_iteration(self, stars, logger, convert_func):
+    def single_iteration(self, stars, logger, convert_funcs, draw_method):
 
         # Perform the fit or compute design matrix as appropriate using just non-reserve stars
         fit_fn = self.model.chisq if self.quadratic_chisq else self.model.fit
@@ -165,10 +165,12 @@ def single_iteration(self, stars, logger, convert_func):
         nremoved = 0  # For this iteration
         use_stars = []  # Just the stars we want to use for fitting.
         all_stars = []  # All the stars (with appropriate flags as necessary)
-        for star in stars:
+        for k, star in enumerate(stars):
             if not star.is_flagged and not star.is_reserve:
                 try:
-                    star = fit_fn(star, logger=logger, convert_func=convert_func)
+                    convert_func = None if convert_funcs is None else convert_funcs[k]
+                    star = fit_fn(star, logger=logger, convert_func=convert_func,
+                                  draw_method=draw_method)
                     use_stars.append(star)
                 except Exception as e:
                     logger.warning("Failed fitting star at %s.", star.image_pos)
@@ -200,6 +202,12 @@ def fit_center(self):
         """
         return self.model._centered
 
+    @property
+    def include_model_centroid(self):
+        """Whether a model that we want to center can have a non-zero centroid during iterations.
+        """
+        return self.model._centered and self.model._model_can_be_offset
+
     def interpolateStarList(self, stars):
         """Update the stars to have the current interpolated fit parameters according to the
         current PSF model.
@@ -225,8 +233,8 @@ def interpolateStar(self, star):
         self.model.normalize(star)
         return star
 
-    def _drawStar(self, star, center=None):
-        return self.model.draw(star, center=center)
+    def _drawStar(self, star):
+        return self.model.draw(star)
 
     def _getRawProfile(self, star):
         return self.model.getProfile(star.fit.get_params(self._num)), self.model._method
diff --git a/piff/singlechip.py b/piff/singlechip.py
index df08635b..7dd5bf89 100644
--- a/piff/singlechip.py
+++ b/piff/singlechip.py
@@ -24,7 +24,7 @@
 from .util import write_kwargs, read_kwargs, make_dtype, adjust_value, run_multi
 
 # Used by SingleChipPSF.fit
-def single_chip_run(chipnum, single_psf, stars, wcs, pointing, convert_func, logger):
+def single_chip_run(chipnum, single_psf, stars, wcs, pointing, convert_funcs, draw_method, logger):
     # Make a copy of single_psf for each chip
     psf_chip = copy.deepcopy(single_psf)
 
@@ -34,7 +34,8 @@ def single_chip_run(chipnum, single_psf, stars, wcs, pointing, convert_func, log
 
     # Run the psf_chip fit function using this stars and wcs (and the same pointing)
     logger.warning("Building solution for chip %s with %d stars", chipnum, len(stars_chip))
-    psf_chip.fit(stars_chip, wcs_chip, pointing, logger=logger, convert_func=convert_func)
+    psf_chip.fit(stars_chip, wcs_chip, pointing, logger=logger, convert_funcs=convert_funcs,
+                 draw_method=draw_method)
 
     return psf_chip
 
@@ -92,7 +93,7 @@ def parseKwargs(cls, config_psf, logger):
 
         return { 'single_psf' : single_psf, 'nproc' : nproc }
 
-    def fit(self, stars, wcs, pointing, logger=None, convert_func=None):
+    def fit(self, stars, wcs, pointing, 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.
@@ -100,9 +101,12 @@ def fit(self, stars, wcs, pointing, logger=None, convert_func=None):
         :param pointing:        A galsim.CelestialCoord object giving the telescope pointing.
                                 [Note: pointing should be None if the WCS is not a CelestialWCS]
         :param logger:          A logger object for logging debug info. [default: None]
-        :param convert_func:    An optional function to apply to the profile being fit before
-                                drawing it onto the image.  This is used by composite PSFs to
-                                isolate the effect of just this model component. [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
+                                to isolate the effect of just this model component.  If provided,
+                                it should be the same length as stars. [default: None]
+        :param draw_method:     The method to use with the GalSim drawImage command. If not given,
+                                use the default method for the PSF model being fit. [default: None]
         """
         logger = galsim.config.LoggerWrapper(logger)
         self.wcs = wcs
@@ -110,7 +114,7 @@ def fit(self, stars, wcs, pointing, logger=None, convert_func=None):
         self.psf_by_chip = {}
 
         chipnums = list(wcs.keys())
-        args = [(chipnum, self.single_psf, stars, wcs, pointing, convert_func)
+        args = [(chipnum, self.single_psf, stars, wcs, pointing, convert_funcs, draw_method)
                 for chipnum in chipnums]
 
         output = run_multi(single_chip_run, self.nproc, raise_except=False,
@@ -139,6 +143,12 @@ def fit_center(self):
         """
         return self.single_psf.fit_center
 
+    @property
+    def include_model_centroid(self):
+        """Whether a model that we want to center can have a non-zero centroid during iterations.
+        """
+        return self.single_psf.include_model_centroid
+
     def interpolateStar(self, star):
         """Update the star to have the current interpolated fit parameters according to the
         current PSF model.
@@ -152,11 +162,11 @@ def interpolateStar(self, star):
         chipnum = star['chipnum']
         return self.psf_by_chip[chipnum].interpolateStar(star)
 
-    def _drawStar(self, star, center=None):
+    def _drawStar(self, star):
         if 'chipnum' not in star.data.properties:
             raise ValueError("SingleChip requires the star to have a chipnum property")
         chipnum = star['chipnum']
-        return self.psf_by_chip[chipnum]._drawStar(star, center=center)
+        return self.psf_by_chip[chipnum]._drawStar(star)
 
     def _getProfile(self, star):
         chipnum = star['chipnum']
diff --git a/piff/sumpsf.py b/piff/sumpsf.py
index 2182209c..6a2ca868 100644
--- a/piff/sumpsf.py
+++ b/piff/sumpsf.py
@@ -167,7 +167,7 @@ def initialize_params(self, stars, logger, default_init=None):
 
         return stars, nremoved
 
-    def single_iteration(self, stars, logger, convert_func):
+    def single_iteration(self, stars, logger, convert_funcs, draw_method):
         nremoved = 0  # For this iteration
 
         # Draw the current models for each component
@@ -187,7 +187,8 @@ def single_iteration(self, stars, logger, convert_func):
                 modified_stars.append(modified_star)
 
             # Fit this component
-            new_stars, nremoved1 = comp.single_iteration(modified_stars, logger, convert_func)
+            new_stars, nremoved1 = comp.single_iteration(modified_stars, logger, convert_funcs,
+                                                         draw_method)
 
             # new_stars now has the updated fit components, but the data parts are corrupted by
             # subtracting the other components during fitting.
@@ -212,6 +213,12 @@ def fit_center(self):
         """
         return any([comp.fit_center for comp in self.components])
 
+    @property
+    def include_model_centroid(self):
+        """Whether a model that we want to center can have a non-zero centroid during iterations.
+        """
+        return any([comp.include_model_centroid for comp in self.components])
+
     def interpolateStarList(self, stars):
         """Update the stars to have the current interpolated fit parameters according to the
         current PSF model.
@@ -238,11 +245,11 @@ def interpolateStar(self, star):
             star = comp.interpolateStar(star)
         return star
 
-    def _drawStar(self, star, center=None):
+    def _drawStar(self, star):
         # Draw each component
         comp_stars = []
         for comp in self.components:
-            comp_star = comp._drawStar(star, center=center)
+            comp_star = comp._drawStar(star)
             comp_stars.append(comp_star)
 
         # Add them up.
diff --git a/tests/test_convolvepsf.py b/tests/test_convolvepsf.py
new file mode 100644
index 00000000..f90e718a
--- /dev/null
+++ b/tests/test_convolvepsf.py
@@ -0,0 +1,735 @@
+# Copyright (c) 2016 by Mike Jarvis and the other collaborators on GitHub at
+# https://github.com/rmjarvis/Piff  All rights reserved.
+#
+# Piff is free software: Redistribution and use in source and binary forms
+# with or without modification, are permitted provided that the following
+# conditions are met:
+#
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the disclaimer given in the accompanying LICENSE
+#    file.
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the disclaimer given in the documentation
+#    and/or other materials provided with the distribution.
+
+import galsim
+import numpy as np
+import piff
+import os
+import fitsio
+from functools import lru_cache
+
+from piff_test_helper import timer
+
+@lru_cache(maxsize=1)
+def make_screens():
+    # Some parameters copied from psf_wf_movie.py in GalSim repo.
+    Ellerbroek_alts = [0.0, 2.58, 5.16, 7.73, 12.89, 15.46]  # km
+    Ellerbroek_weights = [0.652, 0.172, 0.055, 0.025, 0.074, 0.022]
+    Ellerbroek_interp = galsim.LookupTable(Ellerbroek_alts, Ellerbroek_weights,
+                                           interpolant='linear')
+
+    # Use given number of uniformly spaced altitudes
+    alts = np.max(Ellerbroek_alts)*np.arange(6)/5.
+    weights = Ellerbroek_interp(alts)  # interpolate the weights
+    weights /= sum(weights)  # and renormalize
+
+    spd = []  # Wind speed in m/s
+    dirn = [] # Wind direction in radians
+    r0_500 = [] # Fried parameter in m at a wavelength of 500 nm.
+    u = galsim.UniformDeviate(1234)
+    for i in range(6):
+        spd.append(u() * 20)
+        dirn.append(u() * 360*galsim.degrees)
+        r0_500.append(0.1 * weights[i]**(-3./5))
+
+    screens = galsim.Atmosphere(r0_500=r0_500, speed=spd, direction=dirn, altitude=alts, rng=u,
+                                screen_size=102.4, screen_scale=0.1)
+
+    # Add in an optical screen
+    screens.append(galsim.OpticalScreen(diam=8, defocus=0.7, astig1=-0.8, astig2=0.7,
+                                        trefoil1=-0.6, trefoil2=0.5,
+                                        coma1=0.5, coma2=0.7, spher=0.8, obscuration=0.4))
+
+    aper = galsim.Aperture(diam=8, obscuration=0.4)
+
+    return screens, aper
+
+
+@timer
+def test_trivial_convolve1():
+    """Test the trivial case of using a Convolve of 1 component.
+    """
+    # This is essentially the same as test_single_image in test_simple.py
+
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=2)
+    else:
+        logger = piff.config.setup_logger(log_file='output/test_trivial_convolve1.log')
+
+    # Make the image
+    image = galsim.Image(2048, 2048, scale=0.26)
+
+    # Where to put the stars. (Randomish, but make sure no blending.)
+    x_list = [ 123.12, 345.98, 567.25, 1094.94, 924.15, 1532.74, 1743.11, 888.39, 1033.29, 1409.31 ]
+    y_list = [ 345.43, 567.45, 1094.32, 924.29, 1532.92, 1743.83, 888.83, 1033.19, 1409.20, 123.11 ]
+
+    # Draw a Gaussian PSF at each location on the image.
+    sigma = 1.3
+    g1 = 0.23
+    g2 = -0.17
+    psf = galsim.Gaussian(sigma=sigma).shear(g1=g1, g2=g2)
+    targets = []
+    for x,y in zip(x_list, y_list):
+        bounds = galsim.BoundsI(int(x-31), int(x+32), int(y-31), int(y+32))
+        psf.drawImage(image=image[bounds], method='no_pixel', center=(x,y))
+        targets.append(image[bounds])
+    image.addNoise(galsim.GaussianNoise(rng=galsim.BaseDeviate(1234), sigma=1e-6))
+
+    # Write out the image to a file
+    image_file = os.path.join('output','trivial_convolve1_im.fits')
+    image.write(image_file)
+
+    # Write out the catalog to a file
+    dtype = [ ('x','f8'), ('y','f8') ]
+    data = np.empty(len(x_list), dtype=dtype)
+    data['x'] = x_list
+    data['y'] = y_list
+    cat_file = os.path.join('output','trivial_convolve1_cat.fits')
+    fitsio.write(cat_file, data, clobber=True)
+
+    psf_file = os.path.join('output','trivial_convolve1_psf.fits')
+    stamp_size = 48
+    config = {
+        'input': {
+            'image_file_name': image_file,
+            'cat_file_name': cat_file,
+            'stamp_size': stamp_size
+        },
+        'psf': {
+            'type': 'Convolve',
+            'components': [
+                {
+                    'type': 'Simple',
+                    'model': { 'type': 'Gaussian',
+                               'fastfit': True,
+                               'include_pixel': False
+                             },
+                    'interp': { 'type': 'Mean' },
+                }
+            ],
+            'max_iter': 10,
+            'chisq_thresh': 0.2,
+        },
+        'output': { 'file_name': psf_file },
+    }
+
+    piff.piffify(config, logger)
+    psf = piff.read(psf_file)
+
+    assert type(psf) is piff.ConvolvePSF
+    assert len(psf.components) == 1
+    assert type(psf.components[0]) is piff.SimplePSF
+    assert type(psf.components[0].model) is piff.Gaussian
+    assert type(psf.components[0].interp) is piff.Mean
+
+    assert psf.chisq_thresh == 0.2
+    assert psf.max_iter == 10
+
+    for i, star in enumerate(psf.stars):
+        target = targets[i]
+        test_star = psf.drawStar(star)
+        true_params = [ sigma, g1, g2 ]
+        np.testing.assert_almost_equal(test_star.fit.params[0], true_params, decimal=4)
+        # The drawn image should be reasonably close to the target.
+        target = target[star.image.bounds]
+        print('target max diff = ',np.max(np.abs(test_star.image.array - target.array)))
+        np.testing.assert_allclose(test_star.image.array, target.array, atol=1.e-5)
+
+        # test that draw works
+        test_image = psf.draw(x=star['x'], y=star['y'], stamp_size=config['input']['stamp_size'],
+                              flux=star.fit.flux, offset=star.fit.center/image.scale)
+        # This image should be very close to the same values as test_star
+        # Just make sure to compare over the same bounds.
+        b = test_star.image.bounds & test_image.bounds
+        print('draw/drawStar max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array)))
+        np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-9)
+
+    # Outlier is a valid option
+    config['psf']['outliers'] = {
+        'type': 'Chisq',
+        'nsigma': 5,
+        'max_remove': 3,
+    }
+    piff.piffify(config)
+    psf2 = piff.read(psf_file)
+    # This didn't remove anything, so the result is the same.
+    # (Use the last star from the above loop for comparison.)
+    test_star2 = psf2.drawStar(star)
+    np.testing.assert_almost_equal(test_star2.fit.params[0], true_params, decimal=4)
+    assert test_star2.image == test_star.image
+
+    test_image2 = psf2.draw(x=star['x'], y=star['y'], stamp_size=config['input']['stamp_size'],
+                           flux=star.fit.flux, offset=star.fit.center/image.scale)
+    assert test_image2 == test_image
+
+    # Error if components is missing
+    config1 = config.copy()
+    del config1['psf']['components']
+    with np.testing.assert_raises(ValueError):
+        piff.process(config)
+
+@timer
+def test_convolve_optatm():
+    # Let reality be atmospheric + optical PSF with realistic variation.
+    # And model this with a convolution of a Kolmogorov/GP SimplePSF and a fixed Optical PSF.
+
+    if __name__ == '__main__':
+        size = 2048
+        nstars = 200
+        noise = 20
+    else:
+        size = 1024
+        nstars = 10
+        noise = 2
+
+    pixel_scale = 0.2
+    im = galsim.ImageF(size, size, scale=pixel_scale)
+
+    screens, aper = make_screens()
+
+    rng = galsim.BaseDeviate(1234)
+    x = rng.np.uniform(25, size-25, size=nstars)
+    y = rng.np.uniform(25, size-25, size=nstars)
+
+    for k in range(nstars):
+        flux = 100000
+        theta = ((x[k] - size/2) * pixel_scale * galsim.arcsec,
+                 (y[k] - size/2) * pixel_scale * galsim.arcsec)
+
+        psf = screens.makePSF(lam=500, aper=aper, exptime=100, flux=flux, theta=theta)
+        psf.drawImage(image=im, center=(x[k],y[k]), method='phot', rng=rng, add_to_image=True)
+        bounds = galsim.BoundsI(int(x[k]-33), int(x[k]+33), int(y[k]-33), int(y[k]+33))
+
+    # Add a little noise
+    noise = 10
+    im.addNoise(galsim.GaussianNoise(rng=rng, sigma=noise))
+    image_file = os.path.join('output', 'convolveatmpsf_im.fits')
+    im.write(image_file)
+
+    # Write out the catalog to a file
+    dtype = [ ('x','f8'), ('y','f8') ]
+    data = np.empty(len(x), dtype=dtype)
+    data['x'] = x
+    data['y'] = y
+    cat_file = os.path.join('output','convolveatmpsf_cat.fits')
+    fitsio.write(cat_file, data, clobber=True)
+
+    psf_file = os.path.join('output','convolveatmpsf.fits')
+    stamp_size = 48
+    config = {
+        'input': {
+            'image_file_name': image_file,
+            'cat_file_name': cat_file,
+            'stamp_size': 32,
+            'noise': noise**2,
+        },
+        'select': {
+            'max_snr': 1.e6,
+            'max_edge_frac': 0.1,
+            'hsm_size_reject': True,
+        },
+        'psf': {
+            'type': 'Convolve',
+            'components': [
+                {
+                    'model': { 'type': 'Kolmogorov',
+                                'fastfit': True,
+                             },
+                    'interp': { 'type': 'GP',
+                              },
+                },
+                {
+                    'model': { 'type': 'Optical',
+                               'atmo_type': 'None',
+                               'lam': 500,
+                               'diam': 8,
+                               # These are the correct aberrations, not fitted.
+                               'base_aberrations': [0,0,0,0,0.7,-0.8,0.7,0.5,0.7,-0.6,0.5,0.8],
+                               'obscuration': 0.4,
+                             },
+                    'interp': { 'type': 'Mean', },
+                }
+            ],
+            'outliers': {
+                'type': 'Chisq',
+                'nsigma': 5,
+                'max_remove': 3,
+            }
+        },
+        'output': {
+            'file_name': psf_file,
+       },
+    }
+
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=1)
+    else:
+        logger = piff.config.setup_logger(log_file='output/test_convolve_optatm.log')
+        logger = piff.config.setup_logger(verbose=1)
+
+    psf = piff.process(config, logger)
+
+    if __name__ == '__main__':
+        config['output']['stats'] = [{
+            'type': 'StarImages',
+            'file_name': os.path.join('output','test_convolve_optatm_stars.png'),
+            'nplot': 10,
+            'adjust_stars': True,
+        }]
+        output = piff.Output.process(config['output'], logger=logger)
+        output.write(psf, logger=logger)
+
+    assert type(psf) is piff.ConvolvePSF
+    assert len(psf.components) == 2
+    assert type(psf.components[0]) is piff.SimplePSF
+    assert type(psf.components[0].model) is piff.Kolmogorov
+    assert type(psf.components[0].interp) is piff.GPInterp
+    assert type(psf.components[1]) is piff.SimplePSF
+    assert type(psf.components[1].model) is piff.Optical
+    assert type(psf.components[1].interp) is piff.Mean
+
+    mean_max_rel_err = 0
+    mean_mse = 0
+    count = 0
+    for i, star in enumerate(psf.stars):
+        if star.is_flagged:
+            continue
+        test_star = psf.drawStar(star)
+
+        b = star.image.bounds.withBorder(-8)
+        max_diff = np.max(np.abs(test_star.image[b].array - star.image[b].array))
+        max_val = np.max(np.abs(star.image[b].array))
+        #print('max_diff / max_val = ',max_diff, max_val, max_diff / max_val)
+        mse = np.sum((test_star.image[b].array - star.image[b].array)**2) / flux**2
+        #print('mse = ',mse)
+        mean_max_rel_err += max_diff/max_val
+        mean_mse += mse
+        count += 1
+
+    mean_max_rel_err /= count
+    mean_mse /= count
+    print('mean maximum relative error = ',mean_max_rel_err)
+    print('mean mean-squared error = ',mean_mse)
+    assert mean_max_rel_err < 0.04
+    assert mean_mse < 1.5e-5
+
+    # Without the Optical component, it's quite a bit worse.
+    config['psf'] = config['psf']['components'][0]
+    psf = piff.process(config, logger)
+
+    if __name__ == '__main__':
+        config['output']['stats'][0]['file_name'] = os.path.join('output', 'test_convolve_optatm_stars2.png')
+        output = piff.Output.process(config['output'], logger=logger)
+        output.write(psf, logger=logger)
+
+    mean_max_rel_err = 0
+    mean_mse = 0
+    count = 0
+    for i, star in enumerate(psf.stars):
+        if star.is_flagged:
+            continue
+        test_star = psf.drawStar(star)
+
+        b = star.image.bounds.withBorder(-8)
+        max_diff = np.max(np.abs(test_star.image[b].array - star.image[b].array))
+        max_val = np.max(np.abs(star.image[b].array))
+        #print('max_diff / max_val = ',max_diff / max_val)
+        mse = np.sum((test_star.image[b].array - star.image[b].array)**2) / flux**2
+        #print('mse = ',mse)
+        mean_max_rel_err += max_diff/max_val
+        mean_mse += mse
+        count += 1
+
+    mean_max_rel_err /= count
+    mean_mse /= count
+    print('mean maximum relative error = ',mean_max_rel_err)
+    print('mean mean-squared error = ',mean_mse)
+    assert mean_max_rel_err > 0.06
+    assert mean_mse > 6.e-5
+
+
+@timer
+def test_convolve_pixelgrid():
+    # Same as test_optatm, but use a PixelGrid for one of the components
+
+    if __name__ == '__main__':
+        size = 2048
+        nstars = 200
+        noise = 20
+    else:
+        size = 1024
+        nstars = 10
+        noise = 2
+
+    pixel_scale = 0.2
+    im = galsim.ImageF(size, size, scale=pixel_scale)
+
+    screens, aper = make_screens()
+
+    rng = galsim.BaseDeviate(1234)
+    x = rng.np.uniform(25, size-25, size=nstars)
+    y = rng.np.uniform(25, size-25, size=nstars)
+
+    for k in range(nstars):
+        flux = 100000
+        theta = ((x[k] - size/2) * pixel_scale * galsim.arcsec,
+                 (y[k] - size/2) * pixel_scale * galsim.arcsec)
+
+        psf = screens.makePSF(lam=500, aper=aper, exptime=100, flux=flux, theta=theta)
+        psf.drawImage(image=im, center=(x[k],y[k]), method='phot', rng=rng, add_to_image=True)
+        bounds = galsim.BoundsI(int(x[k]-33), int(x[k]+33), int(y[k]-33), int(y[k]+33))
+
+    # Add a little noise
+    im.addNoise(galsim.GaussianNoise(rng=rng, sigma=noise))
+    image_file = os.path.join('output', 'convolveatmpsf_im.fits')
+    im.write(image_file)
+
+    # Write out the catalog to a file
+    dtype = [ ('x','f8'), ('y','f8') ]
+    data = np.empty(len(x), dtype=dtype)
+    data['x'] = x
+    data['y'] = y
+    cat_file = os.path.join('output','convolveatmpsf_cat.fits')
+    fitsio.write(cat_file, data, clobber=True)
+
+    psf_file = os.path.join('output','convolveatmpsf.fits')
+    config = {
+        'input': {
+            'image_file_name': image_file,
+            'cat_file_name': cat_file,
+            'stamp_size': 32,
+            'noise': noise**2,
+        },
+        'select': {
+            'max_snr': 1.e6,
+            'max_edge_frac': 0.1,
+            'hsm_size_reject': True,
+        },
+        'psf': {
+            'type': 'Convolve',
+            'max_iter': 5,
+            'components': [
+                {
+                    'model': { 'type': 'PixelGrid',
+                               'scale': pixel_scale,
+                               'size': 17,
+                             },
+                    'interp': { 'type': 'Polynomial',
+                                'order': 1,
+                              },
+                },
+                {
+                    'model': { 'type': 'Optical',
+                               'atmo_type': 'None',
+                               'lam': 500,
+                               'diam': 8,
+                               # These are the correct aberrations, not fitted.
+                               'base_aberrations': [0,0,0,0,0.7,-0.8,0.7,0.5,0.7,-0.6,0.5,0.8],
+                               'obscuration': 0.4,
+                             },
+                    'interp': { 'type': 'Mean', },
+                }
+            ],
+            'outliers': {
+                'type': 'Chisq',
+                'nsigma': 5,
+                'max_remove': 3,
+            }
+        },
+        'output': {
+            'file_name': psf_file,
+       },
+    }
+
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=1)
+    else:
+        logger = piff.config.setup_logger(log_file='output/test_convolve_pixelgrid.log')
+        logger = piff.config.setup_logger(verbose=1)
+
+    psf = piff.process(config, logger)
+
+    if __name__ == '__main__':
+        config['output']['stats'] = [{
+            'type': 'StarImages',
+            'file_name': os.path.join('output','test_convolve_pixelgrid_stars.png'),
+            'nplot': 10,
+            'adjust_stars': True,
+        }]
+        output = piff.Output.process(config['output'], logger=logger)
+        output.write(psf, logger=logger)
+
+    assert type(psf) is piff.ConvolvePSF
+    assert len(psf.components) == 2
+    assert type(psf.components[0]) is piff.SimplePSF
+    assert type(psf.components[0].model) is piff.PixelGrid
+    assert type(psf.components[0].interp) is piff.Polynomial
+    assert type(psf.components[1]) is piff.SimplePSF
+    assert type(psf.components[1].model) is piff.Optical
+    assert type(psf.components[1].interp) is piff.Mean
+
+    mean_max_rel_err = 0
+    mean_mse = 0
+    count = 0
+    for i, star in enumerate(psf.stars):
+        if star.is_flagged:
+            continue
+        test_star = psf.drawStar(star)
+
+        b = star.image.bounds.withBorder(-8)
+        max_diff = np.max(np.abs(test_star.image[b].array - star.image[b].array))
+        max_val = np.max(np.abs(star.image[b].array))
+        #print('max_diff / max_val = ',max_diff, max_val, max_diff / max_val)
+        mse = np.sum((test_star.image[b].array - star.image[b].array)**2) / flux**2
+        #print('mse = ',mse)
+        mean_max_rel_err += max_diff/max_val
+        mean_mse += mse
+        count += 1
+
+    mean_max_rel_err /= count
+    mean_mse /= count
+    print('mean maximum relative error = ',mean_max_rel_err)
+    print('mean mean-squared error = ',mean_mse)
+    assert mean_max_rel_err < 0.04
+    assert mean_mse < 2.e-5
+
+
+@timer
+def test_nested_convolve():
+    # Test that ConvolvePSF can be nested inside another composite PSF (in particular one that
+    # would use a convert_func).
+    # This test is essentially the same as test_convolve_optatm, but with an additional convolution
+    # with a DeltaFunction.
+    # It should work whether this is done the straightforward way with a 3-component Convolve
+    # or with two nested Convolves.
+
+    #if __name__ == '__main__':
+    if False:
+        size = 2048
+        nstars = 200
+        noise = 20
+    else:
+        size = 1024
+        nstars = 10
+        noise = 2
+
+    pixel_scale = 0.2
+    im = galsim.ImageF(size, size, scale=pixel_scale)
+
+    screens, aper = make_screens()
+
+    rng = galsim.BaseDeviate(1234)
+    x = rng.np.uniform(25, size-25, size=nstars)
+    y = rng.np.uniform(25, size-25, size=nstars)
+
+    for k in range(nstars):
+        flux = 100000
+        theta = ((x[k] - size/2) * pixel_scale * galsim.arcsec,
+                 (y[k] - size/2) * pixel_scale * galsim.arcsec)
+
+        psf = screens.makePSF(lam=500, aper=aper, exptime=100, flux=flux, theta=theta)
+        psf.drawImage(image=im, center=(x[k],y[k]), method='phot', rng=rng, add_to_image=True)
+        bounds = galsim.BoundsI(int(x[k]-33), int(x[k]+33), int(y[k]-33), int(y[k]+33))
+
+    # Add a little noise
+    noise = 10
+    im.addNoise(galsim.GaussianNoise(rng=rng, sigma=noise))
+    image_file = os.path.join('output', 'convolve_nested_im.fits')
+    im.write(image_file)
+
+    # Write out the catalog to a file
+    dtype = [ ('x','f8'), ('y','f8') ]
+    data = np.empty(len(x), dtype=dtype)
+    data['x'] = x
+    data['y'] = y
+    cat_file = os.path.join('output','convolve_nested_cat.fits')
+    fitsio.write(cat_file, data, clobber=True)
+
+    # First do it as a three-component convolution.
+    psf_file = os.path.join('output','convolve_nested.fits')
+    stamp_size = 48
+    config = {
+        'input': {
+            'image_file_name': image_file,
+            'cat_file_name': cat_file,
+            'stamp_size': 32,
+            'noise': noise**2,
+        },
+        'select': {
+            'max_snr': 1.e6,
+            'max_edge_frac': 0.1,
+            'hsm_size_reject': True,
+        },
+        'psf': {
+            'type': 'Convolve',
+            'components': [
+                {
+                    'model': { 'type': 'Kolmogorov',
+                               'fastfit': True,
+                             },
+                    'interp': { 'type': 'GP',
+                              },
+                },
+                {
+                    'model': { 'type': 'Optical',
+                               'atmo_type': 'None',
+                               'lam': 500,
+                               'diam': 8,
+                               # These are the correct aberrations, not fitted.
+                               'base_aberrations': [0,0,0,0,0.7,-0.8,0.7,0.5,0.7,-0.6,0.5,0.8],
+                               'obscuration': 0.4,
+                             },
+                    'interp': { 'type': 'Mean', },
+                },
+                {
+                    'model': { 'type': 'Gaussian',
+                               'init': 'delta',
+                             },
+                    'interp': { 'type': 'Mean', },
+                },
+            ],
+            'outliers': {
+                'type': 'Chisq',
+                'nsigma': 5,
+                'max_remove': 3,
+            }
+        },
+        'output': {
+            'file_name': psf_file,
+       },
+    }
+
+    if __name__ == '__main__':
+        logger = piff.config.setup_logger(verbose=1)
+    else:
+        logger = piff.config.setup_logger(log_file='output/test_convolve_optatm.log')
+        logger = piff.config.setup_logger(verbose=1)
+
+    psf = piff.process(config, logger)
+
+    assert type(psf) is piff.ConvolvePSF
+    assert len(psf.components) == 3
+    assert type(psf.components[0]) is piff.SimplePSF
+    assert type(psf.components[0].model) is piff.Kolmogorov
+    assert type(psf.components[0].interp) is piff.GPInterp
+    assert type(psf.components[1]) is piff.SimplePSF
+    assert type(psf.components[1].model) is piff.Optical
+    assert type(psf.components[1].interp) is piff.Mean
+    assert type(psf.components[2]) is piff.SimplePSF
+    assert type(psf.components[2].model) is piff.Gaussian
+    assert type(psf.components[2].interp) is piff.Mean
+
+    mean_max_rel_err = 0
+    mean_mse = 0
+    count = 0
+    for i, star in enumerate(psf.stars):
+        if star.is_flagged:
+            continue
+        test_star = psf.drawStar(star)
+
+        b = star.image.bounds.withBorder(-8)
+        max_diff = np.max(np.abs(test_star.image[b].array - star.image[b].array))
+        max_val = np.max(np.abs(star.image[b].array))
+        mse = np.sum((test_star.image[b].array - star.image[b].array)**2) / flux**2
+        mean_max_rel_err += max_diff/max_val
+        mean_mse += mse
+        count += 1
+
+    mean_max_rel_err /= count
+    mean_mse /= count
+    print('mean maximum relative error = ',mean_max_rel_err)
+    print('mean mean-squared error = ',mean_mse)
+    assert mean_max_rel_err < 0.04
+    assert mean_mse < 1.5e-5
+
+    # Repeat using nested Convolves
+    config['psf'] = {
+        'type': 'Convolve',
+        'components': [
+            {
+                'type': 'Convolve',
+                'components': [
+                    {
+                        'model': { 'type': 'Kolmogorov',
+                                   'fastfit': True,
+                                 },
+                        'interp': { 'type': 'GP',
+                                  },
+                    },
+                    {
+                        'model': { 'type': 'Optical',
+                                   'atmo_type': 'None',
+                                   'lam': 500,
+                                   'diam': 8,
+                                   # These are the correct aberrations, not fitted.
+                                   'base_aberrations': [0,0,0,0,0.7,-0.8,0.7,0.5,0.7,-0.6,0.5,0.8],
+                                   'obscuration': 0.4,
+                                 },
+                        'interp': { 'type': 'Mean', },
+                    },
+                ],
+            },
+            {
+                'model': { 'type': 'Gaussian' },
+                # Note: default init is delta, so don't need to specify it explicitily.
+                'interp': { 'type': 'Mean', },
+            },
+        ],
+    }
+
+    psf = piff.process(config, logger)
+
+    assert type(psf) is piff.ConvolvePSF
+    assert len(psf.components) == 2
+    assert type(psf.components[0]) is piff.ConvolvePSF
+    assert type(psf.components[0].components[0]) is piff.SimplePSF
+    assert type(psf.components[0].components[0].model) is piff.Kolmogorov
+    assert type(psf.components[0].components[0].interp) is piff.GPInterp
+    assert type(psf.components[0].components[1]) is piff.SimplePSF
+    assert type(psf.components[0].components[1].model) is piff.Optical
+    assert type(psf.components[0].components[1].interp) is piff.Mean
+    assert type(psf.components[1]) is piff.SimplePSF
+    assert type(psf.components[1].model) is piff.Gaussian
+    assert type(psf.components[1].interp) is piff.Mean
+
+    mean_max_rel_err = 0
+    mean_mse = 0
+    count = 0
+    for i, star in enumerate(psf.stars):
+        if star.is_flagged:
+            continue
+        test_star = psf.drawStar(star)
+
+        b = star.image.bounds.withBorder(-8)
+        max_diff = np.max(np.abs(test_star.image[b].array - star.image[b].array))
+        max_val = np.max(np.abs(star.image[b].array))
+        mse = np.sum((test_star.image[b].array - star.image[b].array)**2) / flux**2
+        mean_max_rel_err += max_diff/max_val
+        mean_mse += mse
+        count += 1
+
+    mean_max_rel_err /= count
+    mean_mse /= count
+    print('mean maximum relative error = ',mean_max_rel_err)
+    print('mean mean-squared error = ',mean_mse)
+    assert mean_max_rel_err < 0.04
+    assert mean_mse < 1.5e-5
+
+
+
+if __name__ == '__main__':
+    test_trivial_convolve1()
+    test_convolve_optatm()
+    test_convolve_pixelgrid()
+    test_nested_convolve()
diff --git a/tests/test_gsobject_model.py b/tests/test_gsobject_model.py
index e7b7c604..1c05d4b4 100644
--- a/tests/test_gsobject_model.py
+++ b/tests/test_gsobject_model.py
@@ -1091,7 +1091,7 @@ def test_fail():
         stars, _ = psf.initialize_params([star3], logger=cl.logger)
         with np.testing.assert_raises(RuntimeError):
             # Raises an error that all stars were flagged
-            psf.single_iteration(stars, logger=cl.logger, convert_func=None)
+            psf.single_iteration(stars, logger=cl.logger, convert_funcs=None, draw_method=None)
     assert "Failed fitting star" in cl.output
     print('5')
 
diff --git a/tests/test_pixel.py b/tests/test_pixel.py
index e3680382..0c7f3f51 100644
--- a/tests/test_pixel.py
+++ b/tests/test_pixel.py
@@ -1703,6 +1703,54 @@ def test_color():
                                   # Anyway, I think the fit is working, just this test doesn't
                                   # seem quite the right thing.
 
+@timer
+def test_convert_func():
+    """Test PixelGrid fitting with a non-trivial convert_func
+    """
+
+    # This is the kind of convert_func that might be used by ConvolvePSF
+    # Start without noise and where the image is exactly representable by the converted PixelGrid.
+
+    optics = galsim.OpticalPSF(lam=500, diam=8,
+                               aberrations=[0.0,0.0,0.0,0.0,0.7,-0.8,0.7,0.5,0.7,-0.6,0.5,0.8],
+                               obscuration=0.4)
+    gauss = galsim.Gaussian(sigma=0.2)
+    delta = galsim.DeltaFunction()
+
+    convert_funcs = [
+        lambda prof: prof,
+        lambda prof: prof.shear(g1=0.2, g2=-0.05),
+        lambda prof: prof.dilate(1.12),
+        lambda prof: prof.shift(dx=0.1, dy=-0.05),
+        lambda prof: galsim.Convolve(prof, delta),
+        lambda prof: galsim.Convolve(prof, gauss),
+        lambda prof: galsim.Convolve(prof, optics),
+        lambda prof: prof + galsim.Convolve(gauss, optics),
+    ]
+
+    size = 5
+    scale = 0.2
+    flux = 1.e5
+    interp = galsim.Lanczos(7)
+    model_im = galsim.Gaussian(sigma=0.3).drawImage(nx=size,ny=size, scale=scale, method='no_pixel')
+    model_im /= model_im.array.sum()
+    model_ii = galsim.InterpolatedImage(model_im, x_interpolant=interp, use_true_center=False)
+
+    for k, convert_func in enumerate(convert_funcs):
+        print('Test convert_func #',k)
+        prof = convert_func(model_ii)
+
+        star = piff.Star.makeTarget(x=0, y=0, u=0, v=0, scale=0.2, stamp_size=32, flux=flux)
+        prof.withFlux(flux).drawImage(image=star.image, center=star.image_pos, method='no_pixel')
+
+        mod = piff.PixelGrid(scale, size, interp)
+        star1 = mod.initialize(star).withFlux(flux=np.sum(star.image.array))
+        star1 = mod.fit(star1, convert_func=convert_func)
+        print('true params = ',model_im.array.ravel())
+        print('fitted params = ',star1.fit.params)
+        np.testing.assert_allclose(star1.fit.params, model_im.array.ravel(), rtol=2.e-3)
+
+
 if __name__ == '__main__':
     #import cProfile, pstats
     #pr = cProfile.Profile()
@@ -1722,6 +1770,7 @@ def test_color():
     test_des2()
     test_var()
     test_color()
+    test_convert_func()
     #pr.disable()
     #ps = pstats.Stats(pr).sort_stats('tottime')
     #ps.print_stats(20)
diff --git a/tests/test_simple.py b/tests/test_simple.py
index 895e4877..f05418b0 100644
--- a/tests/test_simple.py
+++ b/tests/test_simple.py
@@ -310,10 +310,6 @@ def test_single_image():
     assert test_star_nocopy.image.array[1,1] == target_star_copy.image.array[1,1]
     assert test_star_copy.image.array[1,1] == target_star_copy.image.array[1,1]
 
-    test_star_center = psf.model.draw(test_star_copy, copy_image=True, center=(x+1,y+1))
-    np.testing.assert_almost_equal(test_star_center.image.array[1:,1:],
-                                   test_star_copy.image.array[:-1,:-1])
-
     # test that draw works
     test_image = psf.draw(x=target['x'], y=target['y'], stamp_size=config['input']['stamp_size'],
                           flux=target.fit.flux, offset=target.fit.center)
@@ -780,9 +776,11 @@ def test_psf():
     np.testing.assert_raises(NotImplementedError, psf.drawStarList, [star])
     np.testing.assert_raises(NotImplementedError, psf._drawStar, star)
     np.testing.assert_raises(NotImplementedError, psf._getProfile, star)
-    np.testing.assert_raises(NotImplementedError, psf.single_iteration, [star], None, None)
+    np.testing.assert_raises(NotImplementedError, psf.single_iteration, [star], None, None, None)
     with np.testing.assert_raises(NotImplementedError):
         psf.fit_center
+    with np.testing.assert_raises(NotImplementedError):
+        psf.include_model_centroid
 
     # initialize_params doesn't do anything, but works
     stars, nremove = psf.initialize_params([star], None)
diff --git a/tests/test_wcs.py b/tests/test_wcs.py
index b14c2c08..7755b96b 100644
--- a/tests/test_wcs.py
+++ b/tests/test_wcs.py
@@ -375,6 +375,7 @@ def test_single():
     # 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
+    assert psf.include_model_centroid is False
     raw1 = psf._getRawProfile(star)
     raw2 = psf.psf_by_chip[star['chipnum']]._getRawProfile(star)
     assert raw1 == raw2