GriddedPSFModel evaluated PSFs with the same input flux have very different output fluxes

[york-stsci]

This has been tested with a GriddedPSFModel created out of a set of webbpsf Roman WFI apertures, using the sample code below. Note that instead of using webbpsf directly to create the PSF grid, instead the grid is created image-by-image and put into an NDData array after creation. This is because webbpsf currently does not preserve parity for even oversamples, and the simulator that I'm using these models for requires odd parity to build its EPSF models.

[york-stsci]

import numpy as np

from astropy.io import fits
from astropy.nddata import NDData
from photutils.psf import GriddedPSFModel
from webbpsf import wfirst as roman

def create_gridded_psf(fov_pixels, oversample, psf_grid):
    meta = {'oversampling': oversample, 'grid_xypos': []}
    pixel_scale = 0.11
    num_psf = psf_grid*psf_grid
    fov = fov_pixels * oversample
    if fov%2 == 0:
        fov += 1
    data = np.zeros((num_psf, fov, fov), dtype=np.float64)
    if num_psf == 1:
        psf_xgrid = [2048.]
        psf_ygrid = [2048.]
    else:
        pix_low, pix_high = 4., 4092.
        psf_xgrid = np.linspace(4., 4092., psf_grid)
        psf_ygrid = np.linspace(4., 4092., psf_grid)
    i = 0
    for yp in psf_ygrid:
        for xp in psf_xgrid:
            meta['grid_xypos'].append((xp, yp))
            ins = roman.WFI()
            ins.options['parity'] = 'odd'
            ins.filter = 'F087'
            ins.detector = 'SCA01'
            ins.detector_position = (xp, yp)
            ins.pixelscale = pixel_scale/oversample
            psf_fits = ins.calc_psf(fov_pixels=fov, oversample=1)
            psf = psf_fits["DET_SAMP"].data
            print("PSF at ({},{}) has total flux {}".format(xp, yp, np.sum(psf)))
            print("PSF Shape is {}".format(psf.shape))
            data[i][:] = psf * oversample**2
            i += 1
    return GriddedPSFModel(NDData(data, meta=meta))

grid = create_gridded_psf(65, 10, 5)
half_grid = 32
x_loc = [1234.56, 2198.75, 3114.03, 155.603, 1610.203, 2048.0, 2048.0]
y_loc = [1234.56, 729.14, 2541.89, 1338.08, 1338.04, 1338.04, 1338.08]
fluxes = [1000.0, 1000.0, 1000.0, 39.482, 39.482, 100.0, 100.0]

for (x,y,flux) in zip(x_loc, y_loc, fluxes):
    ix, iy = int(round(x)), int(round(y))
    lx = ix - half_grid
    hx = ix + half_grid + 1
    ly = iy - half_grid
    hy = iy + half_grid + 1
    xg, yg = np.mgrid[lx:hx, ly:hy]
    psf_img = grid.evaluate(xg, yg, flux=flux, x_0=x, y_0=y)
    print("Image at ({},{}) with flux {} has sum {}".format(x, y, flux, np.sum(psf_img)))

[larrybradley]

@york-stsci Many thanks for the code example. I'll try to take a look soon.