nucex

Separates a comet's nucleus from its coma in a space-telescope image. The nucleus is a point source, the coma a power law; both go through the telescope's PSF and are fitted together. nucex returns the nucleus's flux, its uncertainty, and a verdict in plain words on whether to believe it.

Free, and yours to use however it helps.

import nucex
from nucex import GaussianPSF, SyntheticComet

psf = GaussianPSF(3.0)                                    # FWHM 3 px; or ArrayPSF, EffectivePSF
comet = SyntheticComet(shape=(121, 121), psf=psf, gamma=1.0, kc=5.0, scheme="exact")
image, sigma, kn_true = comet.image(eta=0.30, seed=42)    # a known nucleus: 30 % of the light inside 15 px
res = nucex.extract(image, sigma, psf)
print(res.summary())

k_N = 199.77 ± 2.8 (stat 1.3, sys 2.5); eta15 = 29.83 %; gamma = 0.998 (inner 0.988); chi2/dof = 1.012
TRUSTED: nucleus share 29.8 %, k_N 72 sigma from zero, model fits
Not measurable: a coma core the PSF hides, holding 1 % of the coma's light inside 1 FWHM, would move k_N by 0.41 %.

k_N is the nucleus's total flux, in the image's units. On a FITS frame:

python -m nucex FRAME_flc.fits --x 2040 --y 1032 --stdpsf STDPSF_WFC3UV_F606W.fits --sky --json out.json

The model

I(x, y) = k_N P(x - x_c) + [ k_C (rho / 1 px)^-gamma (1 + azimuthal terms) ] * P  (+ sky)

with P the PSF and * convolution: the form of Hui & Li (2018)'s Eq. 2, with one slope for every azimuth (theirs may vary with azimuth).

• One fit, every pixel. The blurred coma and the nucleus are fitted together to every pixel inside r_fit, the core included. The classic technique fits the power law to an annulus of the blurred image, then blurs the fit; Hui & Li (2018, Sec. 3.4) trace its bias to that step.
• Exact pixels at the singularity. The mean of rho^-gamma over any rectangle comes from an exact formula, one smooth part of it integrated by Gauss-Legendre quadrature (nucex.coma.corner_integral). Every subpixel near the centre carries that mean, wherever the centre falls.
• Variable projection. The slope and the centre are searched; the fluxes (k_N, k_C, azimuthal amplitudes, sky) are solved exactly at each step. The errors come from the full Jacobian, every parameter free, scaled by the reduced chi^2 when it exceeds 1.
• The inward extrapolation, tested. Two more fits free the coma's slope inside 2 and 3 PSF FWHM. Each moves k_N and widens its error; the larger of the two, move and widening in quadrature, is kn_err_sys, and kn_err adds it to the statistical error in quadrature. inner_breaks=(1, 1.5, 2, 3) tests closer to the core, at a price: of 160 comae that are true power laws (Gaussian FWHM 1.5 and 3 px, γ 0.9 to 1.5, 20 noise draws each), with a nucleus at 30 % of the light, 46 lose TRUSTED (29 %), against 2 with the defaults.
• What no fit can measure, stated. core_ratio is the coma's light inside 1 FWHM over k_N: a coma core the PSF hides, holding 1 % of that light, would move k_N by core_ratio %.

PSFs: GaussianPSF(fwhm); ArrayPSF(data, oversample) for an optical PSF sampled finer than the detector (Tiny Tim's sub=n); EffectivePSF(data, oversample) or EffectivePSF.from_stdpsf(path, x, y) for an effective PSF, the pixel's response already in it, as in STScI's WFC3/UVIS library. An effective PSF is sampled at pixel centres, not integrated over them. The model's subsampling S (default 5) must equal an ArrayPSF's oversample or divide it by an odd factor: a 4× PSF needs S=4 (--S 4 on the command line).

The verdict

Checked in this order, each with the reason in the summary:

verdict	when
NOT TRUSTED	the nucleus holds less than eta_min (default 10 %) of the light inside 15 px, Hui & Li (2018)'s line; or that share is undefined
NOT DETECTED	k_N is under 3 sigma by the noise alone
CAUTION	k_N is under 3 sigma with the inner slope free; or kn_err_sys exceeds 10 % of k_N; or chi^2/dof > 1.5
TRUSTED	none of the above

Tested against

Every test here uses synthetic comets with known nuclei, and gives the fit their true noise maps.

Hui & Li (2018)'s bias map. nucex.hl18 follows the steps of their Sec. 2. Where their text leaves a choice open, ours is named in its docstring: the annulus fit weighted by 1/σ² from the image's noise map (their text gives that fit no weights; Hui, Jewitt, Mutchler, Agarwal & Kim (2026, ApJL 999, L37; arXiv:2601.21569) "assigned a weighting map based on Poisson statistics", following them), each sector's slope bounded to −1…3, 36 sectors smoothed by a median over 5 and joined by linear interpolation in azimuth, the origin at the peak pixel, the χ² region equal to the photometry aperture. nucex.synth builds synthetic comets on their recipe: their noise (Eq. 4), their μ at the singular centre (Eqs. 5-6), their nucleus share inside 15 px, and our choice of subsampling for the truth, 9. On Gaussian PSFs, their 13 curves (Figs. 2, 3 and 6, digitized from the vector paths of their PDF) come out in their order within each family (PSF width, coma slope, model subsampling) at a nucleus share of 1 %, over 8 noise fields, the same 8 for every curve and every share. Not curve for curve. One open choice was calibrated: the photometry aperture (2 FWHM), on one curve. The rest were set once and never tuned, and the other 12 curves stayed out of sample.

Comae that break the model. A compact core the PSF hides passes for nucleus. So does a coma that changes slope inside about 1.3 PSF FWHM: a steepening there makes the nucleus too bright, a flattening too faint. Measured with the defaults on Gaussian PSFs, FWHM 3 px unless noted, slope 1.0 outside the break, 20 noise draws a case:

the coma	nucleus share	k_N's mean error	stated error over the fitted k_N (median)	TRUSTED	of those, truth beyond 2 sigma
steepening, γ 1.0 → 1.25, inside 0.83 FWHM	10 %	+41.7 %	6.2 %	19 of 20	19
steepening, γ 1.0 → 1.25, inside 0.83 FWHM	30 %	+11.2 %	2.0 %	20 of 20	20
steepening, γ 1.0 → 1.15, inside 1.2 FWHM (FWHM 2 px)	30 %	+5.5 %	2.3 %	20 of 20	13
steepening, γ 1.0 → 1.3, inside 1.3 FWHM	30 %	+20.0 %	7.8 %	20 of 20	13
flattening, γ 1.0 → 0.7, inside 1.3 FWHM	30 %	−11.7 %	7.8 %	20 of 20	3
a compact core (FWHM 1.5 px) holding 5 % of the power law's light inside 3 px	10 %	+7.6 %	4.0 %	20 of 20	5

core_ratio sizes this risk: for one planted core (FWHM 1 px) holding 3.03 % of the coma's light inside 1 FWHM, the mean of core_ratio × 3.03 % came within 14 % of the mean error the core caused, over 20 images at each of shares 10 and 30 %. A 10 % share is a floor, not a guarantee.

A PSF slightly wrong in width. Given a Gaussian 2 or 5 % too narrow or too wide (the truth: FWHM 3 px, slope 1.0, 20 noise draws a case), with a nucleus at 30 % of the light, k_N moved by about as much as the width: −2.1, −4.9, +1.6 and +3.9 % for −2, −5, +2 and +5 %. The stated error, over the fitted k_N, grew to cover it (median 5.4 % at −2 %, 16.6 % at −5 %). Of the 240 images, at shares of 3, 10 and 30 %, 76 were TRUSTED, none with the truth beyond 2 sigma; at −5 % none was TRUSTED. Every fit PSF had the true shape.

scripts/validate.py makes the runs behind these numbers, one CSV row per image; scripts/hl18_map.py runs Hui & Li's technique over their grid. The suites on an HST effective PSF read the STDPSF file named by NUCEX_STDPSF.

On real frames

• Give extract the image, its 1-sigma map in the same units, and a PSF sampled finer than the detector, or an effective PSF.
• Clean cosmic rays first, or pass a mask (True = bad). Masked pixels leave the fit, and count as the fitted model's light in the nucleus share and core_ratio. The command line masks DQ bits 4, 16, 32, 256, 512, 2048, 4096 and 8192 (--dq-bad sets them).
• Set r_fit inside the radius where the coma is a power law and above the noise. Fit the sky (sky=True) unless it's already removed.
• azimuthal=1 or 2 adds azimuthal terms to the coma's brightness, one slope at every azimuth; free_coma_center=True lets the coma's centre sit off the nucleus. Check res.fit.residual.

Install and test

In this folder:

pip install '.[fits]'     # numpy and scipy; [fits] adds astropy for FITS files, the command line and STDPSF files
pip install '.[test]'     # adds pytest
pytest                    # seconds; the FITS test skips without astropy
pytest -m slow            # the noisy coverage test alone

Licence and citation

BSD 3-Clause. Cite Hui & Li (2018, PASP 130, 104501) for the technique's reliability study, and Lamy & Toth (1995, A&A 293) and Lamy, Toth & Weaver (1998, A&A 337, 945) for the technique.