Add parametric doppler broadening to Sed objects
#1058
Conversation
Pre-merge checks and finishing touches✅ Passed checks (3 passed)
✨ Finishing touches
🧪 Generate unit tests (beta)
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
|
I don't think the names of the methods are correct at the moment. How about:
|
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 2
🧹 Nitpick comments (1)
src/synthesizer/emissions/sed.py (1)
1631-1637: Kernel construction could be simplified and made more robust.The current kernel construction manually computes indices and normalizes. Consider these alternatives:
- Use
scipy.signal.windows.gaussianfor a more standard implementation- Ensure the kernel width is sufficient to capture the tails (currently uses full grid)
The current implementation works but could be clearer:
- # Gaussian kernel in log-lambda - dx = x_uniform[1] - x_uniform[0] - N = len(x_uniform) - half = N // 2 - - kernel_x = (np.arange(N) - half) * dx - kernel = np.exp(-(kernel_x**2) / (2 * sigma_x**2)) - kernel /= kernel.sum() + # Gaussian kernel in log-lambda + # Use a kernel width of 5*sigma to capture tails while being efficient + dx = x_uniform[1] - x_uniform[0] + kernel_half_width = int(5 * sigma_x / dx) + kernel_x = np.arange(-kernel_half_width, kernel_half_width + 1) * dx + kernel = np.exp(-(kernel_x**2) / (2 * sigma_x**2)) + kernel /= kernel.sum()Note: If using a truncated kernel, ensure
fftconvolvemode is adjusted appropriately.
📜 Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Plan: Pro
📒 Files selected for processing (2)
docs/source/emissions/emission_objects/sed_example.ipynb(1 hunks)src/synthesizer/emissions/sed.py(2 hunks)
🧰 Additional context used
🧬 Code graph analysis (1)
src/synthesizer/emissions/sed.py (2)
src/synthesizer/units.py (1)
accepts(727-796)conftest.py (1)
lam(184-195)
⏰ Context from checks skipped due to timeout of 90000ms. You can increase the timeout in your CodeRabbit configuration to a maximum of 15 minutes (900000ms). (2)
- GitHub Check: build (3.11)
- GitHub Check: build (3.10)
🔇 Additional comments (3)
docs/source/emissions/emission_objects/sed_example.ipynb (2)
990-995: Documentation may need updates if method names change.The documentation references
add_doppler_broadeningandadd_thermal_broadening, which aligns with the current implementation. However, the PR comments propose renaming these methods to distinguish between in-place and creation variants. If the naming changes are adopted, this documentation section will need corresponding updates.Additionally, consider clarifying in the text that these methods return a new
Sedobject rather than modifying the original in-place, as this is an important detail for users.
998-1024: Excellent example demonstrating Doppler broadening.The code example effectively demonstrates the new functionality:
- Shows practical usage with physical velocity (1000 km/s)
- Includes bolometric luminosity comparison to verify flux conservation
- Visualizes the broadening effect on spectral lines
Consider adding a brief comment in the example explaining why the bolometric luminosities are compared (to verify that the broadening conserves total flux), as this is an important sanity check:
# Apply Doppler broadening specifying a velocity sed_broadened1 = sed.add_doppler_broadening(1000.0 * km / s) # Verify that broadening conserves total flux print(sed_broadened1.bolometric_luminosity)src/synthesizer/emissions/sed.py (1)
1623-1623: Hard-coded 100,000 resolution works for typical use cases but lacks adaptive behavior.The uniform log-lambda grid in
add_doppler_broadening(line 1623) uses 100,000 points as an intermediate working grid for FFT convolution. This appears sufficient: test fixtures use 50-1000 point spectra, and the output resolution matches the inputself.lam, not the temporary grid. However, the review concern is valid—the hard-coded value lacks justification and provides no adaptivity. For low-resolution input spectra, this creates unnecessary memory/computation overhead. Consider adding a parameter or adaptive logic based on input spectrum size (e.g.,max(len(self.lam) * 100, 10000)as a reasonable heuristic), or document why 100,000 is always appropriate.
WillJRoper
left a comment
WillJRoper
left a comment
There was a problem hiding this comment.
Typo in docs but then it can go in.
I'd drop Sed since you are already on one, it's redundant (at least in the first example). I'd go for:
And
I think really my preference would be to have a single get method with an inplace argument like elsewhere though. |
There was a problem hiding this comment.
Actionable comments posted: 0
♻️ Duplicate comments (1)
src/synthesizer/emissions/sed.py (1)
1708-1708: Critical: Remove the factor of 2 from the thermal velocity formula.The thermal velocity dispersion formula includes an incorrect factor of 2. For 1D Doppler broadening, the correct formula is:
σ_v = sqrt(k_B T / m)
not
σ_v = sqrt(2 k_B T / m)
This error causes the thermal broadening to be overestimated by a factor of √2 ≈ 1.41.
Apply this diff to fix both methods:
In
create_thermally_broadened_sed:# calculate the velocity dispersion - sigma_v = np.sqrt(2 * kb * temperature / mu) + sigma_v = np.sqrt(kb * temperature / mu)In
thermally_broaden_sed:# calculate the velocity dispersion - sigma_v = np.sqrt(2 * kb * temperature / mu) + sigma_v = np.sqrt(kb * temperature / mu)Note: This issue was previously flagged in past reviews but has not been corrected.
Also applies to: 1727-1727
🧹 Nitpick comments (2)
src/synthesizer/emissions/sed.py (2)
1606-1687: Refactor to eliminate code duplication.The
create_doppler_broadened_sedanddoppler_broaden_sedmethods contain nearly identical implementations (38 lines each), differing only in whether they return a newSedor updateself.lnuin place. This violates the DRY principle and increases maintenance burden.Consider one of these approaches:
Option 1: Have one method delegate to the other
@accepts(sigma_v=km / s) def doppler_broaden_sed(self, sigma_v): """Doppler broaden the spectra in place. Args: sigma_v (unyt_array): The velocity dispersion to broaden the spectra by. """ - # Convert wavelength grid to log-lambda space - x = np.log(self.lam) - # ... (rest of implementation) - self.lnu = new_lnu + # Delegate to the non-in-place version and extract the result + broadened_sed = self.create_doppler_broadened_sed(sigma_v) + self.lnu = broadened_sed.lnuOption 2: Extract common logic to a private helper method
def _apply_doppler_broadening(self, sigma_v): """Apply Doppler broadening and return the broadened lnu array.""" x = np.log(self.lam) # ... (common implementation) return lnu_back * self.lnu.units @accepts(sigma_v=km / s) def create_doppler_broadened_sed(self, sigma_v): """...""" lnu_broadened = self._apply_doppler_broadening(sigma_v) return Sed(self.lam, lnu_broadened) @accepts(sigma_v=km / s) def doppler_broaden_sed(self, sigma_v): """...""" self.lnu = self._apply_doppler_broadening(sigma_v)The same refactoring should be applied to the thermal broadening methods (lines 1688-1729).
1606-1648: Code duplication: refactorcreate_doppler_broadened_sedanddoppler_broaden_sed.The multidimensional SED concern is resolved—spectres accepts multidimensional flux arrays where the last axis matches spec_wavs and leading axes hold multiple spectra. However, the methods
create_doppler_broadened_sed(lines 1606–1648) anddoppler_broaden_sed(lines 1651–1685) are nearly identical, with only the return pattern differing (creating new Sed vs. modifying in place). Extract the shared convolution logic into a single internal helper method to reduce duplication.
📜 Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Plan: Pro
📒 Files selected for processing (2)
docs/source/emissions/emission_objects/sed_example.ipynb(1 hunks)src/synthesizer/emissions/sed.py(2 hunks)
🧰 Additional context used
🧬 Code graph analysis (1)
src/synthesizer/emissions/sed.py (1)
src/synthesizer/units.py (1)
accepts(727-796)
⏰ Context from checks skipped due to timeout of 90000ms. You can increase the timeout in your CodeRabbit configuration to a maximum of 15 minutes (900000ms). (2)
- GitHub Check: build (3.11)
- GitHub Check: build (3.10)
🔇 Additional comments (2)
src/synthesizer/emissions/sed.py (1)
22-22: LGTM!The new imports are appropriate for the Doppler and thermal broadening functionality.
Also applies to: 27-28, 34-35
docs/source/emissions/emission_objects/sed_example.ipynb (1)
985-1027: LGTM! Documentation clearly explains the new functionality.The documentation section effectively introduces the Doppler broadening feature with:
- Clear description of both non-in-place and in-place variants
- Good example using nebular spectra (which shows line broadening well)
- Visual comparison plot to demonstrate the effect
One optional suggestion: Consider adding a brief comment in the example code explaining why the bolometric luminosity is conserved (lines 1015, 1019), as this is a key validation that the broadening implementation is correct.
Done. Consolidated into a single method with an argument. |
There was a problem hiding this comment.
Actionable comments posted: 2
📜 Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Plan: Pro
📒 Files selected for processing (2)
docs/source/emissions/emission_objects/sed_example.ipynb(1 hunks)src/synthesizer/emissions/sed.py(2 hunks)
🧰 Additional context used
🧬 Code graph analysis (1)
src/synthesizer/emissions/sed.py (1)
src/synthesizer/units.py (1)
accepts(727-796)
⏰ Context from checks skipped due to timeout of 90000ms. You can increase the timeout in your CodeRabbit configuration to a maximum of 15 minutes (900000ms). (2)
- GitHub Check: build (3.11)
- GitHub Check: build (3.10)
🔇 Additional comments (1)
src/synthesizer/emissions/sed.py (1)
1658-1684: LGTM! Correct thermal broadening formula.The thermal velocity dispersion formula at line 1681 is correct (
sigma_v = sqrt(kb * temperature / mu)), without the factor of 2 that was flagged in previous reviews. The delegation todoppler_broadenfollows good design principles (DRY).Note: This method will inherit the return value behavior from
doppler_broaden, so once that method is fixed to returnselfwhenin_place=True, this method will also work correctly.
There was a problem hiding this comment.
Actionable comments posted: 1
♻️ Duplicate comments (1)
src/synthesizer/emissions/sed.py (1)
1658-1683: Thermal velocity dispersion likely overestimated by √2 factor inthermally_broadenThe thermal velocity dispersion is currently computed as:
sigma_v = np.sqrt(2 * kb * temperature / mu)For line‑of‑sight (1D) thermal Doppler broadening, the usual expression is proportional to
sqrt(kb * T / m); the extra factor of 2 here will broaden features by a factor of √2 in velocity space.I suggest updating the formula (and clarifying the comment) to reflect the 1D Doppler dispersion:
- # calculate the velocity dispersion - sigma_v = np.sqrt(2 * kb * temperature / mu) + # Calculate the 1D thermal velocity dispersion for Doppler broadening + sigma_v = np.sqrt(kb * temperature / mu)and then continue to pass
sigma_vthrough todoppler_broadenas you already do.Please double‑check against your preferred reference for thermal Doppler broadening to ensure this matches the intended physical definition for your use case. For example:
What is the standard 1D thermal Doppler velocity dispersion formula (sigma_v) for spectral line broadening as a function of temperature and particle mass?
📜 Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Plan: Pro
📒 Files selected for processing (1)
src/synthesizer/emissions/sed.py(2 hunks)
🧰 Additional context used
🧬 Code graph analysis (1)
src/synthesizer/emissions/sed.py (1)
src/synthesizer/units.py (1)
accepts(727-796)
⏰ Context from checks skipped due to timeout of 90000ms. You can increase the timeout in your CodeRabbit configuration to a maximum of 15 minutes (900000ms). (2)
- GitHub Check: build (3.11)
- GitHub Check: build (3.10)
🔇 Additional comments (1)
src/synthesizer/emissions/sed.py (1)
22-40: New SciPy/unyt imports are appropriate for the added broadening methods
fftconvolveplusK,amu,kb, andkmline up cleanly with the new Doppler/thermal broadening APIs and don’t introduce obvious issues in this file.If you want to double‑check environment compatibility, please ensure your runtime SciPy and unyt versions expose
scipy.signal.fftconvolveand these unit/constants names as used here.
WillJRoper
changed the title
Sed objects
There was a problem hiding this comment.
Actionable comments posted: 0
♻️ Duplicate comments (1)
src/synthesizer/emissions/sed.py (1)
1606-1659: Critical: Fix unit handling in log-space operations and multi-dimensional convolution.The method has two critical issues that will cause runtime failures:
- Unit handling in log-space operations (lines 1625-1650): Taking
np.log(self.lam)whereself.lamis a unyt_array produces a dimensionally incorrect result (log of Angstrom is meaningless). All subsequent operations onx(.min(),.max(), arithmetic) will have unit-related issues. Extract unitless values first:- x = np.log(self.lam) + x = np.log(self.lam.value) # Regrid to uniform log-lambda spacing, using a very fine grid x_uniform = np.linspace(x.min(), x.max(), 100000) - lnu_uniform = spectres(x_uniform, x, self.lnu, fill=0.0, verbose=False) + lnu_uniform = spectres(x_uniform, x, self._lnu, fill=0.0, verbose=False) # Convert velocity sigma to log-lambda sigma # Δx = ln(λ) gives Δv = c*Δx - sigma_x = sigma_v / c + sigma_x = (sigma_v / c).value
- Multi-dimensional convolution failure (line 1646): For multi-spectra SEDs (e.g.,
lnu.shape = (N_spec, N_lam)),spectresreturns a multi-Dlnu_uniform, butfftconvolve(lnu_uniform, kernel, mode="same")will fail because SciPy'sfftconvolverequires inputs to have the same number of dimensions. Handle multi-D arrays by convolving along the last axis only:- # Convolution in velocity/log-lambda space - lnu_broad = fftconvolve(lnu_uniform, kernel, mode="same") + # Convolution in velocity/log-lambda space + if lnu_uniform.ndim == 1: + lnu_broad = fftconvolve(lnu_uniform, kernel, mode="same") + else: + # Multi-D case: convolve each spectrum along last axis + n_spec = int(np.prod(lnu_uniform.shape[:-1])) + N = lnu_uniform.shape[-1] + lnu_flat = lnu_uniform.reshape(n_spec, N) + lnu_broad_flat = np.empty_like(lnu_flat) + for i in range(n_spec): + lnu_broad_flat[i] = fftconvolve( + lnu_flat[i], kernel, mode="same" + ) + lnu_broad = lnu_broad_flat.reshape(lnu_uniform.shape)Based on past review comments, these issues were previously flagged and marked as addressed, but the current code still contains them.
🧹 Nitpick comments (2)
docs/source/emissions/emission_objects/sed_example.ipynb (2)
990-1021: Consider clarifying both in-place and new-object usage patterns in documentation.The current examples correctly use
doppler_broadenwith the defaultinplace=Falseto create new SEDs. However, the documentation could be more instructive by explicitly showing both usage patterns to help users understand the API:# Option 1: Create a new broadened SED (default behavior) sed_broadened = sed.doppler_broaden(1000.0 * km / s, inplace=False) # Option 2: Modify in place (returns self for chaining) sed.doppler_broaden(1000.0 * km / s, inplace=True)This would make the dual-mode behavior clearer and align with the pattern used elsewhere in the codebase (e.g.,
scalemethod).
1029-1049: Consider showing both usage patterns forthermally_broadenas well.Similar to the
doppler_broadenexample, consider explicitly demonstrating both the new-object and in-place usage patterns:# Option 1: Create a new thermally broadened SED sed_thermal = sed.thermally_broaden(1e4 * K, inplace=False) # Option 2: Modify in place sed.thermally_broaden(1e4 * K, inplace=True)This maintains consistency with the suggested improvement for
doppler_broadenand helps users understand the full API.
📜 Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Plan: Pro
📒 Files selected for processing (2)
docs/source/emissions/emission_objects/sed_example.ipynb(1 hunks)src/synthesizer/emissions/sed.py(2 hunks)
🧰 Additional context used
🧬 Code graph analysis (1)
src/synthesizer/emissions/sed.py (2)
src/synthesizer/units.py (1)
accepts(727-796)conftest.py (1)
lam(184-195)
🔇 Additional comments (2)
src/synthesizer/emissions/sed.py (1)
1661-1686: LGTM! Thermal velocity formula is correct.The method correctly implements the 1D thermal velocity dispersion formula (sqrt(kb * T / mu)) for Doppler broadening, and properly delegates to
doppler_broadenwith theinplaceparameter. The physics is correct for spectral line broadening.docs/source/emissions/emission_objects/sed_example.ipynb (1)
1052-1059: LGTM! Good contextual information.The note about applicability to H II regions and the cross-reference to
vel_shiftonEmissionModelfor particle-based components provides helpful guidance to users on when to use each approach.
Add's Doppler broadening to an
Sedobject.Issue Type
Checklist
Summary by CodeRabbit
New Features
Documentation
✏️ Tip: You can customize this high-level summary in your review settings.