Spectral Flux Color Calibration (SFCC)

Physics-grounded colour calibration. SASpro detects stars in the plate-solved frame, matches each one to a Gaia DR3 XP calibrated spectrum (or a Pickles library template as fallback), integrates its spectrum through the telescope's filter/sensor/light-pollution curves to predict the expected R/G, B/G ratios, measures them in the image, fits a correction model, and applies it.

1. Using the tool

The dialog is a two-step workflow with a shared panel of filter/sensor selections between them. Both steps run against the active image window; switching the active window between runs is fine.

Step 1 — Fetch stars & spectra

The White Reference dropdown picks the SED that defines “neutral grey” (default G2V; Vega / A0V is the usual alternative). It only affects the reference integrals S_ref_c, nothing else at this stage.

Click Step 1: Fetch Stars & Spectra from Current View and the pipeline runs:

SEP detection on the frame's luminance. The Auto σ checkbox (default on) runs the waterfall described in §2.1; uncheck and dial the Star detect σ spinbox to force a value. Detections are kept only when their flux-radius is in (0.2, 10] px.
WCS projection of every surviving source to RA/Dec, using the header's own WCS (deprecated RADECSYS and EPOCH keywords are normalized first).
Bulk Gaia match against SASpro's local Gaia library at a 3 arcsec tolerance; each hit carries a Gaia source id for XP lookup later.
SIMBAD region query over a circle covering the field (with retries), to tag matched stars with spectral type, main id, and Johnson B/V/R magnitudes.
Pickles template match for every SIMBAD spectral type (§2.3). Stars without a spectral type but with B and V get their letter class inferred from B−V (§2.3).
For Gaia-matched stars that still have no Pickles template, SASpro computes synthetic Johnson B/V/R from the Gaia XP spectrum via gaiaxpy.generate(system=JKC), caches it, and infers the letter class from the synthetic B−V.
Stars with no main id but a Gaia id get a SIMBAD ids lookup in batches of 200 to resolve a Gaia DR3 source id.

The result is a summary plot in two panels. The status line below shows the raw counts; the plot shows the quality of what you got.

Left: spectral-type histogram. Blue bars are stars with a Pickles template matched from a SIMBAD type; orange bars overlay stars whose letter class came from B−V inference (SIMBAD or Gaia synthetic) because no spectral type was published. A healthy astrophotography field is bell-shaped around F / G / K with tails to A and M; a histogram that's all one letter usually means the SIMBAD query failed or the field is pathologically narrow (star cluster, galactic plane crowding). Blue dominating is good; a lot of orange means SASpro had to fall back to colour-only typing for most of the field.
Right: calibration-source pie. Four mutually-exclusive buckets: Gaia XP (green — the star has its own calibrated spectrum, the gold standard), Pickles via SIMBAD (blue — typed by SIMBAD, approximated with a library SED), Pickles via B−V (orange — typed from colour, same approximation) and no spectrum (grey — dropped at Step 2). More green is strictly better: Gaia XP uses the star's actual spectrum. If grey dominates, your field is outside Gaia DR3's XP coverage or your local Gaia library doesn't have the right spatial group installed — Step 2 will still run, but on fewer stars.

The shared panel — filter/sensor selections

Pick the actual optical chain: R Filter, G Filter, B Filter, Sensor (QE) and up to two LP/Cut Filter curves (e.g. an IR-cut and a light-pollution filter). (None) is a flat-unity curve. These are read from the shipped sasp_data.fits plus a per-user usercustomcurves.fits. Choices persist via QSettings.

Add Custom Filter/Sensor Curve… — import a 2-column CSV (λ_nm, response) or digitize a scanned curve from an image (click bottom-left, bottom-right, top-left on the plot axes, enter λ_min and λ_max, SASpro extracts the darkest-pixel-per-column trace, interpolates anomalies, median-smooths, and writes it to the user FITS as a FILTER or SENSOR extension).
Remove Filter/Sensor Curve… — rewrites usercustomcurves.fits without the chosen extension.

The SASP Viewer — see your optical chain against any star

Click Open SASP Viewer (next to the Step 1 button) to launch a standalone window that plots any stellar SED through any R/G/B filter set, any sensor QE, and up to two LP/cut filters. It doesn't touch the image — it's a free-form exploration of what your optical chain actually sees, and it's where you get a feel for why one filter choice gives you deeper red stars than another, or why that bargain CMOS sensor can't tell a late K from an M.

The dropdowns across the top are the same vocabulary as the main dialog: Star Template (any Pickles SED from O supergiants through late M dwarfs), R / G / B Filter, two optional LP/Cut Filter slots, and Sensor (QE). Change any dropdown and press Plot to redraw.

What the plot shows

One figure, two y-axes:

Left y-axis — stellar flux (erg s⁻¹ cm⁻² Å⁻¹). The black line is the raw SED of the chosen star — a hot O star peaks in the UV and plummets by Hα; a cool M dwarf is nearly dark at B wavelengths and climbs through red and near-IR. The gold lines are the per-channel responses — the SED multiplied by that channel's full system throughput. Where the gold peaks is where that channel's signal actually comes from for this star — often shifted away from the filter's nominal centre by the sensor's QE roll-off.
Right y-axis — relative throughput (0–1). The dashed red / green / blue curves are the system throughputs T_filter × QE × LP₁ × LP₂ alone, independent of the star. These are the instrument; the gold curves are what the instrument does to the star. If the dashed curves look clean and symmetric but the gold curves are chopped off at one side, the sensor's QE is choking that end of the band.
Title — Vega-relative magnitudes. SASpro computes each channel's magnitude relative to Vega (an A0V) through the same chain and prints them above the plot (e.g. R→SII: 3.42 | G→OIII: 4.11). These are system magnitudes you can sanity-check against a known calibrated star; a large asymmetry (one channel much fainter than the others for a G2V) is a clue that your filter balance is off.

Things worth doing with it

See why cool stars are blue-poor. Pick M5V, leave an ordinary R+G+B+IR-cut chain selected: the gold B response is a flat line against the floor while R has a tall peak. That's why red giants look blue-starved before any colour correction.
Compare two filter brands on the same sensor. Pick a stable SED (G2V), switch between an Antlia and a Chroma R filter, and watch how the gold curve's shape and position move. The area under the gold curve is the signal, and the magnitude in the title turns that into a number you can compare side-by-side.
Preview an LP filter before you buy. Add a dual-narrowband or broad LP filter in one of the LP slots; the dashed system curves get carved out at the stopbands, and the gold curves show which stellar flux survives. For broadband work this immediately tells you how much star signal you're giving up for the sky-glow suppression.
Sanity-check a custom curve you just imported. After Add Custom Filter/Sensor Curve… runs, open the viewer and plot it against an A0V. If the dashed curve has spikes, dropouts, or a wrong-side slope, the digitize/CSV import captured something wrong and you can re-run before trusting it in Step 2.
Explore where your green channel actually sits. Change the star from early B to late M in sequence; the gold G response shifts right with cooler stars because the SED is tilting toward red through the same green filter. This is a direct intuition for why a per-channel colour correction needs to depend on the star's own colour — which is what Step 2 does.

Step 2 — Run color calibration

With Step 1 done and the chain picked, click Step 2: Run Color Calibration. SASpro re-detects stars on the live frame, matches them within 3 px to the Step-1 star list, caps the sample to the brightest 500 by SEP flux, performs annulus-background-subtracted aperture photometry per channel, computes the model (§2.4) and applies it (§2.5).

Background Neutralization — runs the standalone BN tool first (so the colour fit isn't fighting a channel pedestal) and again on the final image.
3×3 Color Matrix Correction — hidden by default; when enabled, fits an extra cross-channel matrix after the per-channel scaling (§2.6).
A dialog reports the model chosen (slope-only / affine / quadratic), the star count split by source, and the net R and B scalings at x = 1.

Reading the before/after residuals plot

Two side-by-side panels, both plotting the fractional residual (measured_ratio / expected_ratio) − 1 on the y-axis against the expected ratio on the x-axis. Red circles are the R/G channel, blue squares are B/G. The shaded bands behind each colour are the Q1–Q3 interquartile range — the central half of the scatter.

BEFORE (left panel) — the raw situation. A horizontal cloud centred on zero means the camera is already close to right. A sloped trend (points high on the left, low on the right or vice versa) says the camera's response doesn't scale linearly with the true colour — precisely what SFCC is fixing. A big vertical offset (whole cloud shifted away from zero) is a flat colour cast, trivially corrected by any of the three models.
AFTER (right panel) — the same stars after the correction. Both colours should collapse to a tight horizontal band around zero. If they do, calibration worked. If a residual tilt or curvature remains, the chosen model couldn't fit the shape — usually because the camera's response has a sharp spectral feature the quadratic can't track (narrowband interference, uncoated filter wobble). A big vertical spread with no trend means noisy photometry (short exposure, undersampled stars) rather than a bad model.
Shaded Q1–Q3 band shrinking vs the before panel is the clearest sign the fit removed systematic colour error. If it doesn't shrink, the fit found nothing it could straighten — your frame was already colour-flat, or the model is wrong for your data.

Gradient extraction (Beta, hidden)

A separate Run Gradient Extraction button and method combo (poly2, poly3, rbf) are hidden by default but available on code paths that enable them. It uses the same spectra and photometry to measure a per-channel spectrophotometric gradient and remove it — see §2.7.

The result pane shows three side-by-side panels, one per channel, each a divergent red/blue heatmap of the fitted gradient surface in magnitudes. Panel titles report the peak magnitude for that channel so you can tell at a glance how aggressive the correction was — a 0.03 mag peak is a cosmetic tweak, a 0.3 mag peak was clamped (§2.7) because the fit tried to go further. Yellow dots overlay the star positions used for that channel's fit after sigma-clipping; sparse corners mean the surface is extrapolating there and should be treated with caution. The three channels usually show similar shape (same corner darker, same edge brighter) with different amplitudes — that's a real optical gradient. If the three channels show unrelated patterns, the fit is chasing noise and poly2 is probably the right method to retry with.

2. Technical details

RGB only; no narrowband. The core idea: use each star's real spectrum as a photometric standard. SIMBAD + Pickles can approximate the spectrum from a type label; Gaia XP gives you the actual observed spectrum of that star. SFCC uses XP where it has it and falls back to a hierarchical Pickles match where it doesn't.

Third-party data used. SFCC relies on Gaia DR3 XP calibrated spectra (ESA/DPAC) via gaiaxpy, the Pickles 1998 stellar flux library, SIMBAD (CDS Strasbourg), and SEP (Source Extractor as a Python library). The calibration pipeline, the three-model fit, the pivot-at-median application, and all SASpro-specific heuristics below are implemented in SASpro.

2.1 Star detection — Auto-σ waterfall

SEP runs on the mean of the three channels. Rather than fixing a detection threshold:

\sigma\in\{100,50,25,12,8,5,3\}\,,\;\text{stop at first}\;N\geq 1200

If no rung reaches 1200 the lowest σ is kept so there's still something to work with. Detections are filtered by flux-radius to (0.2, 10] px to drop galaxies/plate flaws. Only sep.extract is retried per rung — the one-off background estimation runs once.

2.2 System throughput & expected integrals

The optical chain combines multiplicatively on a shared 3000–11000 Å 1-Å grid:

T^{c}_{\mathrm{sys}}(\lambda)=T^{c}_{\mathrm{filter}}(\lambda)\cdot QE(\lambda)\cdot T_{\mathrm{LP}_1}(\lambda)\cdot T_{\mathrm{LP}_2}(\lambda)

For each star, the expected integrated flux through channel c is

S^{c}_{\mathrm{star}}=\int f_{\mathrm{star}}(\lambda)\,T^{c}_{\mathrm{sys}}(\lambda)\,d\lambda

where f_star(λ) is either the Gaia XP calibrated spectrum (interpolated onto the grid) or the matched Pickles SED. Integrals are cached per (source_id, system_key) in an SQLite store so repeat runs are instant — system_key is an MD5 over the wavelength grid and the three throughputs so any filter change triggers a recompute. Reference integrals S_ref_c are the same formula applied to the chosen white reference (default G2V).

2.3 Spectral-type → Pickles template

SIMBAD spectral type (e.g. K2III) is matched against the available Pickles extensions by letter class, subclass, and luminosity — exact match first, then same letter + same luminosity + nearest subclass, then same letter + any luminosity + nearest subclass, finally all templates in that letter. When SIMBAD gives no type but SIMBAD or Gaia-synthetic Johnson magnitudes are available, the letter is inferred from B−V:

letter = B if B-V < 0.00 A if 0.00 ≤ B-V < 0.30 F if 0.30 ≤ B-V < 0.58 G if 0.58 ≤ B-V < 0.81 K if 0.81 ≤ B-V < 1.40 M if B-V ≥ 1.40

and the match runs again on the letter alone (so a B−V-only star picks, say, the main-sequence K5V rather than K5III).

2.4 Photometry & the three-model fit

Step-1 catalogue stars are re-matched to Step-2 SEP detections within 3 px and capped to the brightest 500. For each kept star the aperture radii scale with the SEP semi-major axis a:

r=\mathrm{clip}(2.5\,a,\,2,\,12)\,,\;r_{\mathrm{in}}=\mathrm{clip}(3r,6,40)\,,\;r_{\mathrm{out}}=\mathrm{clip}(5r,r_{\mathrm{in}}{+}2,60)

Each channel's background-subtracted flux comes from an annulus-estimated pedestal:

F_c=\sum_{\mathrm{ap}}I_c-\mu^{c}_{\mathrm{bg}}\cdot A_{\mathrm{ap}}\,,\qquad \mu^{c}_{\mathrm{bg}}=\frac{\sum_{\mathrm{ann}}I_c}{A_{\mathrm{ann}}}

Photometry runs in parallel across up to 8 threads; stars with any non-positive channel are dropped. The result is the two measured ratios per star compared against the predicted ratios from §2.2:

\mathrm{meas}_{R/G}=\frac{F_R}{F_G}\,,\quad \mathrm{exp}_{R/G}=\frac{S^{R}_{\mathrm{star}}}{S^{G}_{\mathrm{star}}}

SASpro then fits three models independently for the R and B channels and keeps the best per channel. The slope-only solution is closed-form:

m=\frac{\sum_i\mathrm{meas}_i\,\mathrm{exp}_i}{\sum_i\mathrm{meas}_i^{2}}

and affine / quadratic are standard least-squares and polynomial fits:

\mathrm{exp}=m\cdot\mathrm{meas}+b\quad(\text{least-squares})

\mathrm{exp}=a\cdot\mathrm{meas}^{2}+b\cdot\mathrm{meas}+c

Each is scored by the fractional-residual RMS:

\mathrm{RMS}_{\mathrm{frac}}=\sqrt{\frac{1}{N}\sum_i\!\left(\frac{\mathrm{pred}_i}{\mathrm{exp}_i}-1\right)^{\!2}}

and the winner minimises the sum over both channels:

\min_{M\in\{\mathrm{slope},\mathrm{affine},\mathrm{quad}\}}\;\left(\mathrm{RMS}^{R}_{\mathrm{frac}}+\mathrm{RMS}^{B}_{\mathrm{frac}}\right)

A fractional loss is used, not absolute, because stars span orders of magnitude in flux and residuals should be measured relative to the expected value. Quadratic is only entered when there are ≥ 6 valid stars.

2.5 Applying the correction — pivot at median

With the winning per-channel polynomials, the per-pixel correction factors are the polynomials evaluated at the measured pixel ratio (x_R = I_R/I_G, x_B = I_B/I_G) and clipped to a safe dynamic range:

m_c(x)=\mathrm{clip}\!\left(a_c x^{2}+b_c x+c_c,\;0.25,\;4.0\right)

And applied with a pivot at the channel median rather than at zero, so the overall channel brightness is preserved and only the colour balance moves:

\mathrm{out}_c=p_c+(I_c-p_c)\cdot m_c(x_c)\,,\quad p_c=\mathrm{median}(I_c)

Green is the reference channel and is left unchanged. The clip at [0.25, 4] prevents extreme ratios on bright or dark pixels from driving runaway gains (e.g. a saturated red core where G→0).

2.6 Optional 3×3 color matrix

Off by default and hidden in the current UI. Rarely needed once per-channel scaling is applied; mostly useful for sensors with heavy cross-channel bleed. Not recommended for narrowband.

In ratio space, with G normalised to 1, each star gives a measurement vector and a target vector:

M · [F_R/F_G, 1, F_B/F_G]^T ≈ [S_R/S_G, 1, S_B/S_G]^T

SASpro solves this by least-squares row-by-row (per output channel), then forces the G row to [0, 1, 0] so green maps to itself exactly and the matrix stays near-identity. 3σ-MAD clipping on residuals removes outliers, and the solution is rejected if any of:

condition number κ(M) > 50,
|M_ii| > 3 for any diagonal entry, or
|M_ij| > 1.5 for any off-diagonal entry.

Applied per pixel in ratio space with G held fixed, then a per-pixel negative fix (shift up by min(channels)) and a luminance-pivoted ceiling compression keep the output in gamut without hard-clipping any one channel.

2.7 Gradient extraction (Beta)

Removes a spatially-varying colour gradient per channel using the same star catalogue. For each star and channel the differential magnitude is computed (no annulus — raw aperture only, since the gradient is the pedestal we're after), then field-median-subtracted so only the spatial variation remains:

\Delta m_c(x,y)=-2.5\log_{10}\!\left(\frac{F_c}{S^{c}_{\mathrm{star}}}\right)\;-\;\mathrm{median}(\Delta m_c)

A 2.5σ-MAD clip removes outliers, then a 2D surface (poly2, poly3, or RBF with a thin-plate-spline kernel) is fit per channel, re-centred, and clamped to ±0.3 mag so a bad fit can't blow up the image. The correction is multiplicative in flux:

\mathrm{scale}_c(x,y)=10^{-0.4\,\Delta m_c(x,y)}\,,\quad \mathrm{out}_c=I_c/\mathrm{scale}_c

At least 10 stars are required post-clipping, and at least 5 per channel to proceed with the fit.

2.8 DR4 roadmap (not shipped)

The current pipeline trusts the manufacturer's published QE curve. That curve is a room-temperature bare-die measurement; a real system's response is the product of sensor QE, cover glass, AR coatings, telescope throughput, atmosphere and airmass. With Gaia DR4's ~10× denser XP coverage, the plan is to drop the QE combobox entirely and solve the effective system response R(λ) directly from the stellar population — the same approach Gaia uses internally for its own photometric calibration (Carrasco et al. 2021, Riello et al. 2021). Filters stay as inputs; QE becomes a solved output, with the matrix of star / channel observations providing the overconstrained system.