Negative flux of Roman PSF when using calculateFlux() #1346
Description
Encounter negative flux when using the calculateFlux() function to get the flux of an object containing the Roman PSF (with a specified wcs). This problem has been discussed with @arunkannawadi. Here is the code to reproduce the problem:
import galsim
from galsim import roman
import numpy as np
wcs_dict = galsim.roman.getWCS(world_pos=galsim.CelestialCoord(ra=-70*galsim.degrees, dec=-50*galsim.degrees), PA=1.*galsim.degrees, SCAs=range(1, 18), PA_is_FPA=True)
sed = galsim.SED( galsim.LookupTable( [1000, 26000], [1, 1], interpolant='linear' ),
wave_type='Angstrom', flux_type='fphotons' )
flux = 1
SCA = 1
filter_name = 'H158'
bp= roman.getBandpasses(AB_zeropoint=True)[filter_name]
wcs = wcs_dict[SCA]
SCA_pos_center = galsim.PositionD (2044, 2044)
roman_psf = roman.getPSF(SCA, filter_name, pupil_bin=8,wcs = wcs, SCA_pos = SCA_pos_center)
roman_psf_nowcs = roman.getPSF(17, filter_name, pupil_bin=8, SCA_pos = SCA_pos_center)
point = ( galsim.DeltaFunction() * sed ).withFlux( flux, bp )
print("point flux:", point.calculateFlux(bp))
psf = galsim.Convolve(point, roman_psf)
psf_nowcs = galsim.Convolve(point, roman_psf_nowcs)
print("Flux (with wcs):", psf.calculateFlux(bandpass=bp))
print("Flux (no wcs):", psf_nowcs.calculateFlux(bandpass=bp))
stamp = galsim.Image(1000,1000)
psf.drawImage(bp, method="no_pixel", wcs=wcs,
use_true_center=True, image=stamp)
print("Flux image (WCS specified):", np.sum(stamp.array))
stamp = galsim.Image(1000,1000)
psf.drawImage(bp, method="no_pixel",
use_true_center=True, image=stamp)
print("Flux image (no WCS specified):", np.sum(stamp.array))
local_wcs = wcs.local(SCA_pos_center)
print("WCS jacobian determinant:", local_wcs._det)
This outputs:
point flux: 1.0
Flux (with wcs): -1.0
Flux (no wcs): 1.0
Flux image (WCS specified): 0.9079357
Flux image (no WCS specified): 0.9881175
WCS jacobian determinant: -0.01185563608937022
The problem is coming from the fact that the determinant of the WCS jacobian is negative. Internally within galsim the jacobian determinant gets multiplied by SED in the ChromaticTransformation class, specifically in the sed function:
def sed(self):
sed = self._original.sed * self._flux_ratio
if self._redshift is not None:
sed = sed.atRedshift(self._redshift)
# Need to account for non-unit determinant jacobian in normalization.
if hasattr(self._jac, '__call__'):
@np.vectorize
def detjac(w):
return (np.linalg.det(np.asarray(self._jac(w)).reshape(2,2)))
sed *= detjac
elif self._jac is not None:
sed *= (np.linalg.det(np.asarray(self._jac).reshape(2,2)))
return sed
Having the sed multiply the absolute value of the determinant rather than simply the determinant seems to fix the issue.