Conversion between wavelength and frequency giving discrepant results

[tddesjardins]

Testing some synthetic photometry on Vega, I'm getting different results between pysynphot and synphot. I've verified that both pysynphot and synphot are using the same Vega spectrum and the bandpass ("acs,wfc1,f606w") is using the same files and plotting the filter throughput curves they look identical.

import pysynphot as psyn
import synphot as syn
import stsynphot as stsyn
import matplotlib.pyplot as plt
import numpy as np
import os
from astropy import units as u

# Set up bandpass
bandpass = 'acs,wfc1,f606w'
pband = psyn.ObsBandpass(bandpass)
sband = stsyn.band(bandpass)

# Set up vega spectra
pvega = psyn.FileSpectrum(os.path.join(os.environ['PYSYN_CDBS'], 'calspec', 'alpha_lyr_stis_008.fits'))
svega = stsyn.Vega

# Make observations
pvega_obs = psyn.Observation(pvega, pband)
svega_obs = syn.Observation(svega, sband, binset=sband.binset)

# Get fluxes
pflux = pvega_obs.effstim('flam')
sflux = svega_obs.effstim(syn.units.FLAM)

pflux_jy = pvega_obs.effstim('jy')
sflux_jy = svega_obs.effstim(u.Jy)

I've done this several different ways now, and find the following differences compared to the pysynphot result:

Type	PySynphot	Synphot	% Difference
Count rate	36395078438.946	36395792271.484	0.00196
FLAM	2.848e-9	2.848e-9	1.4531e-14
Jy	3331.797	3251.440	2.4118
FNU	3.332e-20	3.251e-20	2.4118
STmag	0.264	0.264	2.3166e-12
ABmag	0.093	0.120	28.4093

Changing the obsmode (trying ACS/WFC1 F435W, F606W, F814W, and WFC3/UVIS1 F218W, F814W) gives varying results but the frequency-based values are seemingly always off by ~1-2.5% (except ABmag, which has jumped between 2–28%).

Unless I'm doing something wrong with the my commands for calculating the flux in each unit (do I need to do an extra step when using frequency-based units?), then it seems like there's something wrong with the conversion between wavelength and frequency in one of the two packages, but it's not clear which one. @pllim suggested opening this issue on the synphot repo so that it can be investigated (and also suggested @bhilbert4 might be interested).

[bhilbert4]

Interesting. This could be affecting my NIRCam flux calibration values (although for my purposes 2.5% uncertainty is probably fine). I'll be following this to see what comes of it. And I'll try something similar for NIRCam tomorrow to see if I find any differences.

[pllim]

In IRAF SYNPHOT, it only allowed calculation in PHOTLAM or FLAM. I vaguely remember encountering this problem of applying effstim to arbitrary units before. I need to refresh my memory of how I handled it in the refactoring... Thanks for reporting this!

[pllim]

Like Tyler and I discussed offline, the formula for effstim (see https://synphot.readthedocs.io/en/latest/synphot/formulae.html#effective-stimulus) only takes account of F_lambda integration across lambda. If you change F_lambda to F_nu in that formula, it would be wrong unless you also include the lambda to nu conversion factor. Again, need some thinking! 🤔

[pllim]

To add to the excitement, if you convert effstim calculated in FLAM directly to Jy, it doesn't match pflux_jy nor sflux_jy.

>>> sflux.to(u.Jy, u.spectral_density(svega_obs.effective_wavelength()))
<Quantity 3209.37237987 Jy>

[tddesjardins]

That I'm not too surprised about because it's probably having to assume something about the spectral shape when it does that, rather than directly measuring it from the spectrum. Detail about the input spectrum has been lost.

[pllim]

So, after some attempts at calculus, I arrived at this formula for calculating effstim in frequency space (feel free to double check!):

(INT F_nu * P_nu * -(c^2 / nu^3) d_nu ) / (INT P_nu * -(c^2 / nu^3) d_nu)

Hence, the following attempt at calculating effstim in frequency space using what @tddesjardins had above:

>>> from synphot import units
>>> nu = svega_obs.binset.to(u.Hz, u.spectral())
>>> fac = -1 * units.C * units.C / (nu * nu * nu)
>>> fac
<Quantity [-3.33564095e-10, -3.34565789e-10, -3.35569485e-10, ...,
           -4.43731687e-07, -4.43852738e-07, -4.43973811e-07] Angstrom2 / (Hz3 s2)>
>>> influx = svega_obs(nu, flux_unit=u.Jy)  # F_nu * P_nu
>>> num = np.trapz(influx * fac, x=nu)
>>> den = np.trapz(fac * sband(nu), x=nu)
>>> effstim_jy = num / den  # Theoretical effstim calculated in freq space
>>> effstim_jy
 <Quantity 3251.4387289 Jy>
>>> effstim_jy - sflux_jy  # yay?
<Quantity -0.00152408 Jy>

Based on this, I would like to say that synphot gives correct result but I would be happy if you can confirm or deny this claim!

[pllim]

More thoughts here: https://gist.github.com/pllim/ddbe9f0c4bdeb920b0501215fb583c72

Would be interested to hear what you think of this, @tddesjardins and @bhilbert4 .

[bhilbert4]

What you say in the gist makes sense to me. Somewhere in effstim when working with PHOTLAM values, something is wrong. The units of the observation being put into effstim don't seem to matter, but and only when you request the answer in PHOTLAM does the error show up.

[tddesjardins]

I think this was a nice check and indicates that PySynphot is the thing that's incorrect. This is both reassuring AND troubling at the same time!

[bhilbert4]

Maybe it's time to give pysynphot a gold watch and congratulate it on its retirement? :)

[pllim]

And just for completeness, effstim implementation between pysynphot and synphot is completely different. I cannot say I understand what pysynphot attempted to do in all details.

[pllim]

So... Should we close this issue?

[tddesjardins]

If we're confident that synphot is correct, then I say we close it, but maybe an issue should be opened on pysynphot that references this so we don't forget about it.

[pllim]

Done at spacetelescope/pysynphot#96

[pllim]

After further investigation by @tddesjardins and Susana Deustua, it turns out the bug is in synphot after all. More below:

>>> import pysynphot as S
>>> sp_pysyn = S.Vega
>>> bp_pysyn = S.ObsBandpass('acs,wfc1,f555w')
>>> obs_pysyn = S.Observation(sp_pysyn, bp_pysyn)
>>> obs_pysyn.effstim('flam')
3.786643843663875e-09
>>> obs_pysyn.effstim('fnu')
3.630073134871726e-20

>>> from astropy import constants as const
>>> import synphot
>>> import stsynphot
>>> sp_stsyn = stsynphot.Vega
>>> bp_stsyn = stsynphot.band('acs,wfc1,f555w')
>>> obs_stsyn = synphot.Observation(sp_stsyn, bp_stsyn, binset=bp_stsyn.binset)
>>> obs_stsyn.effstim('flam')
<Quantity 3.78664384e-09 FLAM>
>>> obs_stsyn.effstim('fnu')  # Wrong
<Quantity 3.60801141e-20 FNU>
>>> w = bp_stsyn.pivot()
>>> # According to Susana, this one is the correct way and agree with Pysynphot
>>> (obs_stsyn.effstim('flam') * w * w / const.c.cgs).to(synphot.units.FNU)
<Quantity 3.63007313e-20 FNU>

[pllim]

Might also have discrepancies with effstim('vegamag') although effstim('obmag') seems fine. Using Tyler's filters and the two stars, max discrepancy for both ABMAG and VEGAMAG is around 0.05 mag.

p.s. I'll post a Gist of my findings once I can understand more.

[pllim]

See #166 (the promised Gist is also linked over there).

[pllim]

I now think there are bugs on both sides... 🤣 😭 spacetelescope/pysynphot#103

[pllim]

Conclusion: Susana derived effstim and a few other properties from the first principles. From discussions with @deustua and @tddesjardins today, the following is the correct way to fix this:

1. Wavelength in the effstim formula is the bandpass pivot wavelength, not the observation pivot wavelength. In fact, there is no such thing as the latter.
2. effstim should be calculated using the formula in the unit of FLAM. Then it is converted to other requested unit using Astropy's unit conversion and bandpass pivot wavelength.

When Susana publishes her work, it will be mentioned in the doc of this package as a reference.