Closed MetinSa closed 4 months ago
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?
Hello and I apologize for any inconvenience. Yes, synphot dropped support for Python 3.9 but only recently via https://github.com/spacetelescope/synphot_refactor/pull/387 . You should still be able to get synphot < 1.4
in Python 3.9 .
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).
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?
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!
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?
@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!
@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:)
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?
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.
Just for my own record, looks like you are only interested in this part of the blackbody curve.
Right, for this particular bandpass, but it should be able to handle arbitrary frequency/wavelength ranges and temperatures between ~50-500K
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
.
Also, in numpy 2.0, np.trapz
has been deprecated in favor of scipy.integrate.trapezoid
. FYI.
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!
Yes, the blackbody model might work too. Good luck!
from https://github.com/pyOpenSci/software-submission/issues/161:
Instead of manually performing bandpass integrations with the custom
Bandpass
type, we should use the Astropy affiliatedsynphot
package. This would also provide better integration with Astropy.