|
| 1 | +#include <Rcpp.h> |
| 2 | +#include <cmath> |
| 3 | +using namespace Rcpp; |
| 4 | + |
| 5 | +// disp_stars_cpp: Rcpp implementation of the stellar dispersion convolution. |
| 6 | +// |
| 7 | +// Equivalent to the R loop in SFHfunc()/SFHburst() when disp_stars=TRUE: |
| 8 | +// for(i in seq_along(grid)){ |
| 9 | +// z_seq_log = log10(1 + grid[i]*z_disp) |
| 10 | +// new_wave_log = wave_log + z_seq_log |
| 11 | +// new_lum = lum_log - z_seq_log |
| 12 | +// new_lum = 10^approx(x=new_wave_log, y=new_lum, xout=wave_log, |
| 13 | +// rule=2, yleft=new_lum[1], yright=new_lum[n])$y |
| 14 | +// lum_conv = lum_conv + weights[i] * new_lum |
| 15 | +// } |
| 16 | +// return(lum_conv * res) |
| 17 | +// |
| 18 | +// Inputs: |
| 19 | +// wave_log - log10 wavelengths (monotone increasing, length n) |
| 20 | +// lum_log - log10 luminosity (length n) |
| 21 | +// z_disp - wavelength-dependent velocity dispersion in z units (length n) |
| 22 | +// grid - Gaussian quadrature grid points, e.g. seq(-range, range, by=res) |
| 23 | +// weights - dnorm(grid) |
| 24 | +// res - grid spacing (multiplied into final result) |
| 25 | +// |
| 26 | +// [[Rcpp::export(".disp_stars_cpp")]] |
| 27 | +NumericVector disp_stars_cpp(NumericVector wave_log, |
| 28 | + NumericVector lum_log, |
| 29 | + NumericVector z_disp, |
| 30 | + NumericVector grid, |
| 31 | + NumericVector weights, |
| 32 | + double res) { |
| 33 | + int n = wave_log.size(); |
| 34 | + int ng = grid.size(); |
| 35 | + |
| 36 | + if (lum_log.size() != n) stop("lum_log must have same length as wave_log"); |
| 37 | + if (z_disp.size() != n) stop("z_disp must have same length as wave_log"); |
| 38 | + if (weights.size() != ng) stop("weights must have same length as grid"); |
| 39 | + |
| 40 | + NumericVector out(n, 0.0); |
| 41 | + // Reusable buffers for the shifted wavelength / luminosity grids |
| 42 | + NumericVector x_src(n); |
| 43 | + NumericVector y_src(n); |
| 44 | + |
| 45 | + for (int gi = 0; gi < ng; gi++) { |
| 46 | + double g = grid[gi]; |
| 47 | + double w = weights[gi]; |
| 48 | + |
| 49 | + // Build shifted grids: x_src = wave_log + log10(1 + g*z_disp) |
| 50 | + // y_src = lum_log - log10(1 + g*z_disp) |
| 51 | + for (int j = 0; j < n; j++) { |
| 52 | + double val = 1.0 + g * z_disp[j]; |
| 53 | + if (val <= 0.0) { |
| 54 | + stop("1 + grid[i]*z_disp <= 0: cannot take log10. " |
| 55 | + "Reduce veldisp or the grid range."); |
| 56 | + } |
| 57 | + double lv = std::log10(val); |
| 58 | + x_src[j] = wave_log[j] + lv; |
| 59 | + y_src[j] = lum_log[j] - lv; |
| 60 | + } |
| 61 | + |
| 62 | + // Linear interpolation of y_src(x_src) onto wave_log. |
| 63 | + // Uses rule=2: extrapolate with boundary values y_src[0] / y_src[n-1]. |
| 64 | + // Two-pointer sweep works because wave_log is monotone increasing and |
| 65 | + // x_src remains monotone for the small shifts typical of veldisp. |
| 66 | + // lo is intentionally reset to 0 for each grid point gi: since wave_log |
| 67 | + // is scanned left-to-right each iteration, lo must start fresh. |
| 68 | + int lo = 0; |
| 69 | + for (int j = 0; j < n; j++) { |
| 70 | + double xq = wave_log[j]; |
| 71 | + double y_interp; |
| 72 | + |
| 73 | + if (xq <= x_src[0]) { |
| 74 | + y_interp = y_src[0]; |
| 75 | + } else if (xq >= x_src[n - 1]) { |
| 76 | + y_interp = y_src[n - 1]; |
| 77 | + } else { |
| 78 | + // Advance lo so that x_src[lo] <= xq < x_src[lo+1] |
| 79 | + while (lo < n - 2 && x_src[lo + 1] <= xq) { |
| 80 | + lo++; |
| 81 | + } |
| 82 | + double t = (xq - x_src[lo]) / (x_src[lo + 1] - x_src[lo]); |
| 83 | + y_interp = y_src[lo] + t * (y_src[lo + 1] - y_src[lo]); |
| 84 | + } |
| 85 | + |
| 86 | + out[j] += w * std::pow(10.0, y_interp); |
| 87 | + } |
| 88 | + } |
| 89 | + |
| 90 | + // Scale by the grid spacing |
| 91 | + for (int j = 0; j < n; j++) { |
| 92 | + out[j] *= res; |
| 93 | + } |
| 94 | + |
| 95 | + return out; |
| 96 | +} |
0 commit comments