From 096ff1443374c03da260437e1e587f0fb8280220 Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Fri, 25 Sep 2026 22:30:34 -0700 Subject: [PATCH 01/12] Add failing test for Halo.fake_image --- tests/test_galhalo.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/tests/test_galhalo.py b/tests/test_galhalo.py index 45858bd..8e85375 100644 --- a/tests/test_galhalo.py +++ b/tests/test_galhalo.py @@ -167,3 +167,27 @@ def test_halo_io(test): assert np.all(percent_diff(new_halo.taux,test.taux) < 0.01) assert np.all(percent_diff(new_halo.norm_int.flatten(),test.norm_int.flatten()) < 0.01) assert new_halo.lam.unit == 'keV' + +def create_arf(tmp_path, outfile, evals, effarea) -> str: + """ + Create a dummy ARF file for testing + """ + # import numpy as np + from astropy.io import fits + from astropy.table import Table + + t = Table() + t['ENERG_LO'] = evals[:-1] + t['ENERG_HI'] = evals[1:] + t['SPECRESP'] = np.full(len(evals) - 1, effarea) + + hdu_table = fits.BinTableHDU.from_columns(t.as_array(), name='SPECRESP') + hdu_primary = fits.PrimaryHDU() + hdul = fits.HDUList([hdu_primary, hdu_table]) + hdul.writeto(tmp_path / outfile, overwrite=True) + return str(tmp_path / outfile) + +def test_fake_image(tmp_path): + new_halo = Halo(EVALS, THVALS) + arf = create_arf(tmp_path, 'test_arf.fits', EVALS, 100.0) + new_halo.fake_image(arf, src_flux=FABS, exposure=1e4, pix_scale=1.0, num_pix=[32, 32]) \ No newline at end of file From 4c388b457b97f7f5a4d405a96c43ac403884344d Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Sat, 26 Sep 2026 00:04:36 -0700 Subject: [PATCH 02/12] Fix unit handling in Halo.fake_image fake_image raised AttributeError on every call: it branched on self.lam_unit, an attribute removed in 949a0b6 when Halo migrated to astropy Quantities. Two further errors were masked behind it: self.lam was multiplied by u.angstrom despite already carrying a unit, and lmin/lmax were compared against unitless floats. self.lam carries its own unit, so the unit-string branching is unnecessary. Convert explicitly with the spectral equivalency instead, which handles any input unit rather than just two. Adds the first test coverage for fake_image. --- src/xdust/halos/halo.py | 13 ++++--------- tests/test_galhalo.py | 32 ++++++++++++++++++++++++-------- 2 files changed, 28 insertions(+), 17 deletions(-) diff --git a/src/xdust/halos/halo.py b/src/xdust/halos/halo.py index 11bbe74..8bc910d 100644 --- a/src/xdust/halos/halo.py +++ b/src/xdust/halos/halo.py @@ -403,23 +403,18 @@ def fake_image(self, arf, src_flux, exposure, arf = InterpolatedUnivariateSpline(arf_x, arf_y, k=1) # Source counts to use for each energy bin - if self.lam_unit == 'angs': - ltemp = self.lam * u.angstrom - ltemp_kev = ltemp.to(u.keV, equivalencies=u.spectral()).value - arf_temp = arf(ltemp_kev)[::-1] - src_counts = src_flux * arf_temp * exposure - else: - src_counts = src_flux * arf(self.lam) * exposure + lam_kev = self.lam.to(u.keV, equivalencies=u.spectral()).value + src_counts = src_flux * arf(lam_kev) * exposure # Decide which energy indexes to use if lmin is None: imin = 0 else: - imin = min(np.arange(len(self.lam))[self.lam >= lmin]) + imin = min(np.arange(len(self.lam))[self.lam >= lmin * self.lam.unit]) if lmax is None: iend = len(self.lam) else: - iend = max(np.arange(len(self.lam))[self.lam <= lmax]) + iend = max(np.arange(len(self.lam))[self.lam <= lmax * self.lam.unit]) #iend = imax #if imax < 0: diff --git a/tests/test_galhalo.py b/tests/test_galhalo.py index 8e85375..36094c7 100644 --- a/tests/test_galhalo.py +++ b/tests/test_galhalo.py @@ -2,6 +2,8 @@ import numpy as np from scipy.integrate import trapezoid as trapz import astropy.units as u +from astropy.io import fits +from astropy.table import Table from xdust.halos import * from xdust import grainpop @@ -172,22 +174,36 @@ def create_arf(tmp_path, outfile, evals, effarea) -> str: """ Create a dummy ARF file for testing """ - # import numpy as np - from astropy.io import fits - from astropy.table import Table - t = Table() t['ENERG_LO'] = evals[:-1] t['ENERG_HI'] = evals[1:] t['SPECRESP'] = np.full(len(evals) - 1, effarea) - hdu_table = fits.BinTableHDU.from_columns(t.as_array(), name='SPECRESP') + hdu_table = fits.BinTableHDU(t.as_array(), name='SPECRESP') hdu_primary = fits.PrimaryHDU() hdul = fits.HDUList([hdu_primary, hdu_table]) hdul.writeto(tmp_path / outfile, overwrite=True) return str(tmp_path / outfile) -def test_fake_image(tmp_path): - new_halo = Halo(EVALS, THVALS) +@pytest.fixture +def halo(): + new_halo = galhalo.UniformGalHalo(EVALS, THVALS) + new_halo.calculate(GPOP) + return new_halo + +@pytest.fixture +def arf(tmp_path): arf = create_arf(tmp_path, 'test_arf.fits', EVALS, 100.0) - new_halo.fake_image(arf, src_flux=FABS, exposure=1e4, pix_scale=1.0, num_pix=[32, 32]) \ No newline at end of file + return arf + +def test_fake_image(halo, arf): + image = halo.fake_image(arf, src_flux=FABS, exposure=1e4, pix_scale=1.0, num_pix=[8, 16]) + assert image.shape == (16, 8) + assert np.all(image >= 0.0) + assert np.all(image == np.floor(image)) + assert image.sum() > 0 + +def test_fake_image_lmin_lmax(halo, arf): + restricted = halo.fake_image(arf, src_flux=FABS, exposure=1e4, pix_scale=1.0, num_pix=[8, 16], lmin=0.5, lmax=5.0) + unrestricted = halo.fake_image(arf, src_flux=FABS, exposure=1e4, pix_scale=1.0, num_pix=[8, 16]) + assert 0 < restricted.sum() < unrestricted.sum() From a0d3630255f68cf73ccda61c59bf01d86ae56ac0 Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Mon, 28 Sep 2026 15:18:38 -0700 Subject: [PATCH 03/12] Add failing test for UnitConversionError: 's / m' and 's' (time) are not convertible --- tests/test_galhalo.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/test_galhalo.py b/tests/test_galhalo.py index 36094c7..d68311a 100644 --- a/tests/test_galhalo.py +++ b/tests/test_galhalo.py @@ -207,3 +207,10 @@ def test_fake_image_lmin_lmax(halo, arf): restricted = halo.fake_image(arf, src_flux=FABS, exposure=1e4, pix_scale=1.0, num_pix=[8, 16], lmin=0.5, lmax=5.0) unrestricted = halo.fake_image(arf, src_flux=FABS, exposure=1e4, pix_scale=1.0, num_pix=[8, 16]) assert 0 < restricted.sum() < unrestricted.sum() + +def test_time_delay(): + result1 = galhalo.time_delay(100.0, 0.5, 8.0) + assert result1 == pytest.approx(96770, rel=1e-3) + + result2 = galhalo.time_delay(200.0, 0.5, 8.0) + assert result2 == pytest.approx(4 * result1, rel=1e-3) From 04737a730ba0c00122edb1c55ba46559f97006d2 Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Mon, 28 Sep 2026 15:20:57 -0700 Subject: [PATCH 04/12] Fix unit handling in time_delay --- src/xdust/halos/galhalo.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/xdust/halos/galhalo.py b/src/xdust/halos/galhalo.py index a091d02..5823f94 100644 --- a/src/xdust/halos/galhalo.py +++ b/src/xdust/halos/galhalo.py @@ -601,5 +601,4 @@ def time_delay(alpha, x, dkpc): Assumes small angle scattering. """ delta_x = path_diff(alpha, x) - d_cm = dkpc * 1.e3 * u.pc.to('cm') # cm - return delta_x * (d_cm / c.c).to('s').value # seconds + return delta_x * (dkpc * u.kpc / c.c).to('s').value From bfaaa8afc69420f819a170b245f5823bfbe2e8ed Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Mon, 28 Sep 2026 15:36:21 -0700 Subject: [PATCH 05/12] Add a failing test for variable_profile for UnitConversionError --- tests/test_galhalo.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/test_galhalo.py b/tests/test_galhalo.py index d68311a..06ba5cf 100644 --- a/tests/test_galhalo.py +++ b/tests/test_galhalo.py @@ -191,6 +191,13 @@ def halo(): new_halo.calculate(GPOP) return new_halo +@pytest.fixture +def screen_halo(): + new_halo = galhalo.ScreenGalHalo(EVALS, THVALS) + new_halo.calculate(GPOP, x=0.5) + new_halo.calculate_intensity(FABS) + return new_halo + @pytest.fixture def arf(tmp_path): arf = create_arf(tmp_path, 'test_arf.fits', EVALS, 100.0) @@ -214,3 +221,13 @@ def test_time_delay(): result2 = galhalo.time_delay(200.0, 0.5, 8.0) assert result2 == pytest.approx(4 * result1, rel=1e-3) + +def test_variable_profile(screen_halo): + time = np.linspace(0, 100, 10) + lc = np.ones_like(time) + dist = 8.0 + var_profile = screen_halo.variable_profile(time, lc, dist=dist) + assert var_profile.shape == (len(EVALS), len(THVALS)) + assert np.sum(var_profile) > 0 + # 2 * lc and assert the result is ~2× + From bfffa2a0da8011dc05acf0ff784a7fa7b2f2489b Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Mon, 28 Sep 2026 16:08:52 -0700 Subject: [PATCH 06/12] Fix unit handling in variable_profile time_delay documents alpha as a bare float in arcsec, but variable_profile passed self.theta, a Quantity, which broke the comparison in _is_small_angle. Convert at the call site. The intensity accumulator was a plain ndarray while the values written into it carry 1/arcsec^2, so give it norm_int's unit. Return type is now a Quantity; docstring updated to match. --- src/xdust/halos/galhalo.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/xdust/halos/galhalo.py b/src/xdust/halos/galhalo.py index 5823f94..95848c4 100644 --- a/src/xdust/halos/galhalo.py +++ b/src/xdust/halos/galhalo.py @@ -153,7 +153,7 @@ def variable_profile(self, time, lc, dist=8.0, tnow=None): Returns ------- - numpy.ndarray (NE x NTH) [fabs units / arcsec**2] + astropy.units.Quantity (NE x NTH) [fabs units / arcsec**2] Scattering halo intensity as a function of energy and observation angle. """ @@ -171,11 +171,11 @@ def variable_profile(self, time, lc, dist=8.0, tnow=None): theta_rad = self.theta/3600./180.*np.pi ne, ntheta = len(self.lam), len(self.theta) - inten = np.zeros(shape=(ne, ntheta)) + inten = np.zeros(shape=(ne, ntheta)) * self.norm_int.unit lctm = (time-tzero) for i in range(len(self.lam)): - deltat = time_delay(self.theta, self.x, dist) * u.second.to(u.day) + deltat = time_delay(self.theta.to(u.arcsec).value, self.x, dist) * u.second.to(u.day) t = tnow - deltat for j in range(ntheta): inten[i,j] += np.interp(t[j], time, lc * self.norm_int[i,j] * self.fabs[i]) From 0786ca0c053233ff703214d5792db79b68cff651 Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Mon, 28 Sep 2026 16:34:48 -0700 Subject: [PATCH 07/12] Add failing test for fake_variable_image --- tests/test_galhalo.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/tests/test_galhalo.py b/tests/test_galhalo.py index 06ba5cf..185067b 100644 --- a/tests/test_galhalo.py +++ b/tests/test_galhalo.py @@ -226,8 +226,20 @@ def test_variable_profile(screen_halo): time = np.linspace(0, 100, 10) lc = np.ones_like(time) dist = 8.0 - var_profile = screen_halo.variable_profile(time, lc, dist=dist) - assert var_profile.shape == (len(EVALS), len(THVALS)) - assert np.sum(var_profile) > 0 + var_profile1 = screen_halo.variable_profile(time, lc, dist=dist) + assert var_profile1.shape == (len(EVALS), len(THVALS)) + assert np.sum(var_profile1) > 0 + + var_profile2 = screen_halo.variable_profile(time, 2 * lc, dist=dist) + assert var_profile2.shape == (len(EVALS), len(THVALS)) + assert np.sum(var_profile2) > 0 # 2 * lc and assert the result is ~2× +def test_fake_variable_image(screen_halo, arf): + time = np.linspace(0, 100, 10) + lc = np.ones_like(time) + image = screen_halo.fake_variable_image(time, lc=lc, arf=arf, exposure=1e4, num_pix=[8, 16]) + assert image.shape == (16, 8) + assert np.all(image >= 0.0) + assert np.all(image == np.floor(image)) + assert image.sum() > 0 From 6fac7d87dd48be865edc26ba0fc2ed6d206518dd Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Mon, 28 Sep 2026 16:38:57 -0700 Subject: [PATCH 08/12] Restore var_profile computation in fake_variable_image Commit fe4c562 ("first try of remaking docstrings with claude") removed the line that computes var_profile, leaving the loop below referencing an undefined name. fake_variable_image raised NameError on every call. variable_profile already defaults tnow to time[-1], so the separate time_now block that commit also removed is unnecessary. --- src/xdust/halos/galhalo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/xdust/halos/galhalo.py b/src/xdust/halos/galhalo.py index 95848c4..10f1a21 100644 --- a/src/xdust/halos/galhalo.py +++ b/src/xdust/halos/galhalo.py @@ -237,7 +237,7 @@ def fake_variable_image(self, time, lc, arf, to simulate the number of counts in each pixel. """ # intensity cube (NE x NTH), phot/cm^2/s/arcsec^2 - + var_profile = self.variable_profile(time, lc, dist=dist, tnow=tnow) # Decide which energy indexes to use if lmin is None: imin = 0 From 60d5a218ada6b18f725decf33080f85c5dd9d091 Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Mon, 28 Sep 2026 16:39:57 -0700 Subject: [PATCH 09/12] Fix unit handling in fake_variable_image Same lam_unit branching as Halo.fake_image: the attribute was removed in 949a0b6 during the astropy Quantity migration, but the read survived here. self.lam carries its own unit, so convert with the spectral equivalency and strip to a plain array for the spline. --- src/xdust/halos/galhalo.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/xdust/halos/galhalo.py b/src/xdust/halos/galhalo.py index 10f1a21..973813c 100644 --- a/src/xdust/halos/galhalo.py +++ b/src/xdust/halos/galhalo.py @@ -261,13 +261,7 @@ def fake_variable_image(self, time, lc, arf, arf = InterpolatedUnivariateSpline(arf_x, arf_y, k=1) # Conversion erg -> ct for each energy bin - if self.lam_unit in ['Angs', 'Angstrom', 'angs', 'angstrom']: - ener = self.lam * u.angstrom - elif self.lam_unit in ['kev', 'keV']: - ener = self.lam * u.keV - else: - ener = self.lam * u.Unit(self.lam_unit) - int_conv = arf(ener.to(u.keV, equivalencies=u.spectral())) + int_conv = arf(self.lam.to(u.keV, equivalencies=u.spectral()).value) # cm^2 ct/phot r_asec = radius * pix_scale From af096ae07ff0f10a08e34c6c63d483c37fdb3f92 Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Mon, 28 Sep 2026 16:51:50 -0700 Subject: [PATCH 10/12] Failing test for lmin, lmax UnitConversionError --- tests/test_galhalo.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/test_galhalo.py b/tests/test_galhalo.py index 185067b..9838555 100644 --- a/tests/test_galhalo.py +++ b/tests/test_galhalo.py @@ -243,3 +243,10 @@ def test_fake_variable_image(screen_halo, arf): assert np.all(image >= 0.0) assert np.all(image == np.floor(image)) assert image.sum() > 0 + +def test_fake_variable_image_lmin_lmax(screen_halo, arf): + time = np.linspace(0, 100, 10) + lc = np.ones_like(time) + restricted = screen_halo.fake_variable_image(time, lc=lc, arf=arf, exposure=1e4, num_pix=[8, 16], lmin=0.5, lmax=5.0) + unrestricted = screen_halo.fake_variable_image(time, lc=lc, arf=arf, exposure=1e4, num_pix=[8, 16]) + assert 0 < restricted.sum() < unrestricted.sum() From 0d56d27d98ffa1a512fd3bfa16bafe78ced16f30 Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Mon, 28 Sep 2026 16:58:42 -0700 Subject: [PATCH 11/12] Fix lmin/lmax unit comparison in fake_variable_image lmin and lmax are documented as bare floats in halo.lam units but were compared directly against self.lam, an astropy Quantity, raising UnitConversionError whenever either was passed. Give them lam's unit before comparing, as in Halo.fake_image. --- src/xdust/halos/galhalo.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/xdust/halos/galhalo.py b/src/xdust/halos/galhalo.py index 973813c..88a0284 100644 --- a/src/xdust/halos/galhalo.py +++ b/src/xdust/halos/galhalo.py @@ -238,15 +238,16 @@ def fake_variable_image(self, time, lc, arf, """ # intensity cube (NE x NTH), phot/cm^2/s/arcsec^2 var_profile = self.variable_profile(time, lc, dist=dist, tnow=tnow) + # Decide which energy indexes to use if lmin is None: imin = 0 else: - imin = min(np.arange(len(self.lam))[self.lam >= lmin]) + imin = min(np.arange(len(self.lam))[self.lam >= lmin * self.lam.unit]) if lmax is None: iend = len(self.lam) else: - iend = max(np.arange(len(self.lam))[self.lam <= lmax]) + iend = max(np.arange(len(self.lam))[self.lam <= lmax * self.lam.unit]) # set up image grid xlen, ylen = num_pix From 6db5dc062209037abc84a13700f07d5ad8b7599c Mon Sep 17 00:00:00 2001 From: mandarmulherkar Date: Mon, 28 Sep 2026 17:03:49 -0700 Subject: [PATCH 12/12] Strengthen variable_profile and fake_variable_image tests Assert the linearity of variable_profile in lc, which the previous version computed but never checked. pytest.approx does not handle astropy Quantities, so compare .value. --- tests/test_galhalo.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_galhalo.py b/tests/test_galhalo.py index 9838555..f8a30e6 100644 --- a/tests/test_galhalo.py +++ b/tests/test_galhalo.py @@ -233,12 +233,12 @@ def test_variable_profile(screen_halo): var_profile2 = screen_halo.variable_profile(time, 2 * lc, dist=dist) assert var_profile2.shape == (len(EVALS), len(THVALS)) assert np.sum(var_profile2) > 0 - # 2 * lc and assert the result is ~2× + assert np.sum(var_profile2).value == pytest.approx(2 * np.sum(var_profile1).value) def test_fake_variable_image(screen_halo, arf): time = np.linspace(0, 100, 10) lc = np.ones_like(time) - image = screen_halo.fake_variable_image(time, lc=lc, arf=arf, exposure=1e4, num_pix=[8, 16]) + image = screen_halo.fake_variable_image(time, lc=lc, arf=arf, exposure=1e4, pix_scale=1.0, num_pix=[8, 16]) assert image.shape == (16, 8) assert np.all(image >= 0.0) assert np.all(image == np.floor(image))