Skip to content

Fixes to FOV and implane-utils for correct flat_lamp scaling #1007

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

Closed
janusbrink wants to merge 1 commit into from
+25 −9
Closed

Fixes to FOV and implane-utils for correct flat_lamp scaling #1007

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
28 changes: 21 additions & 7 deletions scopesim/optics/fov.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -605,15 +605,17 @@ def _make_scaled_hdu(self, field: ImageSourceField) -> fits.ImageHDU:
"""Deal with BUNIT and pixel area."""
# Note: Do not scale source data - make a copy first.
field_hdu = field.field.copy() # .field is the HDU (yeah...)
field_hdu.data /= self.pixel_area.value
field_hdu.data /= field.pixel_area.value

# TODO: Check if this scaling is actually correct. How does this
# work with the add_imagehdu_to_imagehdu below? Isn't that
# supposed to conserve flux? Test carefully!!
if field.is_bunit_spatially_differential:
logger.debug("differential bunit...")
# Field is in (PHOTLAM) arcsec-2, need to scale by pixarea
field_hdu.data *= field.pixel_area.value
field_hdu.data *= self.pixel_area.value
else:
logger.debug("binned bunit...")
# Pixel area doesn't cancel out, need to convert
new_bunit = field.bunit / u.arcsec**2
field_hdu.header["BUNIT"] = new_bunit.to_string("fits")
Expand Down Expand Up @@ -751,6 +753,7 @@ def _make_cubefields(self):
image = np.sum(field.data, axis=0) * PHOTLAM/u.arcsec**2
image = (image * self.area *
spectral_bin_width).to(u.ph/u.s/u.arcsec**2)
logger.debug("2D FOV make_cubefields: image.mean() = %f", image.mean().value)
# FIXME: This might create a 2D ImageHDU with a 3D header, not ideal...
yield fits.ImageHDU(data=image, header=field.header)

Expand All @@ -763,6 +766,7 @@ def _make_imagefields(self, fov_waveset, bin_widths):
* yield image to be added to canvas image
"""
for field in self._get_image_fields():
logger.debug("2D FOV make_imagefields: field.data.sum() = %f", field.data.sum())
field_hdu = self._make_scaled_hdu(field)

# Reevaluate spectrum onto FOV waveset
Expand All @@ -773,6 +777,7 @@ def _make_imagefields(self, fov_waveset, bin_widths):

# Rescale 2D flux weight map by integrated flux from spectrum
field_hdu.data *= flux # ph s-1 arcsec-2
logger.debug("2D FOV make_imagefields: field_hdu.data.mean() = %f ph s-1 arcsec-2", field_hdu.data.mean())
yield field_hdu

def _make_tablefields(self, fov_waveset, bin_widths, use_photlam=False):
Expand Down Expand Up @@ -851,11 +856,13 @@ def make_hdu(self):
canvas_image_hdu = imp_utils.add_imagehdu_to_imagehdu(
tmp_hdu,
canvas_image_hdu,
conserve_flux=True,
conserve_flux=False,
spline_order=self.spline_order)
logger.debug("2D FOV make_hdu: canvas_image_hdu.data.mean() = %f ph s-1 arcsec-2", canvas_image_hdu.data.mean())

# TODO: Move this further along at some point...
canvas_image_hdu.data *= self.pixel_area.value # eliminate arcsec-2
logger.debug("After scaling by pixel area: canvas_image_hdu.data.mean() = %f ph s-1 pix-1", canvas_image_hdu.data.mean())

for flux, weight, x, y in self._make_tablefields(
fov_waveset, bin_widths):
Expand Down Expand Up @@ -932,7 +939,8 @@ def _make_cubefields(self, fov_waveset):
field_data = field_interp(fov_waveset.value)

# Pixel scale conversion
field_data *= field.pixel_area / self.pixel_area
# field_data *= field.pixel_area / self.pixel_area ## not needed in flux density cubes??
logger.debug("3D FOV make_cubefields: field_data.mean() = %f PHOTLAM arcsec-2", field_data.mean())
field_hdu = fits.ImageHDU(data=field_data, header=field.header)
yield field_hdu

Expand All @@ -957,12 +965,14 @@ def _make_imagefields(self, fov_waveset, spline_order):
canvas_image_hdu = imp_utils.add_imagehdu_to_imagehdu(
field_hdu,
canvas_image_hdu,
spline_order=spline_order)
spline_order=spline_order,
conserve_flux=False,
)

spec = field.spectrum(fov_waveset)

# 2D * 1D -> 3D
field_cube = canvas_image_hdu.data[None, :, :] * spec[:, None, None]
logger.debug("3D FOV make_imagefields: field_cube.mean() = %f", field_cube.mean().value)
yield field_cube.value

def _make_tablefields(self, fov_waveset):
Expand Down Expand Up @@ -1113,7 +1123,9 @@ def make_hdu(self):
canvas_cube_hdu = imp_utils.add_imagehdu_to_imagehdu(
field_hdu,
canvas_cube_hdu,
spline_order=self.spline_order)
spline_order=self.spline_order,
conserve_flux=False,
)

canvas_cube_hdu.data = sum(self._make_imagefields(
fov_waveset, self.spline_order),
Expand All @@ -1132,6 +1144,8 @@ def make_hdu(self):
canvas_cube_hdu.data *= 1e4 # ph/s/AA/arcsec2 --> ph/s/um/arcsec2
canvas_cube_hdu.header["BUNIT"] = "ph s-1 um-1 arcsec-2"

logger.debug("3D FOV make_hdu: canvas_cube_hdu.data.mean() = %f", canvas_cube_hdu.data.mean())

return canvas_cube_hdu # [ph s-1 um-1 (arcsec-2)]

def plot_data(self):
Expand Down
6 changes: 4 additions & 2 deletions scopesim/optics/image_plane_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -590,7 +590,9 @@ def rescale_imagehdu(imagehdu: fits.ImageHDU, pixel_scale: float | u.Quantity,
new_im = np.nan_to_num(new_im, copy=False)
sum_new = np.sum(new_im)
if sum_new != 0:
new_im *= sum_orig / sum_new
flux_factor = sum_orig / sum_new
logger.debug("flux factor = %f", flux_factor)
new_im *= flux_factor
elif sum_orig != 0:
logger.warning(
"rescale_imagehdu: all input flux (%g) was lost in resampling; "
Expand Down Expand Up @@ -835,7 +837,7 @@ def add_imagehdu_to_imagehdu(image_hdu: fits.ImageHDU,
spline_order=spline_order,
conserve_flux=conserve_flux)
# TODO: Perhaps add separately formatted WCS logger?
# logger.debug("fromrescale %s", WCS(new_hdu.header, key=canvas_wcs.wcs.alt))
logger.debug("fromrescale %s", WCS(new_hdu.header, key=canvas_wcs.wcs.alt))
new_hdu = reorient_imagehdu(new_hdu,
wcs_suffix=canvas_wcs.wcs.alt,
spline_order=spline_order,
Expand Down
Loading