Skip to content

Do better with waveranges #936

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
hugobuddel wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Do better with waveranges #936

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7171eae
Round waverange to 8 decimal places to prevent floating point errors
hugobuddel Jun 1, 2026
602e9cd
Include endpoint of spectrum in prepare_source
hugobuddel Jun 1, 2026
f536277
Fix test
hugobuddel Jun 1, 2026
6b8bcbf
Properly recalculate fov_waveset
hugobuddel Jun 1, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions scopesim/optics/fov.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1165,8 +1165,10 @@ def extract_range_from_spectrum(spectrum, waverange):
assert isinstance(spectrum, SourceSpectrum), (
f"spectrum must be of type synphot.SourceSpectrum: {type(spectrum)}")

wave_min, wave_max = quantify(waverange, u.um).to(u.AA).value
spec_waveset = spectrum.waveset.to(u.AA).value
# Round to 8 decimal places, because the conversion to u.AA will cause
# irrelevant floating point differences.
wave_min, wave_max = quantify(waverange, u.um).to(u.AA).value.round(8)
spec_waveset = spectrum.waveset.to(u.AA).value.round(8)
mask = (spec_waveset > wave_min) * (spec_waveset < wave_max)

# FIXME: Why did I comment this out in 2023? Seems useful to have...
Expand Down
15 changes: 13 additions & 2 deletions scopesim/optics/optical_train.py
View file
Open in desktop
Comment thread
teutoburg marked this conversation as resolved.
Original file line number Diff line number Diff line change
Expand Up @@ -340,7 +340,12 @@ def prepare_source(self, source):
wave_max = max(fov.meta["wave_max"] for fov in self.fov_manager.fovs)
wave_unit = u.Unit(from_currsys("!SIM.spectral.wave_unit", self.cmds))
dwave = from_currsys("!SIM.spectral.spectral_bin_width", self.cmds) # Not a quantity
fov_waveset = np.arange(wave_min.value, wave_max.value, dwave) * wave_unit
# Include the last point by adding dwave, because arange excludes the last point.
fov_waveset = np.arange(wave_min.value, wave_max.value + 0.9 * dwave, dwave)
# Ensure the start and end values are exact.
fov_waveset[0] = wave_min.value
fov_waveset[-1] = wave_max.value
fov_waveset *= wave_unit
fov_waveset = fov_waveset.to(u.um)

field.spectra[ispec] = SourceSpectrum(Empirical1D,
Expand Down Expand Up @@ -382,11 +387,17 @@ def prepare_source(self, source):
cube.header["CUNIT2"] = "deg"

# Put on fov wavegrid
# TODO: This is the same code as earlier in the function; refactor.
wave_min = min(fov.meta["wave_min"] for fov in self.fov_manager.fovs)
wave_max = max(fov.meta["wave_max"] for fov in self.fov_manager.fovs)
wave_unit = u.Unit(from_currsys("!SIM.spectral.wave_unit", self.cmds))
dwave = from_currsys("!SIM.spectral.spectral_bin_width", self.cmds) # Not a quantity
fov_waveset = np.arange(wave_min.value, wave_max.value, dwave) * wave_unit
# Include the last point by adding dwave, because arange excludes the last point.
fov_waveset = np.arange(wave_min.value, wave_max.value + 0.9 * dwave, dwave)
# Ensure the start and end values are exact.
fov_waveset[0] = wave_min.value
fov_waveset[-1] = wave_max.value
fov_waveset *= wave_unit
fov_waveset = fov_waveset.to(u.um)

if (wave.to(u.um).min() > fov_waveset.max() or
Expand Down
5 changes: 4 additions & 1 deletion scopesim/tests/tests_optics/test_fov_utls.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,10 @@ def test_extracts_the_wave_range_needed(self):
new_spec = extract_range_from_spectrum(spec, waverange)

assert len(new_spec.waverange) == 2
assert new_spec.waverange[0] == 1.98 * u.um
# Cannot test for equality, because this does not hold in astropy.units:
# 1.98 * u.um == 19800.0 * u.AA
# assert new_spec.waverange[0] == 1.98 * u.um
assert abs(new_spec.waverange[0] == 1.98 * u.um) < 1e-10 * u.AA
Copy link
Copy Markdown
Contributor

@teutoburg teutoburg Jun 1, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't quite understand this line. You're checking for the absolute values of the equality check? Wouldn't that always be a bool?

Copy link
Copy Markdown
Collaborator Author

@hugobuddel hugobuddel Jun 3, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oops, that should be a - sign in there, not an ==

Copy link
Copy Markdown
Contributor

@teutoburg teutoburg Jun 3, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then it looks like something where we should really use numpy.testing instead...

assert new_spec(1.98 * u.um).value == approx(12.8)

@pytest.mark.parametrize(("endpoint", "msg"),
Expand Down
Loading