Donut Models

Models for fitting single or multiple out-of-focus (donut) stellar images.

Base class

danish.donut_model.BaseDonutModel

Base class for donut models.

Parameters:

Name	Type	Description	Default
factory	DonutFactoryBase		required
bkg_order	int	Order of polynomial background model to use. If -1, no background.	required
npix	int	Number of pixels along image edge. Must be odd.	required
seed	int	Random seed for use when creating noisy donut images.	required
atm_mode	str	Atmospheric PSF parameterization. 'fwhm' for scalar FWHM (default), 'ixx' for second moment tensor (Ixx, Ixy, Iyy) in arcsec^2.	'fwhm'

Source code in danish/donut_model.py

class BaseDonutModel:
    """Base class for donut models.

    Parameters
    ----------
    factory : DonutFactoryBase
    bkg_order : int
        Order of polynomial background model to use.  If -1, no background.
    npix : int
        Number of pixels along image edge. Must be odd.
    seed : int
        Random seed for use when creating noisy donut images.
    atm_mode : str, optional
        Atmospheric PSF parameterization.  'fwhm' for scalar FWHM (default),
        'ixx' for second moment tensor (Ixx, Ixy, Iyy) in arcsec^2.
    """
    # I think Ixx might formally diverge for a Kolmogorov profile, but we need
    # something to convert between FWHM and Ixx/Iyy for the atmospheric kernel,
    # So we will just use the constant factor here.  This is computed from
    # galsim.Kolmogorov(fwhm=1).calculateMomentRadius().  So we have:
    #   sigma = KOLM_CONST * fwhm
    #   Ixx = Iyy = sigma^2 = (KOLM_CONST * fwhm)^2
    # and
    #   fwhm = sigma / KOLM_CONST = sqrt(Ixx) / KOLM_CONST
    # which we will generalize to
    #   fwhm = sqrt((Ixx + Iyy) / 2) / KOLM_CONST

    KOLM_CONST = 0.651

    def __init__(
        self, factory, npix, seed, bkg_order, atm_mode='fwhm', loss_fn=None,
    ):
        assert npix % 2 == 1, "npix must be odd"
        assert atm_mode in ('fwhm', 'ixx'), "atm_mode must be 'fwhm' or 'ixx'"

        self.factory = factory
        self.npix = npix
        self.no2 = (npix-1)//2
        self.gsrng = galsim.BaseDeviate(seed)
        # arcseconds per pixel
        # This is only used for the atmospheric part of the model, where we
        # ignore distortion and draw a circular profile in image coordinates.
        self.sky_scale = (
            3600*np.rad2deg(1/factory.focal_length)*factory.pixel_scale
        )
        self.bkg_order = bkg_order
        self.nbkg = (bkg_order+1)*(bkg_order+2)//2
        self.atm_mode = atm_mode
        self.loss_fn = loss_fn if loss_fn is not None else chi2_loss

    @property
    def natm(self):
        """Number of atmospheric parameters (1 for fwhm, 3 for ixx)."""
        return 1 if self.atm_mode == 'fwhm' else 3

    @lru_cache(maxsize=1000)
    def _atm_kernel(self, dx, dy, Ixx, Ixy, Iyy):
        """Compute atmospheric kernel.

        Parameters
        ----------
        dx, dy : float
            Offset in arcseconds.
        Ixx, Ixy, Iyy : float, optional
            Second moment tensor in arcsec^2.

        Returns
        -------
        array of float
            Atmospheric kernel array.
        """
        fwhm = np.sqrt((Ixx + Iyy) / 2) / self.KOLM_CONST
        e1 = (Ixx - Iyy) / (Ixx + Iyy)
        e2 = 2*Ixy / (Ixx + Iyy)
        obj = galsim.Kolmogorov(fwhm=fwhm).shear(e1=e1, e2=e2).shift(dx, dy)
        img = obj.drawImage(nx=self.no2, ny=self.no2, scale=self.sky_scale)
        return img.array

    @lru_cache(maxsize=1000)
    def _opt_kernel(
        self,
        zk,
        thx, thy,
        x_offset=None, y_offset=None,
    ):
        """Compute optical kernel.

        Parameters
        ----------
        zk : array of float
            Zernike coefficients.
        thx, thy : float
            Field angle in radians.
        x_offset, y_offset : galsim.zernike.Zernike, optional
            Additional pupil distortion coefficients.

        Returns
        -------
        array of float
            Optical kernel array.
        """
        result = self.factory.image(
            aberrations=zk,
            x_offset=x_offset, y_offset=y_offset,
            thx=thx, thy=thy,
            npix=self.npix,
        )
        return result

    @lru_cache(maxsize=1000)
    def _bkg(
        self,
        bkg
    ):
        """Compute background model.

        Parameters
        ----------
        bkg : array of float
            Background polynomial coefficients.

        Returns
        -------
        array of float
            Background model array.
        """
        no2 = self.no2
        zkbkg = galsim.zernike.Zernike([0]+list(bkg), R_outer=no2)
        result = zkbkg(*np.mgrid[-no2:no2+1, -no2:no2+1][::-1])
        return result

    @lru_cache(maxsize=100)
    def _model(
        self,
        flux, dx, dy, Ixx, Ixy, Iyy, zk, bkg,
        thx, thy,
        x_offset=None, y_offset=None,
        sky_level=None
    ):
        """Compute model image.

        Parameters
        ----------
        flux : float
            Flux level at which to set image.
        dx, dy : float
            Offset in arcseconds.
        Ixx, Ixy, Iyy : float
            Second moment tensor in arcsec^2.
        zk : array of float
            Zernike coefficients.
        bkg : array of float
            Background polynomial coefficients.
        thx, thy : float
            Field angle in radians.
        x_offset, y_offset : galsim.zernike.Zernike, optional
            Additional pupil distortion coefficients.
        sky_level : float, optional
            Sky level to use when adding Poisson noise to image.

        Returns
        -------
        array of float
            Model image array.
        """
        atm = self._atm_kernel(dx, dy, Ixx, Ixy, Iyy)
        opt = self._opt_kernel(zk, thx, thy, x_offset=x_offset, y_offset=y_offset)
        arr = convolve(opt, atm)
        arr *= flux

        if bkg:
            arr += self._bkg(bkg)

        if sky_level is not None:
            pn = galsim.PoissonNoise(self.gsrng, sky_level=sky_level)
            img = galsim.Image(arr)
            img.addNoise(pn)
            arr = img.array

        return arr

    def _chi(
        self,
        data,
        model,
        var
    ):
        """Compute chi = (data - model)/error.

        Parameters
        ----------
        data : array of float
            Observed image data.
        model : array of float
            Model image data.
        var : array of float
            Variance of the observed data.

        Returns
        -------
        array of float
            Chi values for each pixel.
        """
        result = self.loss_fn(data, model, var).ravel()
        return result

natm property

natm

Number of atmospheric parameters (1 for fwhm, 3 for ixx).

Single-donut model

danish.donut_model.SingleDonutModel

Bases: BaseDonutModel

Model individual donuts using single Zernike offsets to a reference single Zernike series.

Parameters:

Name	Type	Description	Default
factory	DonutFactoryBase		required
bkg_order	int	Order of polynomial background model to use. If -1, no background.	-1
z_ref	array of float	Constant reference Zernike coefs to add to fitted coefficients. [0] is ignored, [1] is piston, [4] is defocus, etc.	None
x_offset	Zernike	Additional pupil distortion coefficients.	None
y_offset	Zernike	Additional pupil distortion coefficients.	None
z_terms	sequence of int	Which Zernike coefficients to include in the fit. E.g., [4,5,6,11] will fit defocus, astigmatism, and spherical.	()
thx	float	Field angle in radians.	None
thy	float	Field angle in radians.	None
npix	int	Number of pixels along image edge. Must be odd.	181
seed	int	Random seed for use when creating noisy donut images with this class.	57721

Source code in danish/donut_model.py

class SingleDonutModel(BaseDonutModel):
    """Model individual donuts using single Zernike offsets to a reference
    single Zernike series.

    Parameters
    ----------
    factory : DonutFactoryBase
    bkg_order : int, optional
        Order of polynomial background model to use.  If -1, no background.
    z_ref : array of float
        Constant reference Zernike coefs to add to fitted coefficients.
        [0] is ignored, [1] is piston, [4] is defocus, etc.
    x_offset, y_offset : galsim.zernike.Zernike, optional
        Additional pupil distortion coefficients.
    z_terms : sequence of int
        Which Zernike coefficients to include in the fit.  E.g.,
        [4,5,6,11] will fit defocus, astigmatism, and spherical.
    thx, thy : float
        Field angle in radians.
    npix : int
        Number of pixels along image edge.  Must be odd.
    seed : int
        Random seed for use when creating noisy donut images with this class.
    """
    def __init__(
        self,
        factory, *,
        bkg_order=-1,
        z_ref=None,
        x_offset=None, y_offset=None,
        z_terms=(),
        thx=None, thy=None,
        npix=181,
        seed=57721,
        loss_fn=None,
    ):
        super().__init__(
            factory, npix, seed, bkg_order, loss_fn=loss_fn,
        )
        self.z_ref = z_ref
        self.x_offset = x_offset
        self.y_offset = y_offset
        self.z_terms = z_terms
        self.thx = thx
        self.thy = thy

    def model(
        self,
        flux, dx, dy, fwhm, z_fit, *,
        bkg=(),
        sky_level=None,
    ):
        """Compute donut model image.

        Parameters
        ----------
        flux : float
            Flux level at which to set image.
        dx, dy : float
            Offset in arcseconds.
        fwhm : float
            Full width half maximum of Kolmogorov kernel.
        z_fit : sequence of float
            Zernike perturbations.
        bkg : sequence of float
            Background polynomial coefficients.
        sky_level : float
            Sky level to use when adding Poisson noise to image.

        Returns
        -------
        img : array of float
            Model image.
        """
        Ixx = Iyy = (fwhm * self.KOLM_CONST) ** 2
        Ixy = 0.0
        zk = np.array(self.z_ref)
        for i, term in enumerate(self.z_terms):
            zk[term] += z_fit[i]
        return self._model(
            flux, dx, dy, Ixx, Ixy, Iyy,
            tuple(zk), tuple(bkg),
            self.thx, self.thy,
            x_offset=self.x_offset, y_offset=self.y_offset,
            sky_level=sky_level
        )

    def pack_params(self, *, flux, dx, dy, fwhm, z_fit, bkg=()):
        """Pack parameters into a single tuple for optimization.

        Parameters
        ----------
        flux : float
            Flux level at which to set image.
        dx, dy : float
            Offset in arcseconds.
        fwhm : float
            Full width half maximum of Kolmogorov kernel.
        z_fit : array of float
            Zernike perturbations.
        bkg : array of float
            Background polynomial coefficients.

        Returns
        -------
        tuple of float
            Packed parameters.
        """
        return (flux, dx, dy, fwhm, *z_fit, *bkg)

    def unpack_params(self, params):
        """
        Unpack parameters from the single tuple used for optimization.

        Parameters
        ----------
        params : sequence of float
            Packed parameters tuple as returned by `pack_params`.

        Returns
        -------
        dict
            Dictionary with keys 'flux', 'dx', 'dy', 'fwhm', 'z_fit', and 'bkg'.
        """
        flux, dx, dy, fwhm, *rest = params
        out = dict(
            flux=flux,
            dx=dx,
            dy=dy,
            fwhm=fwhm,
        )
        split = len(rest) - self.nbkg
        out["z_fit"] = rest[:split]
        out["bkg"] = rest[split:]
        return out

    def chi(
        self, params, data, var
    ):
        """Compute chi = (data - model)/error.

        The error is modeled as sqrt(model + var).

        Parameters
        ----------
        params : sequence of float
            Order is: (flux, dx, dy, fwhm, *z_fit, *bkg)
        data : array of float (npix, npix)
            Image against which to compute chi.
        var : float or array of float (npix, npix)
            Variance of the sky only.  Do not include Poisson contribution of
            the signal, as this will be added from the current model.

        Returns
        -------
        chi : array of float
            Flattened array of chi residuals.
        """
        mod = self.model(**self.unpack_params(params))
        return self._chi(data, mod, var)

    def jac(
        self, params, data, var
    ):
        """Compute jacobian d(chi)/d(param).

        Parameters
        ----------
        params : sequence of float
            Order is: (flux, dx, dy, fwhm, *z_fit, *bkg)
        data : array of float
            Image against which to compute chi.
        var : float or array of float (npix, npix)
            Variance of the sky only.  Do not include Poisson contribution of
            the signal, as this will be added from the current model.

        Returns
        -------
        jac : array of float
            Jacobian array d(chi)/d(param).  First dimension is pixels, second
            dimension is param.
        """
        out = np.zeros((self.npix**2, len(params)))
        chi0 = self.chi(params, data, var)

        flux = params[0]
        step = [
            1e-3 * flux,  # flux
            0.01,         # dx
            0.01,         # dy
            0.01          # fwhm
        ]
        step += [1e-8]*len(self.z_terms)  # Zernikes (in meters)
        step += [0.01]*self.nbkg

        for i in range(len(params)):
            params1 = np.array(params)
            params1[i] += step[i]
            chi1 = self.chi(params1, data, var)
            out[:, i] = (chi1-chi0)/step[i]

        return out

model

model(flux, dx, dy, fwhm, z_fit, *, bkg=(), sky_level=None)

Compute donut model image.

Parameters:

Name	Type	Description	Default
flux	float	Flux level at which to set image.	required
dx	float	Offset in arcseconds.	required
dy	float	Offset in arcseconds.	required
fwhm	float	Full width half maximum of Kolmogorov kernel.	required
z_fit	sequence of float	Zernike perturbations.	required
bkg	sequence of float	Background polynomial coefficients.	()
sky_level	float	Sky level to use when adding Poisson noise to image.	None

Returns:

Name	Type	Description
img	array of float	Model image.

Source code in danish/donut_model.py

def model(
    self,
    flux, dx, dy, fwhm, z_fit, *,
    bkg=(),
    sky_level=None,
):
    """Compute donut model image.

    Parameters
    ----------
    flux : float
        Flux level at which to set image.
    dx, dy : float
        Offset in arcseconds.
    fwhm : float
        Full width half maximum of Kolmogorov kernel.
    z_fit : sequence of float
        Zernike perturbations.
    bkg : sequence of float
        Background polynomial coefficients.
    sky_level : float
        Sky level to use when adding Poisson noise to image.

    Returns
    -------
    img : array of float
        Model image.
    """
    Ixx = Iyy = (fwhm * self.KOLM_CONST) ** 2
    Ixy = 0.0
    zk = np.array(self.z_ref)
    for i, term in enumerate(self.z_terms):
        zk[term] += z_fit[i]
    return self._model(
        flux, dx, dy, Ixx, Ixy, Iyy,
        tuple(zk), tuple(bkg),
        self.thx, self.thy,
        x_offset=self.x_offset, y_offset=self.y_offset,
        sky_level=sky_level
    )

pack_params

pack_params(*, flux, dx, dy, fwhm, z_fit, bkg=())

Pack parameters into a single tuple for optimization.

Parameters:

Name	Type	Description	Default
flux	float	Flux level at which to set image.	required
dx	float	Offset in arcseconds.	required
dy	float	Offset in arcseconds.	required
fwhm	float	Full width half maximum of Kolmogorov kernel.	required
z_fit	array of float	Zernike perturbations.	required
bkg	array of float	Background polynomial coefficients.	()

Returns:

Type	Description
tuple of float	Packed parameters.

Source code in danish/donut_model.py

def pack_params(self, *, flux, dx, dy, fwhm, z_fit, bkg=()):
    """Pack parameters into a single tuple for optimization.

    Parameters
    ----------
    flux : float
        Flux level at which to set image.
    dx, dy : float
        Offset in arcseconds.
    fwhm : float
        Full width half maximum of Kolmogorov kernel.
    z_fit : array of float
        Zernike perturbations.
    bkg : array of float
        Background polynomial coefficients.

    Returns
    -------
    tuple of float
        Packed parameters.
    """
    return (flux, dx, dy, fwhm, *z_fit, *bkg)

unpack_params

unpack_params(params)

Unpack parameters from the single tuple used for optimization.

Parameters:

Name	Type	Description	Default
params	sequence of float	Packed parameters tuple as returned by pack_params.	required

Returns:

Type	Description
dict	Dictionary with keys 'flux', 'dx', 'dy', 'fwhm', 'z_fit', and 'bkg'.

Source code in danish/donut_model.py

def unpack_params(self, params):
    """
    Unpack parameters from the single tuple used for optimization.

    Parameters
    ----------
    params : sequence of float
        Packed parameters tuple as returned by `pack_params`.

    Returns
    -------
    dict
        Dictionary with keys 'flux', 'dx', 'dy', 'fwhm', 'z_fit', and 'bkg'.
    """
    flux, dx, dy, fwhm, *rest = params
    out = dict(
        flux=flux,
        dx=dx,
        dy=dy,
        fwhm=fwhm,
    )
    split = len(rest) - self.nbkg
    out["z_fit"] = rest[:split]
    out["bkg"] = rest[split:]
    return out

chi

chi(params, data, var)

Compute chi = (data - model)/error.

The error is modeled as sqrt(model + var).

Parameters:

Name	Type	Description	Default
params	sequence of float	Order is: (flux, dx, dy, fwhm, z_fit, bkg)	required
data	array of float (npix, npix)	Image against which to compute chi.	required
var	float or array of float (npix, npix)	Variance of the sky only. Do not include Poisson contribution of the signal, as this will be added from the current model.	required

Returns:

Name	Type	Description
chi	array of float	Flattened array of chi residuals.

Source code in danish/donut_model.py

def chi(
    self, params, data, var
):
    """Compute chi = (data - model)/error.

    The error is modeled as sqrt(model + var).

    Parameters
    ----------
    params : sequence of float
        Order is: (flux, dx, dy, fwhm, *z_fit, *bkg)
    data : array of float (npix, npix)
        Image against which to compute chi.
    var : float or array of float (npix, npix)
        Variance of the sky only.  Do not include Poisson contribution of
        the signal, as this will be added from the current model.

    Returns
    -------
    chi : array of float
        Flattened array of chi residuals.
    """
    mod = self.model(**self.unpack_params(params))
    return self._chi(data, mod, var)

jac

jac(params, data, var)

Compute jacobian d(chi)/d(param).

Parameters:

Name	Type	Description	Default
params	sequence of float	Order is: (flux, dx, dy, fwhm, z_fit, bkg)	required
data	array of float	Image against which to compute chi.	required
var	float or array of float (npix, npix)	Variance of the sky only. Do not include Poisson contribution of the signal, as this will be added from the current model.	required

Returns:

Name	Type	Description
jac	array of float	Jacobian array d(chi)/d(param). First dimension is pixels, second dimension is param.

Source code in danish/donut_model.py

def jac(
    self, params, data, var
):
    """Compute jacobian d(chi)/d(param).

    Parameters
    ----------
    params : sequence of float
        Order is: (flux, dx, dy, fwhm, *z_fit, *bkg)
    data : array of float
        Image against which to compute chi.
    var : float or array of float (npix, npix)
        Variance of the sky only.  Do not include Poisson contribution of
        the signal, as this will be added from the current model.

    Returns
    -------
    jac : array of float
        Jacobian array d(chi)/d(param).  First dimension is pixels, second
        dimension is param.
    """
    out = np.zeros((self.npix**2, len(params)))
    chi0 = self.chi(params, data, var)

    flux = params[0]
    step = [
        1e-3 * flux,  # flux
        0.01,         # dx
        0.01,         # dy
        0.01          # fwhm
    ]
    step += [1e-8]*len(self.z_terms)  # Zernikes (in meters)
    step += [0.01]*self.nbkg

    for i in range(len(params)):
        params1 = np.array(params)
        params1[i] += step[i]
        chi1 = self.chi(params1, data, var)
        out[:, i] = (chi1-chi0)/step[i]

    return out

Multi-donut models

danish.donut_model.DZMultiDonutModel

Bases: BaseMultiDonutModel

Multi donut model that uses double Zernike coefficients to parameterize the wavefront.

Parameters:

Name	Type	Description	Default
factory	DonutFactoryBase		required
dz_terms	sequence of tuple of int	List of (k, j) indices specifying which double Zernike terms to use.	()
bkg_order	int	Order of the background polynomial to fit. If -1, no background.	required
dz_ref	DoubleZernike	Double Zernike coefficients to use for constructing Single Zernike reference coefficients to use for each modeled donut. Either this kwarg or z_refs must be set.	required
z_refs	array of float	Single Zernike reference coefficients for each donut. First dimension is donut, second dimension is pupil Zernike coefficient.	required
field_radius	float	Field radius in radians. If dz_ref is provided, then this is ignored and the field radius will be inferred from dz_ref.	required
thxs	float	Field angles in radians.	required
thys	float	Field angles in radians.	required
npix	int	Number of pixels along image edge. Must be odd.	required
seed	int	Random seed for use when creating noisy donut images with this class.	required

Source code in danish/donut_model.py

class DZMultiDonutModel(BaseMultiDonutModel):
    """Multi donut model that uses double Zernike coefficients to parameterize the
    wavefront.

    Parameters
    ----------
    factory : DonutFactoryBase
    dz_terms : sequence of tuple of int
        List of (k, j) indices specifying which double Zernike terms to use.
    bkg_order : int, optional
        Order of the background polynomial to fit.  If -1, no background.
    dz_ref : DoubleZernike
        Double Zernike coefficients to use for constructing Single Zernike
        reference coefficients to use for each modeled donut.  Either this kwarg
        or `z_refs` must be set.
    z_refs : array of float
        Single Zernike reference coefficients for each donut.  First dimension
        is donut, second dimension is pupil Zernike coefficient.
    field_radius : float
        Field radius in radians.  If dz_ref is provided, then this is ignored and
        the field radius will be inferred from dz_ref.
    thxs, thys : float
        Field angles in radians.
    npix : int
        Number of pixels along image edge.  Must be odd.
    seed : int
        Random seed for use when creating noisy donut images with this class.
    """
    def __init__(self, *args, dz_terms=(), **kwargs):
        # Want dz_terms
        self.dz_terms = dz_terms
        self.k_max = max([term[0] for term in dz_terms], default=0)
        self.j_max = max([term[1] for term in dz_terms], default=0)
        super().__init__(*args, **kwargs)

    def _get_z_fits(self, dz_fit):
        """Convert this class's wavefront parameterization into individual donut
        Zernike coefficients.

        Parameters
        ----------
        dz_fit : array of float
            Array of double Zernike coefficients for this model.

        Returns
        -------
        z_fits : list of array of float
            List of single Zernike coefficients for each donut.
        """
        arr = np.zeros((self.k_max+1, self.j_max+1))
        for i, (k, j) in enumerate(self.dz_terms):
            arr[k, j] = dz_fit[i]
        dz = galsim.zernike.DoubleZernike(arr, uv_outer=self.field_radius)
        z_fits = []
        for thx, thy in zip(self.thxs, self.thys):
            z_fits.append(dz.xycoef(thx, thy))
        return z_fits

danish.donut_model.DZBasisMultiDonutModel

Bases: BaseMultiDonutModel

Multi donut model that uses a sensitivity matrix to convert mode coefficients into double Zernike coefficients to parameterize the wavefront.

Parameters:

Name	Type	Description	Default
factory	DonutFactoryBase		required
sensitivity	array of float	Sensitivity matrix that converts mode coefficients into double Zernike coefficients. Dimensions are (nmode, k_max+1, j_max+1)	None
bkg_order	int	Order of the background polynomial to fit. If -1, no background.	required
dz_ref	DoubleZernike	Double Zernike coefficients to use for constructing Single Zernike reference coefficients to use for each modeled donut. Either this kwarg or z_refs must be set.	required
z_refs	array of float	Single Zernike reference coefficients for each donut. First dimension is donut, second dimension is pupil Zernike coefficient.	required
field_radius	float	Field radius in radians. If dz_ref is provided, then this is ignored and the field radius will be inferred from dz_ref.	required
thxs	float	Field angles in radians.	required
thys	float	Field angles in radians.	required
npix	int	Number of pixels along image edge. Must be odd.	required
seed	int	Random seed for use when creating noisy donut images with this class.	required

Source code in danish/donut_model.py

class DZBasisMultiDonutModel(BaseMultiDonutModel):
    """Multi donut model that uses a sensitivity matrix to convert mode coefficients
    into double Zernike coefficients to parameterize the wavefront.

    Parameters
    ----------
    factory : DonutFactoryBase
    sensitivity : array of float
        Sensitivity matrix that converts mode coefficients into double Zernike
        coefficients.  Dimensions are (nmode, k_max+1, j_max+1)
    bkg_order : int, optional
        Order of the background polynomial to fit.  If -1, no background.
    dz_ref : DoubleZernike
        Double Zernike coefficients to use for constructing Single Zernike
        reference coefficients to use for each modeled donut.  Either this kwarg
        or `z_refs` must be set.
    z_refs : array of float
        Single Zernike reference coefficients for each donut.  First dimension
        is donut, second dimension is pupil Zernike coefficient.
    field_radius : float
        Field radius in radians.  If dz_ref is provided, then this is ignored and
        the field radius will be inferred from dz_ref.
    thxs, thys : float
        Field angles in radians.
    npix : int
        Number of pixels along image edge.  Must be odd.
    seed : int
        Random seed for use when creating noisy donut images with this class.
    """
    def __init__(self, *args, sensitivity=None, **kwargs):
        if sensitivity is None:
            raise ValueError("Must provide sensitivity")
        self.sensitivity = sensitivity
        self.nmode = self.sensitivity.shape[0]
        self.k_max = self.sensitivity.shape[1]
        self.j_max = self.sensitivity.shape[2]
        super().__init__(*args, **kwargs)

    def _get_z_fits(self, mode_coefs):
        """Convert this class's wavefront parameterization into individual donut
        Zernike coefficients.

        Parameters
        ----------
        mode_coefs : array of float
            Array of mode coefficients for this model.

        Returns
        -------
        z_fits : list of array of float
            List of single Zernike coefficients for each donut.
        """
        arr = np.einsum(
            "i,ikj->kj",
            mode_coefs,
            self.sensitivity
        )
        dz = galsim.zernike.DoubleZernike(arr, uv_outer=self.field_radius)
        z_fits = []
        for thx, thy in zip(self.thxs, self.thys):
            z_fits.append(dz.xycoef(thx, thy))
        return z_fits