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)
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.
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:This outputs:
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
ChromaticTransformationclass, specifically in the sed function:Having the sed multiply the absolute value of the determinant rather than simply the determinant seems to fix the issue.