`wlen` and `flam` give a tabulated model spectrum in wavelength and f_λ
        units. We interpolate linearly over both the model and the bandpass
        since they're both discretely sampled.

        Note that quadratic interpolation is both much slower and can blow up
        fatally in some cases. The latter issue might have to do with really large
        X values that aren't zero-centered, maybe?

        I used to use the quadrature integrator, but Romberg doesn't issue
        complaints the way quadrature did. I should probably acquire some idea
        about what's going on under the hood.

def synphot(self, wlen, flam):
        """`wlen` and `flam` give a tabulated model spectrum in wavelength and f_λ
        units. We interpolate linearly over both the model and the bandpass
        since they're both discretely sampled.

        Note that quadratic interpolation is both much slower and can blow up
        fatally in some cases. The latter issue might have to do with really large
        X values that aren't zero-centered, maybe?

        I used to use the quadrature integrator, but Romberg doesn't issue
        complaints the way quadrature did. I should probably acquire some idea
        about what's going on under the hood.

        """
        from scipy.interpolate import interp1d
        from scipy.integrate import romberg

        d = self._ensure_data()

        mflam = interp1d(wlen, flam,
                         kind='linear',
                         bounds_error=False, fill_value=0)

        mresp = interp1d(d.wlen, d.resp,
                         kind='linear',
                         bounds_error=False, fill_value=0)

        bmin = d.wlen.min()
        bmax = d.wlen.max()

        numer = romberg(lambda x: mresp(x) * mflam(x),
                        bmin, bmax, divmax=20)
        denom = romberg(lambda x: mresp(x),
                        bmin, bmax, divmax=20)
        return numer / denom