From 738deabe64cb9a6a52bd34f19b923b7343753383 Mon Sep 17 00:00:00 2001 From: dnarayanan Date: Wed, 8 Jul 2026 15:53:08 -0400 Subject: [PATCH] grid_integrate: deposit the specific energy spectrum at interaction points Only the scalar specific energy received the partial-path deposit when a photon interacted inside a cell; the frequency-resolved spectrum (specific_energy_sum_spectrum) was only deposited at wall crossings. The binned spectrum was therefore biased low in any cell optically thick enough for photons to interact within a single crossing -- measured at 20% low for tau_V ~ 0.1-1 up to ~7x low for tau_V > 10 on a galaxy simulation snapshot, while the bolometric (scalar) energy was correct. Mirror the wall-crossing spectrum deposit at the interaction point with the partial path length tact, so both arrays receive identical increments and the spectrum's bin-sum closes on the scalar exactly (verified closure = 1.0000 in every cell of an optically thick test model with the equivalent fix on the pre-rewrite branch, vs 0.12-0.97 without it). Co-Authored-By: Claude Fable 5 --- src/grid/grid_propagate_3d.f90 | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/grid/grid_propagate_3d.f90 b/src/grid/grid_propagate_3d.f90 index 3769704c..41549b4b 100644 --- a/src/grid/grid_propagate_3d.f90 +++ b/src/grid/grid_propagate_3d.f90 @@ -195,6 +195,23 @@ subroutine grid_integrate(p,tau_required,tau_achieved) end if end do + ! Deposit the same energy in the frequency-resolved spectrum. + ! Without this, the energy deposited over the partial path to + ! each interaction point is counted in the scalar specific + ! energy but missing from the binned spectrum, which biases + ! the spectrum low in any cell optically thick enough for + ! photons to interact within a single crossing. + if (compute_specific_energy_spectrum) then + idx = minloc(abs(log_nu_bins - log10(p%nu)), DIM=1) + do id=1,n_dust + if(density(p%icell%ic, id) > 0._dp) then + specific_energy_sum_spectrum(p%icell%ic, id, idx) = & + & specific_energy_sum_spectrum(p%icell%ic, id, idx) & + & + tact * p%current_kappa(id) * p%energy + end if + end do + end if + if(debug) write(*,'(" [debug] end grid_integrate")') return