Optics: scale PSF cache radius by parameters

psf_lookup_optimization_plan.md

@@ -66,10 +66,12 @@ flux so that both of these hold:
 - the PSF value at the support boundary is below a configured per-pixel HDR
   error budget.
 
-The cache contains weights through the largest radius required by the selected
-configuration.  A dim image uses a cropped subset; an unusually bright image
-uses a larger subset.  If an event requires a radius beyond the cache maximum,
-use the direct Moffat evaluator for that event and record it in diagnostics.
+The cache contains weights through the radius implied by the selected PSF
+parameters at the reference scaled flux of one.  A dim image uses a cropped
+subset; an unusually bright image uses a larger subset.  If an event requires
+a radius beyond that cache maximum, use the direct Moffat evaluator for that
+event and record it in diagnostics.  The radius is not a fixed pixel cap:
+low-beta Moffat profiles need broad wings even when their FWHM is small.
 This prevents a visible hard cutoff in bright-star wings without silently
 discarding their flux.
 

src/optics.c

@@ -1,5 +1,6 @@
 #include "optics.h"
 
+#include <limits.h>
 #include <math.h>
 #include <stdint.h>
 #include <stdio.h>
@@ -80,7 +81,6 @@ LinearRgb blackbody_to_linear_rgb(double temperature_K)
 
 enum {
     PSF_PHASE_RESOLUTION = 64,
-    PSF_MAX_CACHE_RADIUS_PIXELS = 48,
     PSF_QUADRATURE_ORDER = 4,
 };
 
@@ -195,8 +195,14 @@ int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf)
     const double requested_radius = moffat_support_radius(alpha, psf->moffat_beta, 1.0);
     if (!isfinite(requested_radius) || requested_radius <= 0.0)
         return -1;
-    const int radius = (int)fmin((double)PSF_MAX_CACHE_RADIUS_PIXELS,
-                                 ceil(requested_radius) + 1.0);
+    /* The ordinary cache must cover the support implied by this invocation's
+     * PSF parameters.  In particular, a low beta has broad Moffat wings even
+     * when its FWHM is small; a fixed radius would make every ordinary event
+     * fall back to the much slower direct evaluator.  Brighter events can
+     * still exceed this flux=1 reference support and use that fallback. */
+    if (requested_radius > (double)INT_MAX - 1.0)
+        return -1;
+    const int radius = (int)ceil(requested_radius) + 1;
     const size_t nodes = PSF_PHASE_RESOLUTION + 1u;
     const size_t side = (size_t)radius * 2 + 1u;
     if (nodes > SIZE_MAX / nodes || nodes * nodes > SIZE_MAX / side ||