code: def synphot(self, wlen, flam):
          
          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
docs: `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.