diff --git a/hyperion/model/image.py b/hyperion/model/image.py index 232ea345..e9ffd81a 100644 --- a/hyperion/model/image.py +++ b/hyperion/model/image.py @@ -132,6 +132,17 @@ def unit(self, value): else: raise ValueError("unit should be a string") + @property + def units(self): + """ + The units of the image values (alias for the ``.unit`` property). + """ + return self.unit + + @units.setter + def units(self, value): + self.unit = value + @property def wav(self): """ diff --git a/hyperion/model/sed.py b/hyperion/model/sed.py index 12d2d3e4..e742b127 100644 --- a/hyperion/model/sed.py +++ b/hyperion/model/sed.py @@ -123,6 +123,17 @@ def unit(self, value): else: raise ValueError("unit should be a string") + @property + def units(self): + """ + The units of the SED values (alias for the ``.unit`` property). + """ + return self.unit + + @units.setter + def units(self, value): + self.unit = value + @property def wav(self): """ diff --git a/hyperion/model/tests/test_image.py b/hyperion/model/tests/test_image.py index 81ee17ec..9dc08311 100644 --- a/hyperion/model/tests/test_image.py +++ b/hyperion/model/tests/test_image.py @@ -728,3 +728,189 @@ def test_get_image_stokes(self, stokes): with pytest.raises(ValueError) as exc: self.m2.get_image(stokes=stokes) assert exc.value.args[0] == "Only the Stokes I value was stored for this image" + + +def _inside_observer_direction(theta, phi): + t, p = np.radians(theta), np.radians(phi) + return np.array([np.sin(t) * np.cos(p), np.sin(t) * np.sin(p), np.cos(t)]) + + +def _inside_observer_expected_lonlat(d, theta_v, phi_v): + # Map (longitude, latitude) of a photon arriving from direction d, for an + # observer looking towards (theta_v, phi_v). This is d expressed in the + # local spherical basis (r_hat, phi_hat, -theta_hat) at the viewing + # direction, converted to polar coordinates - i.e. an independent + # reimplementation of the transform in images_peeled.f90. + t, p = np.radians(theta_v), np.radians(phi_v) + st, ct, cp, sp = np.sin(t), np.cos(t), np.cos(p), np.sin(p) + vx = np.dot([st * cp, st * sp, ct], d) # r_hat . d -> map centre / +x + vy = np.dot([-sp, cp, 0.], d) # phi_hat . d -> +longitude + vz = np.dot([-ct * cp, -ct * sp, st], d) # -theta_hat . d -> +latitude + lon = np.degrees(np.arctan2(vy, vx)) + lat = np.degrees(np.arctan2(np.hypot(vx, vy), vz)) - 90. + return lon, lat + + +@pytest.mark.requires_hyperion_binaries +def test_inside_observer_sky_coordinates(tmpdir): + # Regression test for the inside-observer sky-coordinate transform. A point + # source at cartesian position P, seen by an observer at the origin, sends + # photons that arrive from direction d = -P/|P|. On the all-sky map for a + # viewing direction (theta_v, phi_v) the source must appear at the position + # obtained by expressing d in the local spherical frame of the viewing + # direction. A previous implementation was discontinuous through + # theta_v = 90 deg and ignored phi_v; the viewing angles below include + # several cases it placed incorrectly (large phi at theta=90, either side + # of theta=90, and general off-axis directions). + + # Place the source so that its photon arrival direction is +x. + photon_dir = _inside_observer_direction(90., 0.) + P = tuple(-0.5 * photon_dir) + + views = [(90., 0.), (90., 90.), (90., 150.), (60., 0.), (120., 0.), + (89., 0.), (91., 0.), (45., 30.), (135., 200.)] + + m = Model() + m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.]) + s = m.add_point_source() + s.position = P + s.luminosity = 1. + s.temperature = 6000. + i = m.add_peeled_images(sed=False, image=True) + i.set_inside_observer((0., 0., 0.)) + i.set_image_limits(180., -180., -90., 90.) + i.set_image_size(360, 180) + i.set_viewing_angles([v[0] for v in views], [v[1] for v in views]) + i.set_wavelength_range(1, 1., 1000.) + m.set_n_initial_iterations(0) + m.set_n_photons(imaging=20000) + m.set_seed(-1) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + + val = out.get_image(group=0).val[:, :, :, 0] # (n_view, n_y, n_x) + n_y, n_x = val.shape[1], val.shape[2] + xmin, xmax, ymin, ymax = 180., -180., -90., 90. + + for iv, (theta_v, phi_v) in enumerate(views): + img = val[iv] + ys, xs = np.nonzero(img > 0) + assert len(xs) > 0, "source not found for view (%.0f, %.0f)" % (theta_v, phi_v) + w = img[ys, xs] + lon = ((xmin + (xs + 0.5) * (xmax - xmin) / n_x) * w).sum() / w.sum() + lat = ((ymin + (ys + 0.5) * (ymax - ymin) / n_y) * w).sum() / w.sum() + + # (1) Convention-independent invariant: the great-circle distance on the + # map from the centre (0, 0) to the source must equal the true angle + # between the photon direction and the viewing direction (any rigid + # rotation preserves angles). This catches the tearing and the ignored + # phi of the old code. + map_sep = np.degrees(np.arccos(np.clip( + np.cos(np.radians(lat)) * np.cos(np.radians(lon)), -1., 1.))) + true_sep = np.degrees(np.arccos(np.clip( + np.dot(photon_dir, _inside_observer_direction(theta_v, phi_v)), -1., 1.))) + assert abs(map_sep - true_sep) < 2., \ + "view (%.0f, %.0f): map separation %.1f deg != true %.1f deg" % ( + theta_v, phi_v, map_sep, true_sep) + + # (2) Full position matches the independent reference (catches roll and + # reflection errors that preserve the angular separation). + exp_lon, exp_lat = _inside_observer_expected_lonlat(photon_dir, theta_v, phi_v) + dlon = (lon - exp_lon + 180.) % 360. - 180. + assert abs(dlon) < 2. and abs(lat - exp_lat) < 2., \ + "view (%.0f, %.0f): (%.1f, %.1f) != expected (%.1f, %.1f)" % ( + theta_v, phi_v, lon, lat, exp_lon, exp_lat) + + +def test_image_and_sed_units_constructor(): + # Regression test: the ``units`` constructor argument of Image and SED must + # be routed through the validated ``unit`` property, so that it is stored + # (readable via ``.unit``) and invalid values are rejected. Previously it + # was assigned to an unvalidated attribute and left ``.unit`` unset. + from ..sed import SED + for cls in (Image, SED): + obj = cls(nu=np.array([1.e10, 1.e11]), units='Jy') + assert obj.unit == 'Jy' + assert obj.units == 'Jy' + with pytest.raises(ValueError): + cls(nu=np.array([1.e10, 1.e11]), units=42) + + +@pytest.mark.requires_hyperion_binaries +def test_inside_observer_flux_dilution(tmpdir): + # Regression test for the inside-observer 1/d^2 flux dilution. Two point + # sources of equal luminosity at distances d1 and d2 from the observer must + # have observed fluxes in the ratio (d2/d1)**2. A previous implementation + # divided by (d - d_min)**2 instead of d**2, giving the wrong ratio whenever + # a nonzero depth minimum was set. + d1, d2, d_min = 2., 5., 1. + m = Model() + m.set_cartesian_grid([-8., 8.], [-8., 8.], [-8., 8.]) + for pos in [(d1, 0., 0.), (0., d2, 0.)]: + s = m.add_point_source() + s.position = pos + s.luminosity = 1. + s.temperature = 6000. + i = m.add_peeled_images(sed=False, image=True) + i.set_inside_observer((0., 0., 0.)) + i.set_image_limits(180., -180., -90., 90.) + i.set_image_size(360, 180) + i.set_viewing_angles([90.], [0.]) + i.set_wavelength_range(1, 1., 1000.) + i.set_depth(d_min, 20.) # nonzero d_min exposes the bug + m.set_n_initial_iterations(0) + m.set_n_photons(imaging=200000) + m.set_seed(-1) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + + val = out.get_image(group=0).val[0, :, :, 0] # single view, single wavelength + brightest = np.sort(val[val > 0])[::-1] + assert len(brightest) >= 2, "expected two point sources in the image" + # the closer source (d1) is the brighter one + ratio = brightest[0] / brightest[1] + np.testing.assert_allclose(ratio, (d2 / d1) ** 2, rtol=0.08) + + +@pytest.mark.requires_hyperion_binaries +def test_inside_observer_peeloff_optical_depth(tmpdir): + # Regression test for integrating the peeloff optical depth all the way to + # the (inside) observer. A point source at distance d sits in uniform, + # purely-absorbing dust of flat opacity chi and density rho, so its direct + # light is attenuated by exp(-chi*rho*d) over the full path to the observer. + # A previous implementation integrated the optical depth only to d - d_min, + # under-attenuating the source whenever a nonzero depth minimum was set + # (comparing the dust/no-dust flux ratio at fixed distance cancels the + # 1/d^2 dilution, isolating the optical-depth path). + from ...dust import IsotropicDust + + d, chi, rho, d_min = 1.e16, 1., 1.e-16, 0.5e16 + tau_full = chi * rho * d # = 1.0 + + def peak_flux(with_dust): + m = Model() + m.set_cartesian_grid([-2.e16, 2.e16], [-2.e16, 2.e16], [-2.e16, 2.e16]) + if with_dust: + dust = IsotropicDust([3.e9, 3.e16], [0., 0.], [chi, chi]) + dust.set_lte_emissivities(n_temp=10, temp_min=0.1, temp_max=1.e4) + m.add_density_grid(np.ones((1, 1, 1)) * rho, dust) + s = m.add_point_source() + s.position = (d, 0., 0.) + s.luminosity = 1. + s.temperature = 6000. + i = m.add_peeled_images(sed=False, image=True) + i.set_inside_observer((0., 0., 0.)) + i.set_image_limits(180., -180., -90., 90.) + i.set_image_size(360, 180) + i.set_viewing_angles([90.], [0.]) + i.set_wavelength_range(1, 1., 1000.) + i.set_depth(d_min, 3.e16) # nonzero d_min previously truncated the tau path + m.set_n_initial_iterations(0) + m.set_n_photons(imaging=100000) + m.set_seed(-1) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + return out.get_image(group=0).val[0, :, :, 0].max() + + ratio = peak_flux(with_dust=True) / peak_flux(with_dust=False) + np.testing.assert_allclose(ratio, np.exp(-tau_full), rtol=0.05) diff --git a/hyperion/model/tests/test_sed.py b/hyperion/model/tests/test_sed.py index 6e86526e..3f94498d 100644 --- a/hyperion/model/tests/test_sed.py +++ b/hyperion/model/tests/test_sed.py @@ -470,3 +470,34 @@ def test_sed_negative_group(tmpdir): m.write(tmpdir.join(random_id()).strpath) mo = m.run() np.testing.assert_allclose(mo.get_sed(group=-1).val, mo.get_sed(group=0).val) + + +@pytest.mark.requires_hyperion_binaries +def test_sed_uncertainty_sum_of_squares(tmpdir): + # Regression test for the Monte-Carlo uncertainty estimator. For a single + # isotropic point source with no dust, every imaging photon contributes the + # same weight w to the single SED bin, so the sum-of-squares estimator gives + # sigma = sqrt(sum w**2) = sqrt(N) w against a total flux F = N w, i.e. a + # relative uncertainty of exactly 1/sqrt(N). A previous (sign-flipped) + # sample-variance formula over-estimated this by a factor sqrt(2). + n_photons = 10000 + m = Model() + m.set_cartesian_grid([-1., 1.], [-1., 1.], [-1., 1.]) + s = m.add_point_source() + s.luminosity = 1. + s.temperature = 6000. + i = m.add_peeled_images(sed=True, image=False) + i.set_viewing_angles([45.], [45.]) + i.set_aperture_radii(1, 1.e10, 1.e10) + i.set_wavelength_range(1, 0.01, 5000.) + i.set_uncertainties(True) + m.set_n_initial_iterations(0) + m.set_n_photons(imaging=n_photons) + m.set_seed(-1) + m.write(tmpdir.join(random_id()).strpath) + out = m.run(tmpdir.join(random_id()).strpath) + + sed = out.get_sed(uncertainties=True) + flux = np.nansum(sed.val) + sigma = np.sqrt(np.nansum(np.array(sed.unc) ** 2)) + np.testing.assert_allclose(sigma / flux, 1. / np.sqrt(n_photons), rtol=0.03) diff --git a/src/images/image_type.f90 b/src/images/image_type.f90 index 558fb82b..57e51243 100644 --- a/src/images/image_type.f90 +++ b/src/images/image_type.f90 @@ -665,13 +665,12 @@ subroutine image_write(img,group) cube5d = img%sed - if(img%uncertainties) then - where(img%sedn > 1) - cube5de = sqrt((img%sed2 + (img%sed)**2 / img%sedn) / (img%sedn - 1)) * sqrt(img%sedn) - elsewhere - cube5de = 0._dp - end where - end if + ! The uncertainty on the total flux in each bin is estimated with the + ! standard Monte-Carlo sum-of-squares estimator sqrt(sum(x_i**2)), + ! which accounts for both the scatter in the photon weights and the + ! Poisson counting uncertainty. + + if(img%uncertainties) cube5de = sqrt(img%sed2) if(.not.img%use_exact_nu) then cube5d = cube5d / dnunorm @@ -727,13 +726,9 @@ subroutine image_write(img,group) cube6d = img%img - if(img%uncertainties) then - where(img%imgn > 1) - cube6de = sqrt((img%img2 + (img%img)**2 / img%imgn) / (img%imgn - 1)) * sqrt(img%imgn) - elsewhere - cube6de = 0._dp - end where - end if + ! See the SED branch above for the derivation of the uncertainties + + if(img%uncertainties) cube6de = sqrt(img%img2) if(.not.img%use_exact_nu) then cube6d = cube6d / dnunorm diff --git a/src/images/images_peeled.f90 b/src/images/images_peeled.f90 index 45eed6cc..2e535bf4 100644 --- a/src/images/images_peeled.f90 +++ b/src/images/images_peeled.f90 @@ -100,8 +100,8 @@ subroutine peeloff_photon(p_orig, polychromatic) type(photon) :: p real(dp) :: tau integer :: ip,ig,iv,id - type(angle3d_dp) :: a_req, a_sky, a_diff - type(vector3d_dp) :: v_req + type(angle3d_dp) :: a_req, a_view + type(vector3d_dp) :: v_req, v_a, v_sky logical,intent(in) :: polychromatic real(dp) :: x_image, y_image real(dp) :: tmax @@ -151,29 +151,53 @@ subroutine peeloff_photon(p_orig, polychromatic) ! the xmax wall. call place_in_cell(p) + ! The depth range given by d_min and d_max is only used to select + ! which events to include in the image - the optical depth should be + ! integrated all the way to the observer (for inside observers) or to + ! the edge of the grid (for external observers). if(inside_observer(ig)) then dr = p%r-r_peeloff(ig) d = sqrt(dr.dot.dr) - tmax = d - d_min(ig) + tmax = d else d = -(v_req.dot.p%r) - tmax = - d_min(ig) + d + tmax = huge(1._dp) end if if(d < d_min(ig) .or. d > d_max(ig)) cycle if(inside_observer(ig)) then - if(abs(viewing_angles(ip)%cost) .gt. 1.e-10) then - call rotate_angle3d(viewing_angles(ip), angle3d_deg(90._dp, 0._dp), a_diff) - call difference_angle3d(a_diff, -p%a, a_sky) - else - a_sky = p%a - end if - - ! Convert to angles in degrees - x_image = atan2(a_sky%sinp, a_sky%cosp) * rad2deg - y_image = atan2(a_sky%sint, a_sky%cost) * rad2deg - 90._dp + ! Sky position of the photon for an observer looking towards the + ! direction a_view = (theta_v, phi_v). The map is centred on the + ! viewing direction, so we express the photon direction p%a in the + ! local orthonormal frame at a_view given by the spherical basis + ! vectors (r_hat, phi_hat, -theta_hat), writing t, p for theta_v, + ! phi_v: + ! + ! r_hat = ( sin t cos p, sin t sin p, cos t ) (viewing dir) + ! phi_hat = ( -sin p, cos p, 0 ) (local east) + ! -theta_hat = (-cos t cos p, -cos t sin p, sin t ) (local north) + ! + ! The three components of v_sky below are the projections of p%a onto + ! these axes, i.e. v_sky = R p%a where R is the rigid rotation whose + ! rows are the vectors above and which therefore sends a_view to the + ! map centre +x = (theta, phi) = (90, 0). R is orthonormal, so this is + ! exactly continuous in (theta_v, phi_v) with no special cases at the + ! poles (theta_v = 0, 90, 180) - unlike a frame built from the + ! pole-referenced rotate_angle3d / difference_angle3d routines. + a_view = viewing_angles(ip) + call angle3d_to_vector3d(p%a, v_a) + v_sky%x = (v_a%x * a_view%cosp + v_a%y * a_view%sinp) * a_view%sint + v_a%z * a_view%cost + v_sky%y = - v_a%x * a_view%sinp + v_a%y * a_view%cosp + v_sky%z = - (v_a%x * a_view%cosp + v_a%y * a_view%sinp) * a_view%cost + v_a%z * a_view%sint + + ! Convert to map coordinates: longitude is the azimuth of v_sky in the + ! x-y plane, and latitude is (colatitude measured from +z) - 90 deg. + ! Both use atan2 on the Cartesian components, so they stay well-defined + ! when v_sky is near the poles (where sin(colatitude) -> 0). + x_image = atan2(v_sky%y, v_sky%x) * rad2deg + y_image = atan2(sqrt(v_sky%x**2 + v_sky%y**2), v_sky%z) * rad2deg - 90._dp ! Make sure the photon falls inside the image (wrap angles around) x_image = peeled_image(ig)%x_max + modulo(x_image - peeled_image(ig)%x_max, 360._dp) @@ -205,12 +229,11 @@ subroutine peeloff_photon(p_orig, polychromatic) end if end if - ! For inside observer, don't want optical depth to escape grid, just to go to observer! - ! Need to include 1/d^2! - if(.not.killed) then - if(inside_observer(ig)) p%s = p%s / (4._dp * pi * tmax**2._dp) + ! For inside observers, include the 1/d^2 flux dilution + + if(inside_observer(ig)) p%s = p%s / (4._dp * pi * d**2._dp) if(polychromatic) then