Is effstim wrong when calculated in PHOTNU and PHOTLAM?

[pllim]

This does not make sense to me:

>>> import pysynphot as S
>>> S.__version__
'0.9.12'
>>> bp = S.ObsBandpass('acs,wfc1,f555w')
>>> sp = S.Vega
>>> obs = S.Observation(sp, bp)
>>> bp.pivot()
5360.938305731454
>>> obs.effstim()
1024.232011041325
>>> obs.effstim('photnu')
8.744735882086226e+28
>>> S.units.Photlam().ToPhotnu(bp.pivot(), obs.effstim())
9.818819146083892e-09
>>> obs.effstim('jy')
3630.073134871727
>>> S.units.Photlam().ToJy(bp.pivot(), obs.effstim())
3638.3434633950515

As shown above, effstim() calculated in Jy is not that far off compared to effstim() calculated in PHOTLAM and then converted to Jy around bandpass' pivot wavelength. However, that is not the case for PHOTNU. Why?

@deustua and @tddesjardins , since you have been looking into this from the INS side, what do you think? Is this a bug or am I not understanding something about effective stimulus in PHOTNU?

git blame shows that the Photnu.unitResponse() method used in the calculation in question was implemented by @vglaidler 11 years ago and has no changed since.

[pllim]

Looking at this more and referring to @deustua 's "first principle" formula that states effstim in PHOTLAM should be the same as effstim in FLAM converted to PHOTLAM along bandpass' pivot wavelength, I think PHOTLAM values might be problematic as well but at a smaller scale.

>>> from astropy import units as u
>>> import pysynphot as S
>>> S.__version__
'0.9.12'
>>> bp = S.ObsBandpass('acs,wfc1,f555w')
>>> sp = S.Vega
>>> obs = S.Observation(sp, bp)
>>> w = bp.pivot()
>>> effstim_flam = obs.effstim('flam')
>>> effstim_photlam = obs.effstim()
>>> x = (effstim_flam * (u.erg / u.s / u.cm ** 2 / u.AA)).to(u.ph / u.s / u.cm ** 2 / u.AA, u.spectral_density(w * u.AA))
>>> x
<Quantity 1021.92386949 ph / (Angstrom cm2 s)>
>>> print('{:.2%}'.format(abs(x.value - effstim_photlam) / effstim_photlam))
0.23%