Add Rcpp speedup for disp_stars convolution in SFHfunc() and SFHburst()

[Copilot]

The disp_stars=TRUE path in SFHfunc() and SFHburst() was bottlenecked by a pure-R loop over Gaussian grid points with repeated approx() calls per wavelength. This replaces that loop with a compiled C++ implementation.

C++ implementation (src/disp_stars.cpp)

New function disp_stars_cpp() replicates the R algorithm exactly:

• For each grid point: shift log-wavelength grid by log10(1 + grid[i]*z_disp), interpolate back onto the original grid via a monotone two-pointer sweep (equivalent to approx(..., rule=2, yleft=y[1], yright=y[n])), accumulate weights[i] * 10^y_interp
• Multiply by res at the end
• Pre-allocated buffers reused across grid iterations; guards against 1 + grid[i]*z_disp <= 0

Wiring (R/SFH.R)

Both disp_stars blocks reduced from ~10 lines to one call:

lum = .disp_stars_cpp(wave_lum_log, lum_log, z_disp, grid, weights, res)

All surrounding setup (grid/weights, LSF handling, z_disp calculation) is unchanged.

Package plumbing

• src/RcppExports.cpp — registered _ProSpect_disp_stars_cpp
• R/RcppExports.R — added .disp_stars_cpp() wrapper
• DESCRIPTION and NAMESPACE already carried Rcpp; no changes needed

Regression test (tests/test_disp_stars_cpp.R)

Compares C++ output against the original R loop for four (veldisp, range, res) combinations with fixed-seed synthetic spectra, asserting max relative error < 1e-10.

Original prompt

Create a pull request implementing an Rcpp speed-up for the disp_stars convolution in SFHfunc() and SFHburst() in ProSpect.

Context:

• Repository: asgr/ProSpect
• Target file: R/SFH.R (see https://github.com/asgr/ProSpect/blob/master/R/SFH.R)
• disp_stars=TRUE is common and is currently slow due to an R loop over grid and repeated approx() calls.

Requested change:

1. Add an Rcpp implementation of the dispersion step (log-lambda convolution with wavelength-dependent z_disp) that matches the existing algorithm:

   • Inputs: wave_log (log10 wavelengths), lum_log (log10 luminosity), z_disp vector, grid vector, weights vector, res scalar.
   • For each grid point, compute x_src = wave_log + log10(1 + grid[i]*z_disp) and y_src = lum_log - log10(1 + grid[i]*z_disp).
   • Linearly interpolate y_src(x_src) onto wave_log using a monotonic two-pointer sweep (equivalent to approx(..., rule=2, yleft=y_src[1], yright=y_src[n])).
   • Accumulate out += weights[i] * 10^y_interp.
   • Multiply final result by res.
   • Avoid per-iteration allocations by reusing buffers.
   • Add basic input validation (e.g., ensure 1 + grid[i]*z_disp[j] > 0 to avoid NaNs; behavior should be consistent with current R implementation—either error with clear message or allow NaNs).

2. Wire this into SFHfunc() and SFHburst():

   • Replace the inner for(i in seq_along(grid)) loops in the disp_stars blocks with a call to the new C++ function.
   • Keep the surrounding computation (grid/weights, LSF handling, z_disp calculation) the same.

3. Package plumbing:

   • Ensure ProSpect builds with the new compiled code.
   • Add/update DESCRIPTION (Imports/LinkingTo) and NAMESPACE as needed for Rcpp.
   • Add src/ code, and any required Rcpp boilerplate (RcppExports.cpp, RcppExports.R) using standard Rcpp attributes.
   • Add minimal documentation/comments explaining the new function and equivalence.

4. Add a small unit/regression test (if the repo uses testthat) or at least an R script under tests/ verifying that for a fixed seed / synthetic inputs, the Rcpp result matches the old R loop to within a tight tolerance (e.g., relative error <1e-10) for several random z_disp vectors and different grid settings.

Acceptance criteria:

• The package builds and checks.
• For representative inputs, disp_stars=TRUE outputs are numerically consistent with the previous implementation (tolerance-based).
• The new code meaningfully reduces runtime of the dispersion block.

No other functional changes.

This pull request was created from Copilot chat.

---

💡 You can make Copilot smarter by setting up custom instructions, customizing its development environment and configuring Model Context Protocol (MCP) servers. Learn more Copilot coding agent tips in the docs.

[asgr]

Speeds up the star dispersion