Use synphot package to perform bandpass integration

[MetinSa]

from pyOpenSci/software-submission#161:

Instead of manually performing bandpass integrations with the custom Bandpass type, we should use the Astropy affiliated synphot package. This would also provide better integration with Astropy.

[MetinSa]

I noticed that synphot is only available for Python >= 3.10, so I will put this off until ZodiPy drop 3.9 unless there is a work around @pllim? Is it correct that according to https://numpy.org/neps/nep-0029-deprecation_policy.html astropy-affiliated packages should drop python 3.9 already on Apr 05 this year?

[pllim]

Hello and I apologize for any inconvenience. Yes, synphot dropped support for Python 3.9 but only recently via spacetelescope/synphot_refactor#387 . You should still be able to get synphot < 1.4 in Python 3.9 .

[pllim]

Would be happy to discuss more if you are interested. Integration is nothing too fancy. synphot uses trapezoid integration under the curve provided by scipy (used to be np.trapz in older versions). See https://synphot.readthedocs.io/en/latest/api/synphot.spectrum.BaseSpectrum.html#synphot.spectrum.BaseSpectrum.integrate (analytical was added much later and defers to astropy.modeling but not widely used).

[MetinSa]

Thanks it should work with synphot < 1.4.

Yeah that would be nice! I want to use synphot to compute a table of bandpass integrated blackbody emission values where the blackbodies have varying temperature. It should work for either a center frequency or an empirical bandpass provided by the user as two 1D arrays corresponding to weights and frequencies. See get_bandpass_interpolation_table in https://github.com/Cosmoglobe/zodipy/blob/main/zodipy/_bandpass.py.
Not quite sure how I would do this using these models but synphot.models.Empirical1d or synphot.spectrum.BaseSpectrum seem relevant? I also see that there is a synphot.models.Blackbody1D which means that I probably don't need to use any custom b_nu function in ZodiPy either?

[pllim]

Generally speaking, this would give you a blackbody spectrum (before passing through the bandpass):

from astropy import units as u
from synphot import SourceSpectrum
from synphot.models import BlackBodyNorm1D

bb_temp = 5777 * u.K
bb = SourceSpectrum(BlackBodyNorm1D, temperature=bb_temp)

# Sampling frequencies. This samples most of the flux but feel free to extend as needed.
freq = bb.waveset.to(u.Hz, u.spectral())

# Sampled flux
flux_unit = u.MJy
flux = bb(freq, flux_unit=flux_unit)

But to put it through bandpass first before sampling, you have to construct an Observation, as documented in https://synphot.readthedocs.io/en/latest/synphot/observation.html .

Often times, you might also want to renormalize the source spectrum first to some standard you want before passing it through the bandpass using https://synphot.readthedocs.io/en/latest/api/synphot.spectrum.BaseSourceSpectrum.html#synphot.spectrum.BaseSourceSpectrum.normalize .

So a more complete example would be more like this:

from astropy import units as u
from synphot import SourceSpectrum, SpectralElement, Observation, units as syn_units
from synphot.models import BlackBodyNorm1D, Box1D

bb_temp = 5777 * u.K
bb = SourceSpectrum(BlackBodyNorm1D, temperature=bb_temp)

# Say, you know it is 17 VEGAMAG in Johnson V filter
vega = SourceSpectrum.from_vega()  # For unit conversion  
johnson_v = SpectralElement.from_filter('johnson_v')  
bb_norm = bb.normalize(17 * syn_units.VEGAMAG, johnson_v, vegaspec=vega)

# Now pass the normalized blackbody into your bandpass. Using box here but there are many available.
bandpass = SpectralElement(Box1D, amplitude=1, x_0=5000 * u.AA, width=10 * u.AA)
obs = Observation(bb_norm, bandpass)

# https://synphot.readthedocs.io/en/latest/synphot/formulae.html#synphot-formula-effstim
effstim = obs.effstim(flux_unit=u.MJy)

When it comes to an Observation, usually people want the effective stimulus instead of integration under the curve.

Hope this helps!

[MetinSa]

Thanks you, that helps!

I am still having some trouble... Is this the way to go if I want to bandpass integrate for a range of different blackbodies:

bandpass = SpectralElement(Empirical1D, points=frequencies, lookup_table=weights)
for temp in temperatures:
    bb = SourceSpectrum(BlackBody1D(temp))
    obs = Observation(bb, bandpass)
    bp_integrated_emission = obs.integrate()

I am also interested in getting the emission in spectral flux density units, so for instance MJy/sr.
How efficient is this compared to just evaluation b_nu(frequencies) and then bandpass integrating with np.trapz for instance?

[pllim]

@MetinSa , is that a real bandpass or you are using it to represent something else? If you have a minimal example with just numbers and expected answers, I can see what in synphot can best represent your use case. Thanks!

[MetinSa]

@pllim yes the frequency and weights arrays represent a user inputted bandpass in normalized Jy units. I made a gist where i read in the DIRBE 25um bandpass and do the bandpass integration i wish to do in ZodiPy with synphot https://gist.github.com/MetinSa/9732a366fbe932a044819dcaa948df88.

Thank you for being so quick and nice to respond:)

[pllim]

When I run your code, I find that normalized_weights (the actual numbers multiplied to blackbody flux) is all negative. In the context of "bandpass", in which the values should only be between 0 and 1 (inclusive), this is un-physical. Could you please clarify?

[MetinSa]

Ahh, I had forgotton to flip the frequencies and weights arrays after changing to frequency convention so that its in increasing order. Updated the gist now.

[pllim]

Just for my own record, looks like you are only interested in this part of the blackbody curve.

[MetinSa]

Right, for this particular bandpass, but it should be able to handle arbitrary frequency/wavelength ranges and temperatures between ~50-500K

[pllim]

I don't grok the bandpass normalization, but that doesn't matter.

I tried re-constructing what I think you are doing with synphot API but I cannot get the numbers to match exactly. If all you want to do is simple trapezoid integration with some custom normalization, synphot API might be a little of an overkill, because there are a lot of assumptions and checks that went into full-blown synthetic photometry. Furthermore, for historical reason, all the internal calculations are done in wavelength space but looks like you want to calculate in frequency space natively.

However, synphot does have a synphot.blackbody.blackbody_nu function that you can use to replace your b_nu. p.s. I noticed that your b_nu is missing the per-steradian in returned flux unit.

Also, in tracing your calculation, I think you need to get rid of the per-micron in the integrated flux unit (because you already integrated over a range of Hz/micron), so the answer I think should be 2.2052944808204042 W / m2 and not 2.2052944808204042 W / (um m2).

I hope this helps and thanks for considering synphot.

[pllim]

Also, in numpy 2.0, np.trapz has been deprecated in favor of scipy.integrate.trapezoid. FYI.

[MetinSa]

I had the feeling that synphot might be a bit overkill for this purpose yes, but this was useful for me nevertheless. I will simply move over to use astropy.models.physical_model.BlackBody which can handle wavelengths too , instead of the custom b_nu implementation.

Also, in numpy 2.0, np.trapz has been deprecated in favor of scipy.integrate.trapezoid. FYI.

Thats a good catch, thanks @pllim!

[pllim]

Yes, the blackbody model might work too. Good luck!