From 7eff774e887730eb73dfff262155e0a41c0a918a Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Sat, 9 Nov 2024 08:33:08 -0500 Subject: [PATCH 01/31] New units.typing esp. to simplify spectral flux density and radiance annotations. --- sbpy/units/typing.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) create mode 100644 sbpy/units/typing.py diff --git a/sbpy/units/typing.py b/sbpy/units/typing.py new file mode 100644 index 00000000..bd65c71f --- /dev/null +++ b/sbpy/units/typing.py @@ -0,0 +1,14 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +"""sbpy unit and quantity typing.""" + +from typing import Union +import astropy.units as u + +UnitLike = Union[str, u.Unit] +SpectralQuantity = Union[u.Quantity["length"], u.Quantity["frequency"]] +SpectralFluxDensityQuantity = Union[ + u.Quantity["spectral_flux_density"], u.Quantity["spectral_flux_density_wav"] +] +SpectralRadianceQuantity = Union[ + u.Quantity["surface_brightness"], u.Quantity["surface_brightness_wav"] +] From 7756beb9c87fe2c3076dbe861f4d6750b3f96d35 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 08:40:54 -0500 Subject: [PATCH 02/31] New surfaces module. --- docs/index.rst | 7 + docs/sbpy/surfaces.rst | 37 +++++ docs/static/scattering-vectors.svg | 127 +++++++++++++++ sbpy/surfaces/__init__.py | 0 sbpy/surfaces/lambertian.py | 27 ++++ sbpy/surfaces/scattered.py | 121 +++++++++++++++ sbpy/surfaces/surface.py | 207 +++++++++++++++++++++++++ sbpy/surfaces/tests/__init__.py | 0 sbpy/surfaces/tests/test_lambertian.py | 43 +++++ sbpy/surfaces/tests/test_scattered.py | 31 ++++ sbpy/surfaces/tests/test_surface.py | 81 ++++++++++ sbpy/surfaces/thermal.py | 104 +++++++++++++ 12 files changed, 785 insertions(+) create mode 100644 docs/sbpy/surfaces.rst create mode 100644 docs/static/scattering-vectors.svg create mode 100644 sbpy/surfaces/__init__.py create mode 100644 sbpy/surfaces/lambertian.py create mode 100644 sbpy/surfaces/scattered.py create mode 100644 sbpy/surfaces/surface.py create mode 100644 sbpy/surfaces/tests/__init__.py create mode 100644 sbpy/surfaces/tests/test_lambertian.py create mode 100644 sbpy/surfaces/tests/test_scattered.py create mode 100644 sbpy/surfaces/tests/test_surface.py create mode 100644 sbpy/surfaces/thermal.py diff --git a/docs/index.rst b/docs/index.rst index 33c62534..129213d4 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -95,6 +95,13 @@ Activity sbpy/activity/index.rst +Surfaces and Shapes +------------------- + +.. toctree:: :maxdepth: 2 + + sbpy/surfaces + Miscellaneous ------------- diff --git a/docs/sbpy/surfaces.rst b/docs/sbpy/surfaces.rst new file mode 100644 index 00000000..fc615992 --- /dev/null +++ b/docs/sbpy/surfaces.rst @@ -0,0 +1,37 @@ +Surfaces Module (`sbpy.surfaces`) +================================= + +The ``surfaces`` module describes the interaction of electromagnetic radiation with surfaces. Sbpy uses the :math:`(i, e, \phi)` model (angle of incidence, angle of emittance, and phase angle) to describe how absorptance, emittance, and reflectance vary with incoming and outgoing radiation. It has a flexible system that can incorporate any surface scattering model that can be described with these three angles. However, most users will use the built-in surface models. + + +.. figure:: ../static/scattering-vectors.svg + :alt: Diagram of surface scattering and emission vectors + + Sbpy's geometric basis for surface scattering and emission: :math:`n` is the surface normal vector, :math:`r_s` is the radial vector to the light source, and :math:`r_o` is the radial vector to the observer. The angle of incidence (:math:`i`), angle of emittance (:math:`e`), phase angle (:math:`\phi`) are labeled. + +A ``Surface`` as methods to + + +Built-in surface models +----------------------- + +The model `sbpy.surfaces.scattered.LambertianSurfaceScatteredSunlight` is used to observe sunlight scattered from a Lambertian surface (light scattered uniformly in all directions). + +Create an instance of the ``LambertianSurfaceScatteredSunlight`` model, and calculate the absorptance, emittance, and reflectance for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)`. + +.. code:: python + + >>> import astropy.units as u + >>> import matplotlib.pyplot as plt + >>> from sbpy.surfaces.scattered import LambertianSurfaceScatteredSunlight + >>> + >>> surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) + >>> + >>> i, e, phi = [30, 60, 90] * u.deg + >>> surface.absorptance(i) # doctest: +FLOAT_CMP + + >>> surface.emittance(e, phi) # doctest: +FLOAT_CMP + + >>> surface.reflectance(i, e, phi) # doctest: +FLOAT_CMP + + diff --git a/docs/static/scattering-vectors.svg b/docs/static/scattering-vectors.svg new file mode 100644 index 00000000..0e437bea --- /dev/null +++ b/docs/static/scattering-vectors.svg @@ -0,0 +1,127 @@ + + + + + + + + + + + + + + + + + + i + e + φ + n + rs + ro + + diff --git a/sbpy/surfaces/__init__.py b/sbpy/surfaces/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/sbpy/surfaces/lambertian.py b/sbpy/surfaces/lambertian.py new file mode 100644 index 00000000..6b1d4e1d --- /dev/null +++ b/sbpy/surfaces/lambertian.py @@ -0,0 +1,27 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import numpy as np + +import astropy.units as u + +from .surface import Surface + + +class LambertianSurface(Surface): + """Abstract base class for Lambertian surfaces.""" + + @u.quantity_input + def absorptance(self, i: u.physical.angle) -> u.Quantity: + cos_i = np.maximum(self._cos(i), np.zeros_like(i.value)) + return (1 - self.phys["albedo"]) * cos_i + + @u.quantity_input + def emittance(self, e: u.physical.angle, phi: u.physical.angle) -> u.Quantity: + return np.maximum(self._cos(e), np.zeros_like(e.value)) + + @u.quantity_input + def reflectance( + self, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle + ) -> u.Quantity: + cos_i = np.maximum(self._cos(i), np.zeros_like(i.value)) + return self.phys["albedo"] * cos_i * self.emittance(e, phi) diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py new file mode 100644 index 00000000..be335249 --- /dev/null +++ b/sbpy/surfaces/scattered.py @@ -0,0 +1,121 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import numpy as np + +import astropy.units as u + +from .surface import Surface +from .lambertian import LambertianSurface +from ..calib import Sun +from ..units.typing import SpectralQuantity, SpectralRadianceQuantity, UnitLike + + +class ScatteredSunlight(Surface): + """Abstract base class to observe light scattered by a surface illuminated by the Sun.""" + + @u.quantity_input + def radiance( + self, + wave_freq: SpectralQuantity, + i: u.physical.angle, + e: u.physical.angle, + phi: u.physical.angle, + rh: u.physical.length, + unit: UnitLike = "W/(m2 sr um)", + ) -> u.Quantity: + """Observed radiance from a surface. + + + Parameters + ---------- + wave_freq : `astropy.units.Quantity` + Wavelength or frequency at which to evaluate the Sun. Arrays are + evaluated with `sbpy.calib.core.Sun.observe()`. + + i : `~astropy.units.Quantity` + Angle from normal of incident light. + + e : `~astropy.units.Quantity` + Angle from normal of emitted light. + + phi : `~astropy.units.Quantity` + Source-target-observer (phase) angle. + + rh : `~astropy.units.Quantity` + Heliocentric distance. + + unit : `~astropy.units.Unit`, optional + Unit of the return value. + + + Returns + ------- + radiance : `~astropy.units.Quantity` + Observed radiance. + + """ + + sun = Sun.from_default() + if wave_freq.size == 1: + F_i = sun(wave_freq) + else: + F_i = sun.observe(wave_freq) + + F_i /= rh.to_value("au") ** 2 + return (F_i * self.reflectance(i, e, phi) / u.sr).to(unit) + + @u.quantity_input + def radiance_from_vectors( + self, + wave_freq: SpectralQuantity, + n: np.ndarray[3], + rs: u.physical.length, + ro: u.physical.length, + unit: UnitLike = "W/(m2 sr um)", + ) -> u.Quantity: + """Observed radiance from a surface. + + + Parameters + ---------- + wave_freq : `astropy.units.Quantity` + Wavelength or frequency at which to evaluate the Sun. Arrays are + evaluated with `sbpy.calib.core.Sun.observe()`. + + n : `numpy.ndarray` + Surface normal vector. + + rs : `~astropy.units.Quantity` + Radial vector from the surface to the light source. + + ro : `numpy.ndarray` + Radial vector from the surface to the observer. + + rh : `~astropy.units.Quantity` + Heliocentric distance. + + unit : `~astropy.units.Unit`, optional + Unit of the return value. + + + Returns + ------- + radiance : `~astropy.units.Quantity` + Observed radiance. + + """ + rh = np.linalg.norm(rs).to("au") + i, e, phi = self._vectors_to_angles(n, rs, ro) + return self.radiance(wave_freq, i, e, phi, rh, unit=unit) + + __doc__ += Surface.__doc__[Surface.__doc__.index("\n") + 1 :] + + +class LambertianSurfaceScatteredSunlight(LambertianSurface, ScatteredSunlight): + """Sunlight scattered from a Lambertian surface. + + The surface is assumed to be illuminated by the Sun. + + """ + + __doc__ += Surface.__doc__[Surface.__doc__.index("\n") :] diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py new file mode 100644 index 00000000..c6cf6101 --- /dev/null +++ b/sbpy/surfaces/surface.py @@ -0,0 +1,207 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +from typing import Tuple +from abc import ABC, abstractmethod + +import numpy as np +from astropy import units as u + +from ..data.phys import Phys +from ..data.decorators import dataclass_input +from ..units.typing import ( + SpectralFluxDensityQuantity, + SpectralQuantity, + SpectralRadianceQuantity, + UnitLike, +) + + +class Surface(ABC): + """Abstract base class for all small-body surfaces. + + + Parameters + ---------- + phys : Phys + Surface physical parameters, e.g., albedo. + + """ + + @dataclass_input + def __init__(self, phys: Phys): + self.phys = phys + + @staticmethod + def _cos(a: u.physical.angle) -> u.Quantity: + """Use to ensure that cos(90 deg) equals 0.""" + + # handle scalars separately + if a.ndim == 0 and u.isclose(a, 90 * u.deg): + return u.Quantity(0) + + x = np.cos(a) + x[u.isclose(a, 90 * u.deg)] = 0 + + return x + + @abstractmethod + def absorptance(self, i: u.physical.angle) -> u.Quantity: + r"""Absorption of incident light. + + The surface is illuminated by incident flux density, :math:`F_i`, at an + angle of :math:`i`, measured from the surface normal direction. + + + Parameters + ---------- + i : `~astropy.units.Quantity` + Angle from normal of incident light. + + + Returns + ------- + a : `~astropy.units.Quantity` + Surface absorptance. + + """ + + @abstractmethod + def emittance(self, e: u.physical.angle, phi: u.physical.angle) -> u.Quantity: + r"""Emission of light. + + The surface is observed at an angle of :math:`e`, measured from the + surface normal direction, and at a solar phase angle of :math:`phi`. + + + Parameters + ---------- + e : `~astropy.units.Quantity` + Angle from normal of emitted light. + + phi : `~astropy.units.Quantity` + Source-target-observer (phase) angle. + + + Returns + ------- + e : `~astropy.units.Quantity` + Surface emittance. + + """ + + @abstractmethod + def reflectance( + self, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle + ) -> u.Quantity: + r"""Reflectance. + + The surface is illuminated by incident flux density, :math:`F_i`, at an + angle of :math:`i`, and emitted toward an angle of :math:`e`, measured + from the surface normal direction. :math:`\phi` is the + source-target-observer (phase) angle. + + + Parameters + ---------- + i : `~astropy.units.Quantity` + Angle from normal of incident light. + + e : `~astropy.units.Quantity` + Angle from normal of emitted light. + + phi : `~astropy.units.Quantity` + Source-target-observer (phase) angle. + + + Returns + ------- + r : `~astropy.units.Quantity` + Surface reflectance. + + """ + + @abstractmethod + def radiance( + self, + F_i: SpectralFluxDensityQuantity, + i: u.physical.angle, + e: u.physical.angle, + phi: u.physical.angle, + ) -> u.Quantity: + """Observed radiance from a surface. + + + Parameters + ---------- + F_i : `astropy.units.Quantity` + Incident light (spectral flux density). + + i : `~astropy.units.Quantity` + Angle from normal of incident light. + + e : `~astropy.units.Quantity` + Angle from normal of emitted light. + + phi : `~astropy.units.Quantity` + Source-target-observer (phase) angle. + + + Returns + ------- + radiance : `~astropy.units.Quantity` + Observed radiance. + + """ + + @staticmethod + def _vectors_to_angles( + n: np.ndarray[3], + rs: u.physical.length, + ro: np.ndarray[3], + ) -> Tuple[u.physical.angle, u.physical.angle, u.physical.angle]: + n_hat = n / np.linalg.norm(n) + rs_hat = rs / np.linalg.norm(rs) + ro_hat = ro / np.linalg.norm(ro) + + i = u.Quantity(np.arccos(np.dot(n_hat, rs_hat)), "rad") + e = u.Quantity(np.arccos(np.dot(n_hat, ro_hat)), "rad") + phi = u.Quantity(np.arccos(np.dot(rs_hat, ro_hat)), "rad") + + return i, e, phi + + @u.quantity_input + def radiance_from_vectors( + self, + F_i: SpectralFluxDensityQuantity, + n: np.ndarray[3], + rs: u.physical.length, + ro: np.ndarray[3], + ) -> u.Quantity: + """Observed radiance from a surface with geometry defined by vectors. + + Input vectors do not need to be normalized. + + + Parameters + ---------- + F_i : `astropy.units.Quantity` + Incident light (spectral flux density). + + n : `numpy.ndarray` + Surface normal vector. + + rs : `~astropy.units.Quantity` + Radial vector from the surface to the light source. + + ro : `numpy.ndarray` + Radial vector from the surface to the observer. + + + Returns + ------- + radiance : `~astropy.units.Quantity` + Observed radiance. + + """ + + return self.radiance(F_i, *self._vectors_to_angles(n, rs, ro)) diff --git a/sbpy/surfaces/tests/__init__.py b/sbpy/surfaces/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/sbpy/surfaces/tests/test_lambertian.py b/sbpy/surfaces/tests/test_lambertian.py new file mode 100644 index 00000000..685d0699 --- /dev/null +++ b/sbpy/surfaces/tests/test_lambertian.py @@ -0,0 +1,43 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import pytest + +import numpy as np +from astropy.coordinates import Angle +from astropy import units as u + +from ..lambertian import LambertianSurface + + +class TestingSurface(LambertianSurface): + def radiance(): # pragma: no cover + pass + + +@pytest.fixture +def surface(): + return TestingSurface({"albedo": 0.1}) + + +def test_emittance(surface: TestingSurface): + e = Angle([0, 30, 45, 60, 90, 100], "deg") + expected = [1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0] + phi = np.random.rand(len(e)) * 360 * u.deg # independent of phi + result = surface.emittance(e, phi) + assert u.allclose(result, expected) + + +def test_absorptance(surface: TestingSurface): + i = Angle([0, 30, 45, 60, 90, 100], "deg") + expected = 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) + result = surface.absorptance(i) + assert u.allclose(result, expected) + + +def test_reflectance(surface: TestingSurface): + i = Angle([0, 30, 45, 60, 90, 0, 90, 100] * u.deg) + e = Angle([0, 60, 45, 30, 0, 90, 90, 0] * u.deg) + phi = np.random.rand(len(i)) * 360 * u.deg # independent of phi + result = surface.reflectance(i, e, phi) + expected = 0.1 * np.array([1, np.sqrt(3) / 4, 1 / 2, np.sqrt(3) / 4, 0, 0, 0, 0]) + assert u.allclose(result, expected) diff --git a/sbpy/surfaces/tests/test_scattered.py b/sbpy/surfaces/tests/test_scattered.py new file mode 100644 index 00000000..98f0b278 --- /dev/null +++ b/sbpy/surfaces/tests/test_scattered.py @@ -0,0 +1,31 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import numpy as np +from astropy import units as u + +from ...calib import Sun, solar_spectrum +from ..scattered import LambertianSurfaceScatteredSunlight + + +def test_radiance_from_vectors(): + # also tests radiance() + + surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) + + # fake an easy solar spectrum for testing + wave = [0.5, 0.55, 0.6] * u.um + spec = [0.5, 1.0, 1.5] * u.W / (u.m**2 * u.um) + with solar_spectrum.set(Sun.from_array(wave, spec)): + n = np.array([1, 0, 0]) + rs = [1, 1, 0] * u.au + ro = [1, -1, 0] * u.au + + # albedo * F_i / rh**2 * cos(45)**2 + # 0.1 * 1 / 2 / 2 + expected = 0.025 * u.W / (u.m**2 * u.um * u.sr) + result = surface.radiance_from_vectors(0.55 * u.um, n, rs, ro) + assert u.isclose(result, expected) + + # again to test branching to Sun.observe + result = surface.radiance_from_vectors((0.549, 0.55, 0.551) * u.um, n, rs, ro) + assert u.allclose(result[1], expected, rtol=0.01) diff --git a/sbpy/surfaces/tests/test_surface.py b/sbpy/surfaces/tests/test_surface.py new file mode 100644 index 00000000..1299b655 --- /dev/null +++ b/sbpy/surfaces/tests/test_surface.py @@ -0,0 +1,81 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import numpy as np +from astropy.coordinates import Angle +from astropy import units as u + +from ..surface import Surface + + +class TestingSurface(Surface): + def absorptance(self, i): # pragma: no cover + pass + + def emittance(self, e, phi): # pragma: no cover + pass + + def reflectance(self, i: Angle, e: Angle, phi: Angle) -> u.Quantity: + return np.cos(i) * np.cos(e) * np.sin(phi) + + def radiance(self, F_i, i, e, phi): + return F_i * self.reflectance(i, e, phi) / u.sr + + +def test_cos(): + a = Angle([0, 30, 45, 60, 90], "deg") + result = Surface._cos(a) + expected = [1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0] + assert np.allclose(result, expected) + + assert np.isclose(Surface._cos(90 * u.deg), 0) + + +def test_radiance(): + surface = TestingSurface({}) + F_i = 1664 * u.W / (u.m**2 * u.um) + result = surface.radiance(F_i, 30 * u.deg, 30 * u.deg, 60 * u.deg) + assert u.isclose(result, F_i * np.sqrt(27) / 8 / u.sr) + + +def test_vectors_to_angles(): + n = [1, 0, 0] + rs = [1, 1, 0] * u.au + ro = [1, -1, 0] * u.au + i, e, phi = Surface._vectors_to_angles(n, rs, ro) + assert u.isclose(i, 45 * u.deg) + assert u.isclose(e, 45 * u.deg) + assert u.isclose(phi, 90 * u.deg) + + n = [1, 0, 0] + rs = [1 / np.sqrt(3), 1, 0] * u.au + ro = [np.sqrt(3), -1, 0] * u.au + i, e, phi = Surface._vectors_to_angles(n, rs, ro) + assert u.isclose(i, 60 * u.deg) + assert u.isclose(e, 30 * u.deg) + assert u.isclose(phi, 90 * u.deg) + + ro = [1 / np.sqrt(3), -1, 0] * u.au + i, e, phi = Surface._vectors_to_angles(n, rs, ro) + assert u.isclose(e, 60 * u.deg) + assert u.isclose(phi, 120 * u.deg) + + +def test_radiance_from_vectors(): + F_i = 1664 * u.W / (u.m**2 * u.um) + surface = TestingSurface({}) + + n = [1, 0, 0] + rs = [1, 1, 0] * u.au + ro = [1, -1, 0] * u.au + result = surface.radiance_from_vectors(F_i, n, rs, ro) + assert u.isclose(result, F_i / 2 / u.sr) + + n = [1, 0, 0] + rs = [1 / np.sqrt(3), 1, 0] * u.au + ro = [np.sqrt(3), -1, 0] * u.au + result = surface.radiance_from_vectors(F_i, n, rs, ro) + assert u.isclose(result, F_i * np.sqrt(3) / 4 / u.sr) + + ro = [1 / np.sqrt(3), -1, 0] * u.au + result = surface.radiance_from_vectors(F_i, n, rs, ro) + assert u.isclose(result, F_i * np.sqrt(3) / 8 / u.sr) diff --git a/sbpy/surfaces/thermal.py b/sbpy/surfaces/thermal.py new file mode 100644 index 00000000..ea5d3042 --- /dev/null +++ b/sbpy/surfaces/thermal.py @@ -0,0 +1,104 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +from abc import abstractmethod + +import numpy as np +from numpy import pi +from astropy.coordinates import Angle +import astropy.constants as const +import astropy.units as u + +from .surface import Surface +from ..spectroscopy.sources import BlackbodySource +from ..units.typing import SpectralQuantity, SpectralRadianceQuantity + + +class SurfaceThermalEmission(Surface): + """Abstract base class for surface thermal emission. + + The surface is assumed to be illuminated by sunlight. + + + Parameters + ---------- + phys : Phys + Surface physical parameters: + - albedo + - emissivity + - beaming parameter + + """ + + @abstractmethod + def T(self, i: Angle, rh: u.physical.length) -> u.Quantity: + """Surface temperature given incidence angle and heliocentric distance. + + + Parameters + ---------- + i : Angle + Angle from normal of incident light. + + rh : `~astropy.units.Quantity` + Heliocentric distance. + + + Returns + ------- + T : `~astropy.units.Quantity` + Surface temperature. + + """ + + def radiance(self, wave_freq: SpectralQuantity, n, rs, ro) -> u.Quantity: + """Observed radiance from a surface. + + + Parameters + ---------- + wave_freq : `~astropy.units.Quantity` + Wavelength or frequency at which to evaluate the emission. + + n : `numpy.ndarray` + Surface normal vector. + + rs : `~astropy.units.Quantity` + Radial vector from the surface to the light source. + + ro : `~astropy.units.Quantity` + Radial vector from the surface to the observer. + + + Returns + ------- + radiance : `~astropy.units.Quantity` + Observed radiance. + + """ + + n_hat = np.linalg.norm(n.value) + rs_hat = np.linalg.norm(rs.value) + ro_hat = np.linalg.norm(ro.value) + + i = np.arccos(np.dot(n_hat, rs_hat)) + e = np.arccos(np.dot(n_hat, ro_hat)) + phi = np.arccos(np.dot(rs_hat, ro_hat)) + + rh = np.sqrt(np.dot(rs, rs)) + T = self.T(i, rh) + + epsilon = self.phys["emissivity"] + BB = BlackbodySource(T)(wave_freq) / pi + + return epsilon * BB * self.emittance(e, phi) + + +class InstantaneousThermalEquilibrium(SurfaceThermalEmission): + """Abstract base class for a surface in LTE with sunlight.""" + + def T(self, i: Angle, rh: u.physical.length) -> u.Quantity: + sun = const.L / (4 * pi * rh**2) + epsilon = self.phys["emissivity"] + eta = self.phys["eta"] + T = (self.absorptance(i) * sun / epsilon / eta / const.sigma_sb) ** (1 / 4) + return T From 448dc9a86e8ac007b02601a5d73b0545962f10f0 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 09:44:30 -0500 Subject: [PATCH 03/31] Finish top-level description of surfaces --- docs/sbpy/surfaces.rst | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/docs/sbpy/surfaces.rst b/docs/sbpy/surfaces.rst index fc615992..6833ec4e 100644 --- a/docs/sbpy/surfaces.rst +++ b/docs/sbpy/surfaces.rst @@ -1,7 +1,7 @@ Surfaces Module (`sbpy.surfaces`) ================================= -The ``surfaces`` module describes the interaction of electromagnetic radiation with surfaces. Sbpy uses the :math:`(i, e, \phi)` model (angle of incidence, angle of emittance, and phase angle) to describe how absorptance, emittance, and reflectance vary with incoming and outgoing radiation. It has a flexible system that can incorporate any surface scattering model that can be described with these three angles. However, most users will use the built-in surface models. +The ``surfaces`` module describes the interaction of electromagnetic radiation with surfaces. Sbpy uses the :math:`(i, e, \phi)` model (angle of incidence, angle of emittance, and phase angle) to describe how light scatters and emits light. It has a flexible system that can incorporate any surface scattering model that can be described with these three angles. A few built-in surface models are provided. .. figure:: ../static/scattering-vectors.svg @@ -9,7 +9,9 @@ The ``surfaces`` module describes the interaction of electromagnetic radiation w Sbpy's geometric basis for surface scattering and emission: :math:`n` is the surface normal vector, :math:`r_s` is the radial vector to the light source, and :math:`r_o` is the radial vector to the observer. The angle of incidence (:math:`i`), angle of emittance (:math:`e`), phase angle (:math:`\phi`) are labeled. -A ``Surface`` as methods to +A instance of a ``Surface`` will have methods to calculate absorptance, emittance, and reflectance. A radiance method is used to calculate the observed spectral radiance of a surface given incident light. + +Surfaces are expected to require albedo and/or emissivity. Conventions on which property is used and when should be defined by each class. For example, a surface that only calculates reflectance may only require albedo, but one that calculates thermal emission may use the convention of albedo for absorbed sunlight and emissivity for emitted thermal radiation. Built-in surface models @@ -35,3 +37,6 @@ Create an instance of the ``LambertianSurfaceScatteredSunlight`` model, and calc >>> surface.reflectance(i, e, phi) # doctest: +FLOAT_CMP + +Building your own surface models +-------------------------------- From 6fb73184e8491919314e709de0e6be3ed6c0d4db Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 12:58:48 -0500 Subject: [PATCH 04/31] Alternative approach for ScatteredSunlight --- sbpy/surfaces/scattered.py | 53 ++++++++++++++++++++++++-------------- 1 file changed, 34 insertions(+), 19 deletions(-) diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py index be335249..612a4e1b 100644 --- a/sbpy/surfaces/scattered.py +++ b/sbpy/surfaces/scattered.py @@ -7,14 +7,31 @@ from .surface import Surface from .lambertian import LambertianSurface from ..calib import Sun -from ..units.typing import SpectralQuantity, SpectralRadianceQuantity, UnitLike +from ..units.typing import SpectralQuantity, SpectralFluxDensityQuantity, UnitLike -class ScatteredSunlight(Surface): - """Abstract base class to observe light scattered by a surface illuminated by the Sun.""" +class ScatteredLight(Surface): + """Abstract base class to observe light scattered by a surface.""" @u.quantity_input def radiance( + self, + F_i: SpectralFluxDensityQuantity, + i: u.physical.angle, + e: u.physical.angle, + phi: u.physical.angle, + ) -> u.Quantity: + """Observed light reflected from a surface.""" + return F_i * self.reflectance(i, e, phi) / u.sr + + radiance.__doc__ += Surface.radiance.__doc__[Surface.radiance.__doc__.index("\n") :] + + +class ScatteredSunlight(ScatteredLight): + """Abstract base class to observe sunlight scattered by a surface.""" + + @u.quantity_input + def sunlight( self, wave_freq: SpectralQuantity, i: u.physical.angle, @@ -23,7 +40,7 @@ def radiance( rh: u.physical.length, unit: UnitLike = "W/(m2 sr um)", ) -> u.Quantity: - """Observed radiance from a surface. + """Radiance from sunlight scattered by a surface. Parameters @@ -32,6 +49,9 @@ def radiance( Wavelength or frequency at which to evaluate the Sun. Arrays are evaluated with `sbpy.calib.core.Sun.observe()`. + rh : `~astropy.units.Quantity` + Heliocentric distance. + i : `~astropy.units.Quantity` Angle from normal of incident light. @@ -41,9 +61,6 @@ def radiance( phi : `~astropy.units.Quantity` Source-target-observer (phase) angle. - rh : `~astropy.units.Quantity` - Heliocentric distance. - unit : `~astropy.units.Unit`, optional Unit of the return value. @@ -56,16 +73,17 @@ def radiance( """ sun = Sun.from_default() + flux_density_unit = u.Unit(unit) * u.sr if wave_freq.size == 1: - F_i = sun(wave_freq) + F_i = sun(wave_freq, unit=flux_density_unit) else: - F_i = sun.observe(wave_freq) + F_i = sun.observe(wave_freq, unit=flux_density_unit) F_i /= rh.to_value("au") ** 2 - return (F_i * self.reflectance(i, e, phi) / u.sr).to(unit) + return self.radiance(F_i, i, e, phi).to(unit) @u.quantity_input - def radiance_from_vectors( + def sunlight_from_vectors( self, wave_freq: SpectralQuantity, n: np.ndarray[3], @@ -73,7 +91,7 @@ def radiance_from_vectors( ro: u.physical.length, unit: UnitLike = "W/(m2 sr um)", ) -> u.Quantity: - """Observed radiance from a surface. + """Observed sunlight reflected from a surface. Parameters @@ -88,12 +106,9 @@ def radiance_from_vectors( rs : `~astropy.units.Quantity` Radial vector from the surface to the light source. - ro : `numpy.ndarray` + ro : `~astropy.units.Quantity` Radial vector from the surface to the observer. - rh : `~astropy.units.Quantity` - Heliocentric distance. - unit : `~astropy.units.Unit`, optional Unit of the return value. @@ -106,9 +121,9 @@ def radiance_from_vectors( """ rh = np.linalg.norm(rs).to("au") i, e, phi = self._vectors_to_angles(n, rs, ro) - return self.radiance(wave_freq, i, e, phi, rh, unit=unit) + return self.sunlight(wave_freq, rh, i, e, phi, unit=unit) - __doc__ += Surface.__doc__[Surface.__doc__.index("\n") + 1 :] + __doc__ += Surface.__doc__[Surface.__doc__.index("\n") + 1 :] # noqa: E203 class LambertianSurfaceScatteredSunlight(LambertianSurface, ScatteredSunlight): @@ -118,4 +133,4 @@ class LambertianSurfaceScatteredSunlight(LambertianSurface, ScatteredSunlight): """ - __doc__ += Surface.__doc__[Surface.__doc__.index("\n") :] + __doc__ += Surface.__doc__[Surface.__doc__.index("\n") :] # noqa: E203 From 05d00539668b345ffeca33915a9ece620e5c7217 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 13:11:31 -0500 Subject: [PATCH 05/31] Test and bugfix ScatteredLight and ScatteredSunlight --- sbpy/surfaces/scattered.py | 2 +- sbpy/surfaces/tests/test_scattered.py | 23 +++++++++++++++++++---- 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py index 612a4e1b..52231b2e 100644 --- a/sbpy/surfaces/scattered.py +++ b/sbpy/surfaces/scattered.py @@ -34,10 +34,10 @@ class ScatteredSunlight(ScatteredLight): def sunlight( self, wave_freq: SpectralQuantity, + rh: u.physical.length, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle, - rh: u.physical.length, unit: UnitLike = "W/(m2 sr um)", ) -> u.Quantity: """Radiance from sunlight scattered by a surface. diff --git a/sbpy/surfaces/tests/test_scattered.py b/sbpy/surfaces/tests/test_scattered.py index 98f0b278..1211af85 100644 --- a/sbpy/surfaces/tests/test_scattered.py +++ b/sbpy/surfaces/tests/test_scattered.py @@ -7,8 +7,23 @@ from ..scattered import LambertianSurfaceScatteredSunlight -def test_radiance_from_vectors(): - # also tests radiance() +def test_scattered_light(): + + surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) + + F_i = 1 * u.Unit("W/(m2 um)") + n = np.array([1, 0, 0]) + rs = [1, 1, 0] * u.au + ro = [1, -1, 0] * u.au + + # albedo * F_i / rh**2 * cos(45)**2 + # 0.1 * 1 / 2 + expected = 0.05 * u.W / (u.m**2 * u.um * u.sr) + result = surface.radiance_from_vectors(F_i, n, rs, ro) + assert u.isclose(result, expected) + + +def test_scattered_sunlight(): surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) @@ -23,9 +38,9 @@ def test_radiance_from_vectors(): # albedo * F_i / rh**2 * cos(45)**2 # 0.1 * 1 / 2 / 2 expected = 0.025 * u.W / (u.m**2 * u.um * u.sr) - result = surface.radiance_from_vectors(0.55 * u.um, n, rs, ro) + result = surface.sunlight_from_vectors(0.55 * u.um, n, rs, ro) assert u.isclose(result, expected) # again to test branching to Sun.observe - result = surface.radiance_from_vectors((0.549, 0.55, 0.551) * u.um, n, rs, ro) + result = surface.sunlight_from_vectors((0.549, 0.55, 0.551) * u.um, n, rs, ro) assert u.allclose(result[1], expected, rtol=0.01) From db280e6642c9f4fe2182c30e82499ba197c3a68d Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 13:16:10 -0500 Subject: [PATCH 06/31] Rename ScatteredSunlight methods --- sbpy/surfaces/scattered.py | 6 +++--- sbpy/surfaces/tests/test_scattered.py | 6 ++++-- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py index 52231b2e..61bf0d6d 100644 --- a/sbpy/surfaces/scattered.py +++ b/sbpy/surfaces/scattered.py @@ -31,7 +31,7 @@ class ScatteredSunlight(ScatteredLight): """Abstract base class to observe sunlight scattered by a surface.""" @u.quantity_input - def sunlight( + def scattered_sunlight( self, wave_freq: SpectralQuantity, rh: u.physical.length, @@ -83,7 +83,7 @@ def sunlight( return self.radiance(F_i, i, e, phi).to(unit) @u.quantity_input - def sunlight_from_vectors( + def scattered_sunlight_from_vectors( self, wave_freq: SpectralQuantity, n: np.ndarray[3], @@ -121,7 +121,7 @@ def sunlight_from_vectors( """ rh = np.linalg.norm(rs).to("au") i, e, phi = self._vectors_to_angles(n, rs, ro) - return self.sunlight(wave_freq, rh, i, e, phi, unit=unit) + return self.scattered_sunlight(wave_freq, rh, i, e, phi, unit=unit) __doc__ += Surface.__doc__[Surface.__doc__.index("\n") + 1 :] # noqa: E203 diff --git a/sbpy/surfaces/tests/test_scattered.py b/sbpy/surfaces/tests/test_scattered.py index 1211af85..96a211a9 100644 --- a/sbpy/surfaces/tests/test_scattered.py +++ b/sbpy/surfaces/tests/test_scattered.py @@ -38,9 +38,11 @@ def test_scattered_sunlight(): # albedo * F_i / rh**2 * cos(45)**2 # 0.1 * 1 / 2 / 2 expected = 0.025 * u.W / (u.m**2 * u.um * u.sr) - result = surface.sunlight_from_vectors(0.55 * u.um, n, rs, ro) + result = surface.scattered_sunlight_from_vectors(0.55 * u.um, n, rs, ro) assert u.isclose(result, expected) # again to test branching to Sun.observe - result = surface.sunlight_from_vectors((0.549, 0.55, 0.551) * u.um, n, rs, ro) + result = surface.scattered_sunlight_from_vectors( + (0.549, 0.55, 0.551) * u.um, n, rs, ro + ) assert u.allclose(result[1], expected, rtol=0.01) From 848d1d1ccd831ba8ab616ed828cb0b9d3a250b7c Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 13:19:40 -0500 Subject: [PATCH 07/31] Address black and codestyle incompatilibities --- sbpy/surfaces/scattered.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py index 61bf0d6d..b6ea8664 100644 --- a/sbpy/surfaces/scattered.py +++ b/sbpy/surfaces/scattered.py @@ -24,7 +24,9 @@ def radiance( """Observed light reflected from a surface.""" return F_i * self.reflectance(i, e, phi) / u.sr - radiance.__doc__ += Surface.radiance.__doc__[Surface.radiance.__doc__.index("\n") :] + radiance.__doc__ += Surface.radiance.__doc__[ + Surface.radiance.__doc__.index("\n") : # noqa: E203 + ] class ScatteredSunlight(ScatteredLight): From fc01a0b8a6a95967ac0a2a280418dde221f57f32 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 13:24:15 -0500 Subject: [PATCH 08/31] sunlight tests need synphot --- sbpy/surfaces/tests/test_scattered.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/sbpy/surfaces/tests/test_scattered.py b/sbpy/surfaces/tests/test_scattered.py index 96a211a9..f40e9f6b 100644 --- a/sbpy/surfaces/tests/test_scattered.py +++ b/sbpy/surfaces/tests/test_scattered.py @@ -1,5 +1,6 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst +import pytest import numpy as np from astropy import units as u @@ -24,6 +25,7 @@ def test_scattered_light(): def test_scattered_sunlight(): + pytest.importorskip("synphot") surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) From 7d08c0ebe2c28572f7a2167208560bcb8bbcb069 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 10:46:51 -0500 Subject: [PATCH 09/31] Replace Surface._cos with more useful Surface._min_zero_cos --- sbpy/surfaces/lambertian.py | 11 ++++++----- sbpy/surfaces/surface.py | 10 +++++----- sbpy/surfaces/tests/test_surface.py | 12 +++++++----- 3 files changed, 18 insertions(+), 15 deletions(-) diff --git a/sbpy/surfaces/lambertian.py b/sbpy/surfaces/lambertian.py index 6b1d4e1d..313a8aaa 100644 --- a/sbpy/surfaces/lambertian.py +++ b/sbpy/surfaces/lambertian.py @@ -1,7 +1,5 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst -import numpy as np - import astropy.units as u from .surface import Surface @@ -12,16 +10,19 @@ class LambertianSurface(Surface): @u.quantity_input def absorptance(self, i: u.physical.angle) -> u.Quantity: - cos_i = np.maximum(self._cos(i), np.zeros_like(i.value)) + # use self._min_zero_cos(i) which ensures cos(>= 90 deg) = 0 + cos_i = self._min_zero_cos(i) return (1 - self.phys["albedo"]) * cos_i @u.quantity_input def emittance(self, e: u.physical.angle, phi: u.physical.angle) -> u.Quantity: - return np.maximum(self._cos(e), np.zeros_like(e.value)) + # use self._min_zero_cos(e) which ensures cos(>= 90 deg) = 0 + return self._min_zero_cos(e) @u.quantity_input def reflectance( self, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle ) -> u.Quantity: - cos_i = np.maximum(self._cos(i), np.zeros_like(i.value)) + # use self._min_zero_cos(e) which ensures cos(>= 90 deg) = 0 + cos_i = self._min_zero_cos(i) return self.phys["albedo"] * cos_i * self.emittance(e, phi) diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index c6cf6101..4eabaae1 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -32,17 +32,17 @@ def __init__(self, phys: Phys): self.phys = phys @staticmethod - def _cos(a: u.physical.angle) -> u.Quantity: - """Use to ensure that cos(90 deg) equals 0.""" + def _min_zero_cos(a: u.physical.angle) -> u.Quantity: + """Use to ensure that cos(>=90 deg) equals 0.""" # handle scalars separately - if a.ndim == 0 and u.isclose(a, 90 * u.deg): + if a.ndim == 0 and u.isclose(np.abs(a), 90 * u.deg): return u.Quantity(0) x = np.cos(a) - x[u.isclose(a, 90 * u.deg)] = 0 + x[u.isclose(np.abs(a), 90 * u.deg)] = 0 - return x + return np.maximum(x, 0) @abstractmethod def absorptance(self, i: u.physical.angle) -> u.Quantity: diff --git a/sbpy/surfaces/tests/test_surface.py b/sbpy/surfaces/tests/test_surface.py index 1299b655..3a9bff6e 100644 --- a/sbpy/surfaces/tests/test_surface.py +++ b/sbpy/surfaces/tests/test_surface.py @@ -21,13 +21,15 @@ def radiance(self, F_i, i, e, phi): return F_i * self.reflectance(i, e, phi) / u.sr -def test_cos(): - a = Angle([0, 30, 45, 60, 90], "deg") - result = Surface._cos(a) - expected = [1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0] +def test_min_zero_cos(): + a = Angle([-91, -90, 0, 30, 45, 60, 90, 91], "deg") + result = Surface._min_zero_cos(a) + expected = [0, 0, 1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0] assert np.allclose(result, expected) - assert np.isclose(Surface._cos(90 * u.deg), 0) + # test scalars + for i in range(len(a)): + assert np.isclose(Surface._min_zero_cos(a[i]), expected[i]) def test_radiance(): From 5322866357f3e090b39c80458692f1ad1e398895 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:38:31 -0500 Subject: [PATCH 10/31] surface_brightness only available for astropy >= 6.1 --- sbpy/units/typing.py | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/sbpy/units/typing.py b/sbpy/units/typing.py index bd65c71f..e06b3962 100644 --- a/sbpy/units/typing.py +++ b/sbpy/units/typing.py @@ -2,13 +2,27 @@ """sbpy unit and quantity typing.""" from typing import Union +from packaging.version import Version + import astropy.units as u +from astropy import __version__ as astropy_version UnitLike = Union[str, u.Unit] -SpectralQuantity = Union[u.Quantity["length"], u.Quantity["frequency"]] -SpectralFluxDensityQuantity = Union[ - u.Quantity["spectral_flux_density"], u.Quantity["spectral_flux_density_wav"] +SpectralQuantity = Union[ + u.Quantity[u.physical.length], u.Quantity[u.physical.frequency] ] -SpectralRadianceQuantity = Union[ - u.Quantity["surface_brightness"], u.Quantity["surface_brightness_wav"] +SpectralFluxDensityQuantity = Union[ + u.Quantity[u.physical.spectral_flux_density], + u.Quantity[u.physical.spectral_flux_density_wav], ] + +if Version(astropy_version) < Version("6.1"): + SpectralRadianceQuantity = Union[ + u.Quantity[u.physical.spectral_flux_density / u.sr], + u.Quantity[u.physical.spectral_flux_density_wav / u.sr], + ] +else: + SpectralRadianceQuantity = Union[ + u.Quantity[u.physical.surface_brightness], + u.Quantity[u.physical.surface_brightness_wav], + ] From 6498d55e79d107a59143717808365d5f06f2faf7 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:39:15 -0500 Subject: [PATCH 11/31] Remove unused imports. --- sbpy/surfaces/surface.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index 4eabaae1..51d895e0 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -8,12 +8,7 @@ from ..data.phys import Phys from ..data.decorators import dataclass_input -from ..units.typing import ( - SpectralFluxDensityQuantity, - SpectralQuantity, - SpectralRadianceQuantity, - UnitLike, -) +from ..units.typing import SpectralFluxDensityQuantity class Surface(ABC): From 999a38133dabf1254a2da9ed41def77e1c042fa0 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:39:31 -0500 Subject: [PATCH 12/31] Import classes into module __init__ --- sbpy/surfaces/__init__.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/sbpy/surfaces/__init__.py b/sbpy/surfaces/__init__.py index e69de29b..2240d50f 100644 --- a/sbpy/surfaces/__init__.py +++ b/sbpy/surfaces/__init__.py @@ -0,0 +1,7 @@ +from .surface import Surface +from .lambertian import LambertianSurface +from .scattered import ( + ScatteredLight, + ScatteredSunlight, + LambertianSurfaceScatteredSunlight, +) From 8070ecfdd356110c98e0c9fbdd6c42a7fd12e0ee Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:40:03 -0500 Subject: [PATCH 13/31] Complete documentaion. --- docs/sbpy/surfaces.rst | 127 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 121 insertions(+), 6 deletions(-) diff --git a/docs/sbpy/surfaces.rst b/docs/sbpy/surfaces.rst index 6833ec4e..7313b621 100644 --- a/docs/sbpy/surfaces.rst +++ b/docs/sbpy/surfaces.rst @@ -3,7 +3,6 @@ Surfaces Module (`sbpy.surfaces`) The ``surfaces`` module describes the interaction of electromagnetic radiation with surfaces. Sbpy uses the :math:`(i, e, \phi)` model (angle of incidence, angle of emittance, and phase angle) to describe how light scatters and emits light. It has a flexible system that can incorporate any surface scattering model that can be described with these three angles. A few built-in surface models are provided. - .. figure:: ../static/scattering-vectors.svg :alt: Diagram of surface scattering and emission vectors @@ -17,15 +16,13 @@ Surfaces are expected to require albedo and/or emissivity. Conventions on which Built-in surface models ----------------------- -The model `sbpy.surfaces.scattered.LambertianSurfaceScatteredSunlight` is used to observe sunlight scattered from a Lambertian surface (light scattered uniformly in all directions). - -Create an instance of the ``LambertianSurfaceScatteredSunlight`` model, and calculate the absorptance, emittance, and reflectance for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)`. +The model `~sbpy.surfaces.scattered.LambertianSurfaceScatteredSunlight` is used to observe sunlight scattered from a Lambertian surface (light scattered uniformly in all directions). -.. code:: python +Create an instance of the ``LambertianSurfaceScatteredSunlight`` model, and calculate the absorptance, emittance, and reflectance for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)`:: >>> import astropy.units as u >>> import matplotlib.pyplot as plt - >>> from sbpy.surfaces.scattered import LambertianSurfaceScatteredSunlight + >>> from sbpy.surfaces import LambertianSurfaceScatteredSunlight >>> >>> surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) >>> @@ -38,5 +35,123 @@ Create an instance of the ``LambertianSurfaceScatteredSunlight`` model, and calc +Radiance of scattered/emitted light +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +All `~sbpy.surfaces.surface.Surface` models have a :meth:`~sbpy.surfaces.surface.Surface.radiance` method that calculates the radiance of scattered/emitted light, given the spectral flux density of incident light at the surface:: + + >>> F_i = 1000 * u.W / u.m**2 / u.um # incident spectral flux density + >>> surface.radiance(F_i, i, e, phi) # doctest: +FLOAT_CMP + + + +Radiance of scattered sunlight +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +``LambertianSurfaceScatteredSunlight`` is derived from the ``ScatteredSunlight`` class, which provides convenience methods for calculating the radiance of sunlight scattered off the surface:: + + >>> wave = 0.55 * u.um + >>> rh = 1 * u.au # heliocentric distance of the surface + >>> surface.scattered_sunlight(wave, rh, i, e, phi) # doctest: +FLOAT_CMP + + +The solar spectrum is configured using Sbpy's calibration system. See :doc:`calib` for details. + + +Radiance from vectors +^^^^^^^^^^^^^^^^^^^^^ + +As an alternative to using :math:`(i, e, \phi)`, radiance may be calculated using vectors that define the normal direction, radial vector of the light source, and radial vector of the observer:: + + >>> # the following vectors are equivalent to (i, e, phi) = (30, 60, 90) deg + >>> n = [1, 0, 0] + >>> rs = [0.866, 0.5, 0] * u.au + >>> ro = [0.5, -0.866, 0] * u.au + >>> + >>> surface.radiance_from_vectors(F_i, n, rs, ro) # doctest: +FLOAT_CMP + + >>> surface.scattered_sunlight_from_vectors(wave, n, rs, ro) # doctest: +FLOAT_CMP + + +Notice that heliocentric distance was not needed in the call to ``scattered_sunlight_from_vectors`` because it was already accounted for by the ``rs`` vector. + + Building your own surface models -------------------------------- + +Defining your own surface model is typically done by creating a new class based on `~sbpy.surfaces.surface.Surface`, and defining methods for ``absorptance``, ``emittance``, and ``reflectance``. The `~sbpy.surfaces.lambertian.Lambertian` model serves as a good example. For surface scattering problems, most users will combine their class with the `~sbpy.surfaces.scattered.ScatteredLight` or `~sbpy.surfaces.scattered.ScatteredSunlight` classes, which provide the ``radiance`` method to complete the ``Surface`` model. + +Here, we define a new surface model with surface scattering proportional to :math:`\cos^2` based on the `~sbpy.surfaces.scattered.ScatteredSunlight` class:: + + >>> import numpy as np + >>> from sbpy.surfaces import Surface, ScatteredSunlight + >>> + >>> class Cos2SurfaceScatteredSunlight(ScatteredSunlight): + ... """Absorption and scattering proportional to cos**2.""" + ... + ... def absorptance(self, i): + ... return (1 - self.phys["albedo"]) * np.cos(i)**2 + ... + ... def emittance(self, e, phi): + ... return np.cos(e)**2 + ... + ... def reflectance(self, i, e, phi): + ... return self.phys["albedo"] * np.cos(i)**2 * self.emittance(e, phi) + >>> + >>> surface = Cos2SurfaceScatteredSunlight({"albedo": 0.1}) + >>> surface.reflectance(i, e, phi) # doctest: +FLOAT_CMP + + >>> surface.radiance(F_i, i, e, phi) # doctest: +FLOAT_CMP + + >>> surface.scattered_sunlight(wave, rh, i, e, phi) # doctest: +FLOAT_CMP + + +However, if a scattering model will be re-used with other classes, e.g., for scattered light and thermal emission modeling, then the most flexible approach is to base the model on ``Surface`` and have derived classes combine the model with scattering or thermal emission classes:: + + >>> class Cos2Surface(Surface): + ... """Abstract base class for absorption and scattering proportional to cos**2. + ... + ... We document this as an "abstract base class" because the ``radiance`` + ... method required by ``Surface`` is not implemented here, and therefore + ... this class cannot be used directly (i.e., it cannot be instantiated). + ... + ... """ + ... + ... def absorptance(self, i): + ... return (1 - self.phys["albedo"]) * np.cos(i)**2 + ... + ... def emittance(self, e, phi): + ... return np.cos(e)**2 + ... + ... def reflectance(self, i, e, phi): + ... return self.phys["albedo"] * np.cos(i)**2 * self.emittance(e, phi) + >>> + >>> class Cos2SurfaceScatteredSunlightAlt(Cos2Surface, ScatteredSunlight): + ... """Absorption and scattering proportional to cos**2. + ... + ... This class combines the ``absorptance``, ``emittance``, and ``reflectance`` + ... methods from ``Cos2Surface`` with the ``radiance`` and ``scattered_sunlight`` + ... methods of ``ScatteredSunlight``. + ... + ... Parameters + ... ... + ... """ + >>> + >>> class Cos2SurfaceThermalEmission(Cos2Surface, ThermalEmission): # doctest: +SKIP + ... """Absorption and thermal emission proportional to cos**2. + ... + ... This class combines the ``absorptance``, ``emittance``, and ``reflectance`` + ... from ``Cos2Surface`` with the ``radiance`` method of a fictitious + ... ``ThermalEmission`` class. + ... + ... Parameters + ... ... + ... """ + +Thanks to their base classes (``Cos2Surface``, ``ScatteredSunlight``, and ``ThermalEmission``), the ``Cos2SurfaceScatteredSunlightAlt`` and ``Cos2SurfaceThermalEmission`` classes are complete and no additional code is needed, just documentation. + +Reference/API +------------- +.. automodapi:: sbpy.surfaces + :no-heading: + :inherited-members: From 0e55a7f7fe367de74d7690819c106e0bf9537ea0 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:55:08 -0500 Subject: [PATCH 14/31] Mark required packages for new doctests blocks --- docs/sbpy/surfaces.rst | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/docs/sbpy/surfaces.rst b/docs/sbpy/surfaces.rst index 7313b621..de9b9ba1 100644 --- a/docs/sbpy/surfaces.rst +++ b/docs/sbpy/surfaces.rst @@ -21,7 +21,6 @@ The model `~sbpy.surfaces.scattered.LambertianSurfaceScatteredSunlight` is used Create an instance of the ``LambertianSurfaceScatteredSunlight`` model, and calculate the absorptance, emittance, and reflectance for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)`:: >>> import astropy.units as u - >>> import matplotlib.pyplot as plt >>> from sbpy.surfaces import LambertianSurfaceScatteredSunlight >>> >>> surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) @@ -50,6 +49,8 @@ Radiance of scattered sunlight ``LambertianSurfaceScatteredSunlight`` is derived from the ``ScatteredSunlight`` class, which provides convenience methods for calculating the radiance of sunlight scattered off the surface:: +.. doctest-requires:: synphot + >>> wave = 0.55 * u.um >>> rh = 1 * u.au # heliocentric distance of the surface >>> surface.scattered_sunlight(wave, rh, i, e, phi) # doctest: +FLOAT_CMP @@ -63,6 +64,8 @@ Radiance from vectors As an alternative to using :math:`(i, e, \phi)`, radiance may be calculated using vectors that define the normal direction, radial vector of the light source, and radial vector of the observer:: +.. doctest-requires:: synphot + >>> # the following vectors are equivalent to (i, e, phi) = (30, 60, 90) deg >>> n = [1, 0, 0] >>> rs = [0.866, 0.5, 0] * u.au @@ -83,6 +86,8 @@ Defining your own surface model is typically done by creating a new class based Here, we define a new surface model with surface scattering proportional to :math:`\cos^2` based on the `~sbpy.surfaces.scattered.ScatteredSunlight` class:: +.. doctest-requires:: synphot + >>> import numpy as np >>> from sbpy.surfaces import Surface, ScatteredSunlight >>> @@ -106,6 +111,11 @@ Here, we define a new surface model with surface scattering proportional to :mat >>> surface.scattered_sunlight(wave, rh, i, e, phi) # doctest: +FLOAT_CMP +.. for test runs without synphot +.. testsetup:: + + >>> from sbpy.surfaces import Surface, ScatteredSunlight + However, if a scattering model will be re-used with other classes, e.g., for scattered light and thermal emission modeling, then the most flexible approach is to base the model on ``Surface`` and have derived classes combine the model with scattering or thermal emission classes:: >>> class Cos2Surface(Surface): From e2b9aa7ca64b0e3be82e39e8776456dc481cb5b3 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:59:38 -0500 Subject: [PATCH 15/31] Avoid: each t must be a type. Got PhysicalType('angle') --- sbpy/surfaces/surface.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index 51d895e0..08a90707 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -1,6 +1,5 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst -from typing import Tuple from abc import ABC, abstractmethod import numpy as np @@ -153,7 +152,7 @@ def _vectors_to_angles( n: np.ndarray[3], rs: u.physical.length, ro: np.ndarray[3], - ) -> Tuple[u.physical.angle, u.physical.angle, u.physical.angle]: + ) -> tuple: n_hat = n / np.linalg.norm(n) rs_hat = rs / np.linalg.norm(rs) ro_hat = ro / np.linalg.norm(ro) From a74835a1930c247ee278ca5e93f74a9502fd1107 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 13:30:46 -0500 Subject: [PATCH 16/31] Fix incompatilibity with oldestdeps --- sbpy/surfaces/scattered.py | 2 +- sbpy/surfaces/surface.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py index b6ea8664..d23e244a 100644 --- a/sbpy/surfaces/scattered.py +++ b/sbpy/surfaces/scattered.py @@ -88,7 +88,7 @@ def scattered_sunlight( def scattered_sunlight_from_vectors( self, wave_freq: SpectralQuantity, - n: np.ndarray[3], + n: np.ndarray, rs: u.physical.length, ro: u.physical.length, unit: UnitLike = "W/(m2 sr um)", diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index 08a90707..9d437142 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -149,9 +149,9 @@ def radiance( @staticmethod def _vectors_to_angles( - n: np.ndarray[3], + n: np.ndarray, rs: u.physical.length, - ro: np.ndarray[3], + ro: np.ndarray, ) -> tuple: n_hat = n / np.linalg.norm(n) rs_hat = rs / np.linalg.norm(rs) @@ -167,9 +167,9 @@ def _vectors_to_angles( def radiance_from_vectors( self, F_i: SpectralFluxDensityQuantity, - n: np.ndarray[3], + n: np.ndarray, rs: u.physical.length, - ro: np.ndarray[3], + ro: np.ndarray, ) -> u.Quantity: """Observed radiance from a surface with geometry defined by vectors. From 334ce77d214b62ecf60fd091b1296e723c69a282 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 16:23:31 -0500 Subject: [PATCH 17/31] Thermal was not supposed to be checked in. --- sbpy/surfaces/thermal.py | 104 --------------------------------------- 1 file changed, 104 deletions(-) delete mode 100644 sbpy/surfaces/thermal.py diff --git a/sbpy/surfaces/thermal.py b/sbpy/surfaces/thermal.py deleted file mode 100644 index ea5d3042..00000000 --- a/sbpy/surfaces/thermal.py +++ /dev/null @@ -1,104 +0,0 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst - -from abc import abstractmethod - -import numpy as np -from numpy import pi -from astropy.coordinates import Angle -import astropy.constants as const -import astropy.units as u - -from .surface import Surface -from ..spectroscopy.sources import BlackbodySource -from ..units.typing import SpectralQuantity, SpectralRadianceQuantity - - -class SurfaceThermalEmission(Surface): - """Abstract base class for surface thermal emission. - - The surface is assumed to be illuminated by sunlight. - - - Parameters - ---------- - phys : Phys - Surface physical parameters: - - albedo - - emissivity - - beaming parameter - - """ - - @abstractmethod - def T(self, i: Angle, rh: u.physical.length) -> u.Quantity: - """Surface temperature given incidence angle and heliocentric distance. - - - Parameters - ---------- - i : Angle - Angle from normal of incident light. - - rh : `~astropy.units.Quantity` - Heliocentric distance. - - - Returns - ------- - T : `~astropy.units.Quantity` - Surface temperature. - - """ - - def radiance(self, wave_freq: SpectralQuantity, n, rs, ro) -> u.Quantity: - """Observed radiance from a surface. - - - Parameters - ---------- - wave_freq : `~astropy.units.Quantity` - Wavelength or frequency at which to evaluate the emission. - - n : `numpy.ndarray` - Surface normal vector. - - rs : `~astropy.units.Quantity` - Radial vector from the surface to the light source. - - ro : `~astropy.units.Quantity` - Radial vector from the surface to the observer. - - - Returns - ------- - radiance : `~astropy.units.Quantity` - Observed radiance. - - """ - - n_hat = np.linalg.norm(n.value) - rs_hat = np.linalg.norm(rs.value) - ro_hat = np.linalg.norm(ro.value) - - i = np.arccos(np.dot(n_hat, rs_hat)) - e = np.arccos(np.dot(n_hat, ro_hat)) - phi = np.arccos(np.dot(rs_hat, ro_hat)) - - rh = np.sqrt(np.dot(rs, rs)) - T = self.T(i, rh) - - epsilon = self.phys["emissivity"] - BB = BlackbodySource(T)(wave_freq) / pi - - return epsilon * BB * self.emittance(e, phi) - - -class InstantaneousThermalEquilibrium(SurfaceThermalEmission): - """Abstract base class for a surface in LTE with sunlight.""" - - def T(self, i: Angle, rh: u.physical.length) -> u.Quantity: - sun = const.L / (4 * pi * rh**2) - epsilon = self.phys["emissivity"] - eta = self.phys["eta"] - T = (self.absorptance(i) * sun / epsilon / eta / const.sigma_sb) ** (1 / 4) - return T From 4441a94853ae32c998f84df349910fceadef5c77 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 26 Nov 2024 22:49:03 -0500 Subject: [PATCH 18/31] Clarify as bidirectional reflectance --- sbpy/surfaces/surface.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index 9d437142..48b98ffa 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -87,12 +87,13 @@ def emittance(self, e: u.physical.angle, phi: u.physical.angle) -> u.Quantity: def reflectance( self, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle ) -> u.Quantity: - r"""Reflectance. + r"""Bidirectional reflectance. - The surface is illuminated by incident flux density, :math:`F_i`, at an - angle of :math:`i`, and emitted toward an angle of :math:`e`, measured - from the surface normal direction. :math:`\phi` is the - source-target-observer (phase) angle. + The surface is illuminated by incident flux density (irradiance), + :math:`F_i`, at an angle of :math:`i`, and emitted toward an angle of + :math:`e`, measured from the surface normal direction. :math:`\phi` is + the source-target-observer (phase) angle. Both the source and the + emitted light are assumed to be collimated. Parameters From 9f69b2f292fda1ce7688cef8a9b53a38c1c5a643 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 10 Dec 2024 09:50:03 -0500 Subject: [PATCH 19/31] Partial updates for new method names. --- docs/sbpy/surfaces.rst | 5 +++ sbpy/surfaces/lambertian.py | 35 ++++++++++++----- sbpy/surfaces/scattered.py | 23 +---------- sbpy/surfaces/surface.py | 78 +++++++++++++++---------------------- 4 files changed, 65 insertions(+), 76 deletions(-) diff --git a/docs/sbpy/surfaces.rst b/docs/sbpy/surfaces.rst index de9b9ba1..78e390b6 100644 --- a/docs/sbpy/surfaces.rst +++ b/docs/sbpy/surfaces.rst @@ -1,6 +1,11 @@ Surfaces Module (`sbpy.surfaces`) ================================= +.. admonition:: warning + + The surface module is being made available on a preview basis. The API is + subject to change. Feedback on the approach used is welcome. + The ``surfaces`` module describes the interaction of electromagnetic radiation with surfaces. Sbpy uses the :math:`(i, e, \phi)` model (angle of incidence, angle of emittance, and phase angle) to describe how light scatters and emits light. It has a flexible system that can incorporate any surface scattering model that can be described with these three angles. A few built-in surface models are provided. .. figure:: ../static/scattering-vectors.svg diff --git a/sbpy/surfaces/lambertian.py b/sbpy/surfaces/lambertian.py index 313a8aaa..88adeff1 100644 --- a/sbpy/surfaces/lambertian.py +++ b/sbpy/surfaces/lambertian.py @@ -1,28 +1,45 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst +import numpy as np import astropy.units as u from .surface import Surface +from ..units.typing import SpectralFluxDensityQuantity class LambertianSurface(Surface): """Abstract base class for Lambertian surfaces.""" @u.quantity_input - def absorptance(self, i: u.physical.angle) -> u.Quantity: - # use self._min_zero_cos(i) which ensures cos(>= 90 deg) = 0 + def absorption( + self, + F_i: SpectralFluxDensityQuantity, + i: u.physical.angle, + ) -> u.Quantity: + # use self._min_zero_cos(i) to ensure cos(>= 90 deg) = 0 cos_i = self._min_zero_cos(i) - return (1 - self.phys["albedo"]) * cos_i + return F_i * (1 - self.phys["albedo"]) * cos_i @u.quantity_input - def emittance(self, e: u.physical.angle, phi: u.physical.angle) -> u.Quantity: - # use self._min_zero_cos(e) which ensures cos(>= 90 deg) = 0 - return self._min_zero_cos(e) + def emission( + self, + X_e: SpectralFluxDensityQuantity, + e: u.physical.angle, + phi: u.physical.angle, + ) -> u.Quantity: + # use self._min_zero_cos(e) to ensure cos(>= 90 deg) = 0 + cos_e = self._min_zero_cos(e) + return X_e * cos_e / np.pi / u.sr @u.quantity_input def reflectance( - self, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle + self, + F_i: SpectralFluxDensityQuantity, + i: u.physical.angle, + e: u.physical.angle, + phi: u.physical.angle, ) -> u.Quantity: - # use self._min_zero_cos(e) which ensures cos(>= 90 deg) = 0 + # use self._min_zero_cos(theta) to ensure cos(>= 90 deg) = 0 cos_i = self._min_zero_cos(i) - return self.phys["albedo"] * cos_i * self.emittance(e, phi) + cos_e = self._min_zero_cos(e) + return F_i * self.phys["albedo"] * cos_i * cos_e / np.pi diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py index d23e244a..9cfc3753 100644 --- a/sbpy/surfaces/scattered.py +++ b/sbpy/surfaces/scattered.py @@ -10,26 +10,7 @@ from ..units.typing import SpectralQuantity, SpectralFluxDensityQuantity, UnitLike -class ScatteredLight(Surface): - """Abstract base class to observe light scattered by a surface.""" - - @u.quantity_input - def radiance( - self, - F_i: SpectralFluxDensityQuantity, - i: u.physical.angle, - e: u.physical.angle, - phi: u.physical.angle, - ) -> u.Quantity: - """Observed light reflected from a surface.""" - return F_i * self.reflectance(i, e, phi) / u.sr - - radiance.__doc__ += Surface.radiance.__doc__[ - Surface.radiance.__doc__.index("\n") : # noqa: E203 - ] - - -class ScatteredSunlight(ScatteredLight): +class ScatteredSunlight(Surface): """Abstract base class to observe sunlight scattered by a surface.""" @u.quantity_input @@ -82,7 +63,7 @@ def scattered_sunlight( F_i = sun.observe(wave_freq, unit=flux_density_unit) F_i /= rh.to_value("au") ** 2 - return self.radiance(F_i, i, e, phi).to(unit) + return self.reflectance(F_i, i, e, phi).to(unit) @u.quantity_input def scattered_sunlight_from_vectors( diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index 48b98ffa..a6196a5a 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -39,8 +39,12 @@ def _min_zero_cos(a: u.physical.angle) -> u.Quantity: return np.maximum(x, 0) @abstractmethod - def absorptance(self, i: u.physical.angle) -> u.Quantity: - r"""Absorption of incident light. + def absorption( + self, + F_i: SpectralFluxDensityQuantity, + i: u.physical.angle, + ) -> u.Quantity: + r"""Absorption of directional, incident light. The surface is illuminated by incident flux density, :math:`F_i`, at an angle of :math:`i`, measured from the surface normal direction. @@ -48,20 +52,28 @@ def absorptance(self, i: u.physical.angle) -> u.Quantity: Parameters ---------- + F_i : `astropy.units.Quantity` + Incident light (spectral flux density / spectral irradiance). + i : `~astropy.units.Quantity` Angle from normal of incident light. Returns ------- - a : `~astropy.units.Quantity` - Surface absorptance. + F_a : `~astropy.units.Quantity` + Absorbed spectral flux density. """ @abstractmethod - def emittance(self, e: u.physical.angle, phi: u.physical.angle) -> u.Quantity: - r"""Emission of light. + def emission( + self, + X_e: SpectralFluxDensityQuantity, + e: u.physical.angle, + phi: u.physical.angle, + ) -> u.Quantity: + r"""Emission of light from a surface, as seen by a distant observer. The surface is observed at an angle of :math:`e`, measured from the surface normal direction, and at a solar phase angle of :math:`phi`. @@ -69,8 +81,11 @@ def emittance(self, e: u.physical.angle, phi: u.physical.angle) -> u.Quantity: Parameters ---------- + X_e : `astropy.units.Quantity` + Emitted spectral radiance. + e : `~astropy.units.Quantity` - Angle from normal of emitted light. + Observed angle from normal. phi : `~astropy.units.Quantity` Source-target-observer (phase) angle. @@ -78,8 +93,9 @@ def emittance(self, e: u.physical.angle, phi: u.physical.angle) -> u.Quantity: Returns ------- - e : `~astropy.units.Quantity` - Surface emittance. + F_e : `~astropy.units.Quantity` + Spectral flux density / spectral irradiance received by the + observer. """ @@ -96,36 +112,6 @@ def reflectance( emitted light are assumed to be collimated. - Parameters - ---------- - i : `~astropy.units.Quantity` - Angle from normal of incident light. - - e : `~astropy.units.Quantity` - Angle from normal of emitted light. - - phi : `~astropy.units.Quantity` - Source-target-observer (phase) angle. - - - Returns - ------- - r : `~astropy.units.Quantity` - Surface reflectance. - - """ - - @abstractmethod - def radiance( - self, - F_i: SpectralFluxDensityQuantity, - i: u.physical.angle, - e: u.physical.angle, - phi: u.physical.angle, - ) -> u.Quantity: - """Observed radiance from a surface. - - Parameters ---------- F_i : `astropy.units.Quantity` @@ -143,8 +129,8 @@ def radiance( Returns ------- - radiance : `~astropy.units.Quantity` - Observed radiance. + F_r : `~astropy.units.Quantity` + Spectral flux density / spectral irradiance received by the observer. """ @@ -165,14 +151,14 @@ def _vectors_to_angles( return i, e, phi @u.quantity_input - def radiance_from_vectors( + def reflectance_from_vectors( self, F_i: SpectralFluxDensityQuantity, n: np.ndarray, rs: u.physical.length, ro: np.ndarray, ) -> u.Quantity: - """Observed radiance from a surface with geometry defined by vectors. + """Vector based alternative to reflectance(). Input vectors do not need to be normalized. @@ -194,9 +180,9 @@ def radiance_from_vectors( Returns ------- - radiance : `~astropy.units.Quantity` - Observed radiance. + reflectance : `~astropy.units.Quantity` + Reflectance. """ - return self.radiance(F_i, *self._vectors_to_angles(n, rs, ro)) + return self.reflectance(F_i, *self._vectors_to_angles(n, rs, ro)) From 833e066e10738f439e53473bfb3777cb071cbb0b Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Fri, 5 Sep 2025 09:20:13 -0400 Subject: [PATCH 20/31] Simplify to avoid using instance variables --- sbpy/surfaces/lambertian.py | 33 ++-- sbpy/surfaces/scattered.py | 250 ++++++++++++++++--------- sbpy/surfaces/surface.py | 59 +++--- sbpy/surfaces/tests/test_lambertian.py | 50 ++--- sbpy/surfaces/tests/test_scattered.py | 74 ++++---- 5 files changed, 281 insertions(+), 185 deletions(-) diff --git a/sbpy/surfaces/lambertian.py b/sbpy/surfaces/lambertian.py index 88adeff1..4f11d035 100644 --- a/sbpy/surfaces/lambertian.py +++ b/sbpy/surfaces/lambertian.py @@ -3,43 +3,46 @@ import numpy as np import astropy.units as u -from .surface import Surface +from .surface import Surface, min_zero_cos from ..units.typing import SpectralFluxDensityQuantity class LambertianSurface(Surface): - """Abstract base class for Lambertian surfaces.""" + """Lambertian surface absorption, emission, and reflectance.""" + @staticmethod @u.quantity_input def absorption( - self, F_i: SpectralFluxDensityQuantity, + epsilon: float, i: u.physical.angle, ) -> u.Quantity: - # use self._min_zero_cos(i) to ensure cos(>= 90 deg) = 0 - cos_i = self._min_zero_cos(i) - return F_i * (1 - self.phys["albedo"]) * cos_i + # use min_zero_cos(i) to ensure cos(>= 90 deg) = 0 + cos_i = min_zero_cos(i) + return F_i * epsilon * cos_i + @staticmethod @u.quantity_input def emission( - self, X_e: SpectralFluxDensityQuantity, + epsilon: float, e: u.physical.angle, phi: u.physical.angle, ) -> u.Quantity: - # use self._min_zero_cos(e) to ensure cos(>= 90 deg) = 0 - cos_e = self._min_zero_cos(e) - return X_e * cos_e / np.pi / u.sr + # use min_zero_cos(e) to ensure cos(>= 90 deg) = 0 + cos_e = min_zero_cos(e) + return X_e * epsilon * cos_e / np.pi / u.sr + @staticmethod @u.quantity_input def reflectance( - self, F_i: SpectralFluxDensityQuantity, + albedo: float, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle, ) -> u.Quantity: - # use self._min_zero_cos(theta) to ensure cos(>= 90 deg) = 0 - cos_i = self._min_zero_cos(i) - cos_e = self._min_zero_cos(e) - return F_i * self.phys["albedo"] * cos_i * cos_e / np.pi + # use min_zero_cos(theta) to ensure cos(>= 90 deg) = 0 + cos_i = min_zero_cos(i) + cos_e = min_zero_cos(e) + return F_i * albedo * cos_i * cos_e / np.pi / u.sr diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py index 9cfc3753..14509d9d 100644 --- a/sbpy/surfaces/scattered.py +++ b/sbpy/surfaces/scattered.py @@ -1,119 +1,199 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst +# # Licensed under a 3-clause BSD style license - see LICENSE.rst -import numpy as np +# from abc import ABC, abstractmethod -import astropy.units as u +# import numpy as np -from .surface import Surface -from .lambertian import LambertianSurface -from ..calib import Sun -from ..units.typing import SpectralQuantity, SpectralFluxDensityQuantity, UnitLike +# import astropy.units as u +# from .surface import Surface +# from .lambertian import LambertianSurface +# from ..calib import Sun +# from ..units.typing import SpectralQuantity, SpectralFluxDensityQuantity, UnitLike -class ScatteredSunlight(Surface): - """Abstract base class to observe sunlight scattered by a surface.""" - @u.quantity_input - def scattered_sunlight( - self, - wave_freq: SpectralQuantity, - rh: u.physical.length, - i: u.physical.angle, - e: u.physical.angle, - phi: u.physical.angle, - unit: UnitLike = "W/(m2 sr um)", - ) -> u.Quantity: - """Radiance from sunlight scattered by a surface. +# class ScatteredLight(ABC): +# """Abstract base class to observe light scattered by a surface.""" +# @u.quantity_input +# def scattered_light( +# self, +# wave_freq: SpectralQuantity, +# i: u.physical.angle, +# e: u.physical.angle, +# phi: u.physical.angle, +# unit: UnitLike = "W/(m2 sr um)", +# ) -> u.Quantity: +# """Radiance from light scattered by a surface. - Parameters - ---------- - wave_freq : `astropy.units.Quantity` - Wavelength or frequency at which to evaluate the Sun. Arrays are - evaluated with `sbpy.calib.core.Sun.observe()`. - rh : `~astropy.units.Quantity` - Heliocentric distance. +# Parameters +# ---------- - i : `~astropy.units.Quantity` - Angle from normal of incident light. +# wave_freq : `astropy.units.Quantity` +# Wavelength or frequency at which to evaluate the light source. - e : `~astropy.units.Quantity` - Angle from normal of emitted light. +# i : `~astropy.units.Quantity` +# Angle from normal of incident light. - phi : `~astropy.units.Quantity` - Source-target-observer (phase) angle. +# e : `~astropy.units.Quantity` +# Angle from normal of emitted light. - unit : `~astropy.units.Unit`, optional - Unit of the return value. +# phi : `~astropy.units.Quantity` +# Source-target-observer (phase) angle. +# unit : `~astropy.units.Unit`, optional +# Unit of the return value. - Returns - ------- - radiance : `~astropy.units.Quantity` - Observed radiance. - """ +# Returns +# ------- - sun = Sun.from_default() - flux_density_unit = u.Unit(unit) * u.sr - if wave_freq.size == 1: - F_i = sun(wave_freq, unit=flux_density_unit) - else: - F_i = sun.observe(wave_freq, unit=flux_density_unit) +# radiance : `~astropy.units.Quantity` +# Observed radiance. - F_i /= rh.to_value("au") ** 2 - return self.reflectance(F_i, i, e, phi).to(unit) +# """ - @u.quantity_input - def scattered_sunlight_from_vectors( - self, - wave_freq: SpectralQuantity, - n: np.ndarray, - rs: u.physical.length, - ro: u.physical.length, - unit: UnitLike = "W/(m2 sr um)", - ) -> u.Quantity: - """Observed sunlight reflected from a surface. +# @u.quantity_input +# def scattered_light_from_vectors( +# self, +# wave_freq: SpectralQuantity, +# n: np.ndarray, +# rs: u.physical.length, +# ro: u.physical.length, +# unit: UnitLike = "W/(m2 sr um)", +# ) -> u.Quantity: +# """Observed light reflected from a surface. - Parameters - ---------- - wave_freq : `astropy.units.Quantity` - Wavelength or frequency at which to evaluate the Sun. Arrays are - evaluated with `sbpy.calib.core.Sun.observe()`. +# Parameters +# ---------- - n : `numpy.ndarray` - Surface normal vector. +# wave_freq : `astropy.units.Quantity` +# Wavelength or frequency at which to evaluate the light source. - rs : `~astropy.units.Quantity` - Radial vector from the surface to the light source. +# n : `numpy.ndarray` +# Surface normal vector. - ro : `~astropy.units.Quantity` - Radial vector from the surface to the observer. +# rs : `~astropy.units.Quantity` +# Radial vector from the surface to the light source. - unit : `~astropy.units.Unit`, optional - Unit of the return value. +# ro : `~astropy.units.Quantity` +# Radial vector from the surface to the observer. +# unit : `~astropy.units.Unit`, optional +# Unit of the return value. - Returns - ------- - radiance : `~astropy.units.Quantity` - Observed radiance. - """ - rh = np.linalg.norm(rs).to("au") - i, e, phi = self._vectors_to_angles(n, rs, ro) - return self.scattered_sunlight(wave_freq, rh, i, e, phi, unit=unit) +# Returns +# ------- - __doc__ += Surface.__doc__[Surface.__doc__.index("\n") + 1 :] # noqa: E203 +# radiance : `~astropy.units.Quantity` +# Observed radiance. +# """ -class LambertianSurfaceScatteredSunlight(LambertianSurface, ScatteredSunlight): - """Sunlight scattered from a Lambertian surface. - The surface is assumed to be illuminated by the Sun. +# class ScatteredSunlight(ScatteredLight): +# """Observe sunlight scattered by a surface.""" - """ +# @u.quantity_input +# def scattered_light( +# self, +# wave_freq: SpectralQuantity, +# i: u.physical.angle, +# e: u.physical.angle, +# phi: u.physical.angle, +# rh: u.physical.length = 1 * u.au, +# unit: UnitLike = "W/(m2 sr um)", +# ) -> u.Quantity: +# """Radiance from sunlight scattered by a surface. - __doc__ += Surface.__doc__[Surface.__doc__.index("\n") :] # noqa: E203 + +# Parameters +# ---------- + +# wave_freq : `astropy.units.Quantity` +# Wavelength or frequency at which to evaluate the Sun. Arrays are +# evaluated with `sbpy.calib.core.Sun.observe()`. + +# i : `~astropy.units.Quantity` +# Angle from normal of incident light. + +# e : `~astropy.units.Quantity` +# Angle from normal of emitted light. + +# phi : `~astropy.units.Quantity` +# Source-target-observer (phase) angle. + +# rh : `~astropy.units.Quantity` +# Heliocentric distance, default = 1 au. + +# unit : `~astropy.units.Unit`, optional +# Unit of the return value. + + +# Returns +# ------- +# radiance : `~astropy.units.Quantity` +# Observed radiance. + +# """ + +# sun = Sun.from_default() +# flux_density_unit = u.Unit(unit) * u.sr +# if wave_freq.size == 1: +# F_i = sun(wave_freq, unit=flux_density_unit) +# else: +# F_i = sun.observe(wave_freq, unit=flux_density_unit) + +# F_i /= rh.to_value("au") ** 2 +# return self.reflectance(F_i, i, e, phi).to(unit) + +# @u.quantity_input +# def scattered_light_from_vectors( +# self, +# wave_freq: SpectralQuantity, +# n: np.ndarray, +# rs: u.physical.length, +# ro: u.physical.length, +# unit: UnitLike = "W/(m2 sr um)", +# ) -> u.Quantity: +# """Observed sunlight reflected from a surface. + + +# Parameters +# ---------- + +# wave_freq : `astropy.units.Quantity` +# Wavelength or frequency at which to evaluate the Sun. Arrays are +# evaluated with `sbpy.calib.core.Sun.observe()`. + +# n : `numpy.ndarray` +# Surface normal vector. + +# rs : `~astropy.units.Quantity` +# Radial vector from the surface to the light source. + +# ro : `~astropy.units.Quantity` +# Radial vector from the surface to the observer. + +# unit : `~astropy.units.Unit`, optional +# Unit of the return value. + + +# Returns +# ------- + +# radiance : `~astropy.units.Quantity` +# Observed radiance. + +# """ + +# rh = np.linalg.norm(rs).to("au") +# i, e, phi = self._vectors_to_angles(n, rs, ro) +# return self.scattered_sunlight(wave_freq, i, e, phi, rh=rh, unit=unit) + + +# class LambertianSurfaceScatteredSunlight(LambertianSurface, ScatteredSunlight): +# """Sunlight scattered from a Lambertian surface.""" diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index a6196a5a..3fa74049 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -10,38 +10,27 @@ from ..units.typing import SpectralFluxDensityQuantity -class Surface(ABC): - """Abstract base class for all small-body surfaces. - - - Parameters - ---------- - phys : Phys - Surface physical parameters, e.g., albedo. +def min_zero_cos(a: u.physical.angle) -> u.Quantity: + """Use to ensure that cos(>=90 deg) equals 0.""" - """ + # handle scalars separately + if a.ndim == 0 and u.isclose(np.abs(a), 90 * u.deg): + return u.Quantity(0) - @dataclass_input - def __init__(self, phys: Phys): - self.phys = phys + x = np.cos(a) + x[u.isclose(np.abs(a), 90 * u.deg)] = 0 - @staticmethod - def _min_zero_cos(a: u.physical.angle) -> u.Quantity: - """Use to ensure that cos(>=90 deg) equals 0.""" - - # handle scalars separately - if a.ndim == 0 and u.isclose(np.abs(a), 90 * u.deg): - return u.Quantity(0) + return np.maximum(x, 0) - x = np.cos(a) - x[u.isclose(np.abs(a), 90 * u.deg)] = 0 - return np.maximum(x, 0) +class Surface(ABC): + """Abstract base class for all small-body surfaces.""" + @staticmethod @abstractmethod def absorption( - self, F_i: SpectralFluxDensityQuantity, + epsilon: float, i: u.physical.angle, ) -> u.Quantity: r"""Absorption of directional, incident light. @@ -52,9 +41,13 @@ def absorption( Parameters ---------- + F_i : `astropy.units.Quantity` Incident light (spectral flux density / spectral irradiance). + epsilon : float + Emissivity of the surface. + i : `~astropy.units.Quantity` Angle from normal of incident light. @@ -66,10 +59,11 @@ def absorption( """ + @staticmethod @abstractmethod def emission( - self, X_e: SpectralFluxDensityQuantity, + epsilon: float, e: u.physical.angle, phi: u.physical.angle, ) -> u.Quantity: @@ -84,6 +78,9 @@ def emission( X_e : `astropy.units.Quantity` Emitted spectral radiance. + epsilon : float + Emissivity of the surface. + e : `~astropy.units.Quantity` Observed angle from normal. @@ -99,9 +96,14 @@ def emission( """ + @staticmethod @abstractmethod def reflectance( - self, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle + F_i: SpectralFluxDensityQuantity, + albedo: float, + i: u.physical.angle, + e: u.physical.angle, + phi: u.physical.angle, ) -> u.Quantity: r"""Bidirectional reflectance. @@ -117,6 +119,9 @@ def reflectance( F_i : `astropy.units.Quantity` Incident light (spectral flux density). + albedo : float + Surface albedo. + i : `~astropy.units.Quantity` Angle from normal of incident light. @@ -180,8 +185,8 @@ def reflectance_from_vectors( Returns ------- - reflectance : `~astropy.units.Quantity` - Reflectance. + F_r : `~astropy.units.Quantity` + Spectral flux density / spectral irradiance received by the observer. """ diff --git a/sbpy/surfaces/tests/test_lambertian.py b/sbpy/surfaces/tests/test_lambertian.py index 685d0699..323bf476 100644 --- a/sbpy/surfaces/tests/test_lambertian.py +++ b/sbpy/surfaces/tests/test_lambertian.py @@ -9,35 +9,43 @@ from ..lambertian import LambertianSurface -class TestingSurface(LambertianSurface): - def radiance(): # pragma: no cover - pass - - -@pytest.fixture -def surface(): - return TestingSurface({"albedo": 0.1}) +def test_absorption(): + F_i = 1 * u.Jy + epsilon = 0.9 + i = Angle([0, 30, 45, 60, 90, 100], "deg") + expected = 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) * u.Jy + result = LambertianSurface.absorption(F_i, epsilon, i) + assert u.allclose(result, expected) -def test_emittance(surface: TestingSurface): +def test_emission(): + X_e = 1 * u.Jy + epsilon = 0.9 e = Angle([0, 30, 45, 60, 90, 100], "deg") - expected = [1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0] + expected = ( + 0.9 + * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) + / np.pi + * u.Jy + / u.sr + ) phi = np.random.rand(len(e)) * 360 * u.deg # independent of phi - result = surface.emittance(e, phi) - assert u.allclose(result, expected) - - -def test_absorptance(surface: TestingSurface): - i = Angle([0, 30, 45, 60, 90, 100], "deg") - expected = 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) - result = surface.absorptance(i) + result = LambertianSurface.emission(X_e, epsilon, e, phi) assert u.allclose(result, expected) -def test_reflectance(surface: TestingSurface): +def test_reflectance(): + F_i = 1 * u.Jy + albedo = 0.1 i = Angle([0, 30, 45, 60, 90, 0, 90, 100] * u.deg) e = Angle([0, 60, 45, 30, 0, 90, 90, 0] * u.deg) phi = np.random.rand(len(i)) * 360 * u.deg # independent of phi - result = surface.reflectance(i, e, phi) - expected = 0.1 * np.array([1, np.sqrt(3) / 4, 1 / 2, np.sqrt(3) / 4, 0, 0, 0, 0]) + result = LambertianSurface.reflectance(F_i, albedo, i, e, phi) + expected = ( + 0.1 + * np.array([1, np.sqrt(3) / 4, 1 / 2, np.sqrt(3) / 4, 0, 0, 0, 0]) + / np.pi + * u.Jy + / u.sr + ) assert u.allclose(result, expected) diff --git a/sbpy/surfaces/tests/test_scattered.py b/sbpy/surfaces/tests/test_scattered.py index f40e9f6b..7824b465 100644 --- a/sbpy/surfaces/tests/test_scattered.py +++ b/sbpy/surfaces/tests/test_scattered.py @@ -1,50 +1,50 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst +# # Licensed under a 3-clause BSD style license - see LICENSE.rst -import pytest -import numpy as np -from astropy import units as u +# import pytest +# import numpy as np +# from astropy import units as u -from ...calib import Sun, solar_spectrum -from ..scattered import LambertianSurfaceScatteredSunlight +# from ...calib import Sun, solar_spectrum +# from ..scattered import LambertianSurfaceScatteredSunlight -def test_scattered_light(): +# def test_scattered_light(): - surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) +# surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) - F_i = 1 * u.Unit("W/(m2 um)") - n = np.array([1, 0, 0]) - rs = [1, 1, 0] * u.au - ro = [1, -1, 0] * u.au +# F_i = 1 * u.Unit("W/(m2 um)") +# n = np.array([1, 0, 0]) +# rs = [1, 1, 0] * u.au +# ro = [1, -1, 0] * u.au - # albedo * F_i / rh**2 * cos(45)**2 - # 0.1 * 1 / 2 - expected = 0.05 * u.W / (u.m**2 * u.um * u.sr) - result = surface.radiance_from_vectors(F_i, n, rs, ro) - assert u.isclose(result, expected) +# # albedo * F_i / rh**2 * cos(45)**2 +# # 0.1 * 1 / 2 +# expected = 0.05 * u.W / (u.m**2 * u.um * u.sr) +# result = surface.radiance_from_vectors(F_i, n, rs, ro) +# assert u.isclose(result, expected) -def test_scattered_sunlight(): - pytest.importorskip("synphot") +# def test_scattered_sunlight(): +# pytest.importorskip("synphot") - surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) +# surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) - # fake an easy solar spectrum for testing - wave = [0.5, 0.55, 0.6] * u.um - spec = [0.5, 1.0, 1.5] * u.W / (u.m**2 * u.um) - with solar_spectrum.set(Sun.from_array(wave, spec)): - n = np.array([1, 0, 0]) - rs = [1, 1, 0] * u.au - ro = [1, -1, 0] * u.au +# # fake an easy solar spectrum for testing +# wave = [0.5, 0.55, 0.6] * u.um +# spec = [0.5, 1.0, 1.5] * u.W / (u.m**2 * u.um) +# with solar_spectrum.set(Sun.from_array(wave, spec)): +# n = np.array([1, 0, 0]) +# rs = [1, 1, 0] * u.au +# ro = [1, -1, 0] * u.au - # albedo * F_i / rh**2 * cos(45)**2 - # 0.1 * 1 / 2 / 2 - expected = 0.025 * u.W / (u.m**2 * u.um * u.sr) - result = surface.scattered_sunlight_from_vectors(0.55 * u.um, n, rs, ro) - assert u.isclose(result, expected) +# # albedo * F_i / rh**2 * cos(45)**2 +# # 0.1 * 1 / 2 / 2 +# expected = 0.025 * u.W / (u.m**2 * u.um * u.sr) +# result = surface.scattered_sunlight_from_vectors(0.55 * u.um, n, rs, ro) +# assert u.isclose(result, expected) - # again to test branching to Sun.observe - result = surface.scattered_sunlight_from_vectors( - (0.549, 0.55, 0.551) * u.um, n, rs, ro - ) - assert u.allclose(result[1], expected, rtol=0.01) +# # again to test branching to Sun.observe +# result = surface.scattered_sunlight_from_vectors( +# (0.549, 0.55, 0.551) * u.um, n, rs, ro +# ) +# assert u.allclose(result[1], expected, rtol=0.01) From aad7ae305f54f1d767514adddebcb115bad5a35d Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Sun, 14 Dec 2025 08:56:05 -0500 Subject: [PATCH 21/31] Check in before stripping flux density. --- docs/sbpy/surfaces.rst | 206 ++++++++++--------------- sbpy/surfaces/__init__.py | 5 - sbpy/surfaces/lambertian.py | 74 ++++++++- sbpy/surfaces/surface.py | 174 +++++++++++++++------ sbpy/surfaces/tests/test_lambertian.py | 40 +++-- sbpy/surfaces/tests/test_surface.py | 48 +----- sbpy/units/typing.py | 10 +- 7 files changed, 315 insertions(+), 242 deletions(-) diff --git a/docs/sbpy/surfaces.rst b/docs/sbpy/surfaces.rst index 78e390b6..72f6d123 100644 --- a/docs/sbpy/surfaces.rst +++ b/docs/sbpy/surfaces.rst @@ -1,114 +1,126 @@ Surfaces Module (`sbpy.surfaces`) ================================= +Introduction +------------ + .. admonition:: warning The surface module is being made available on a preview basis. The API is - subject to change. Feedback on the approach used is welcome. + subject to change. Feedback on the approach is welcome. -The ``surfaces`` module describes the interaction of electromagnetic radiation with surfaces. Sbpy uses the :math:`(i, e, \phi)` model (angle of incidence, angle of emittance, and phase angle) to describe how light scatters and emits light. It has a flexible system that can incorporate any surface scattering model that can be described with these three angles. A few built-in surface models are provided. +The ``surfaces`` module describes the interaction of electromagnetic radiation with surfaces. Sbpy uses the :math:`(i, e, \phi)` model (angle of incidence, angle of emittance, and phase angle) to describe how light scatters and emits light. It has a flexible system that can incorporate any surface scattering model that can be described with these three angles. .. figure:: ../static/scattering-vectors.svg :alt: Diagram of surface scattering and emission vectors Sbpy's geometric basis for surface scattering and emission: :math:`n` is the surface normal vector, :math:`r_s` is the radial vector to the light source, and :math:`r_o` is the radial vector to the observer. The angle of incidence (:math:`i`), angle of emittance (:math:`e`), phase angle (:math:`\phi`) are labeled. -A instance of a ``Surface`` will have methods to calculate absorptance, emittance, and reflectance. A radiance method is used to calculate the observed spectral radiance of a surface given incident light. - -Surfaces are expected to require albedo and/or emissivity. Conventions on which property is used and when should be defined by each class. For example, a surface that only calculates reflectance may only require albedo, but one that calculates thermal emission may use the convention of albedo for absorbed sunlight and emissivity for emitted thermal radiation. +An implementation of the ``Surface`` model has methods to calculate electromagnetic absorption, emission, and bidirectional reflectance. -Built-in surface models ------------------------ +Getting Started +--------------- -The model `~sbpy.surfaces.scattered.LambertianSurfaceScatteredSunlight` is used to observe sunlight scattered from a Lambertian surface (light scattered uniformly in all directions). +Currently a Lambertian surface model is implemented. A Lambertian surface scatters and emits light uniformly in all directions. -Create an instance of the ``LambertianSurfaceScatteredSunlight`` model, and calculate the absorptance, emittance, and reflectance for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)`:: +Create an instance of the ``LambertianSurface`` model, and calculate the absorption and bidirectional reflectance for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)` and an albedo of 0.1 (emissivity of 0.9), given incident spectral flux density of 100 W m:superscript:`-2` μm:superscript:`-1`:: >>> import astropy.units as u - >>> from sbpy.surfaces import LambertianSurfaceScatteredSunlight - >>> - >>> surface = LambertianSurfaceScatteredSunlight({"albedo": 0.1}) + >>> from sbpy.surfaces import LambertianSurface >>> + >>> surface = LambertianSurface() >>> i, e, phi = [30, 60, 90] * u.deg - >>> surface.absorptance(i) # doctest: +FLOAT_CMP - - >>> surface.emittance(e, phi) # doctest: +FLOAT_CMP - - >>> surface.reflectance(i, e, phi) # doctest: +FLOAT_CMP - - - -Radiance of scattered/emitted light -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -All `~sbpy.surfaces.surface.Surface` models have a :meth:`~sbpy.surfaces.surface.Surface.radiance` method that calculates the radiance of scattered/emitted light, given the spectral flux density of incident light at the surface:: - - >>> F_i = 1000 * u.W / u.m**2 / u.um # incident spectral flux density - >>> surface.radiance(F_i, i, e, phi) # doctest: +FLOAT_CMP - - + >>> albedo = 0.1 + >>> epsilon = 1 - albedo + >>> F_i = 100 * u.W / u.m**2 / u.um + >>> + >>> surface.absorption(F_i, epsilon, i=i) # doctest: +FLOAT_CMP + -Radiance of scattered sunlight -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Now calculate 10 μm emission from the same surface, assuming it is at a temperature of 200 K:: -``LambertianSurfaceScatteredSunlight`` is derived from the ``ScatteredSunlight`` class, which provides convenience methods for calculating the radiance of sunlight scattered off the surface:: - -.. doctest-requires:: synphot - - >>> wave = 0.55 * u.um - >>> rh = 1 * u.au # heliocentric distance of the surface - >>> surface.scattered_sunlight(wave, rh, i, e, phi) # doctest: +FLOAT_CMP - - -The solar spectrum is configured using Sbpy's calibration system. See :doc:`calib` for details. - - -Radiance from vectors -^^^^^^^^^^^^^^^^^^^^^ + >>> from sbpy.spectroscopy.sources import BlackbodySource + >>> + >>> bb = BlackbodySource(200 * u.K) + >>> I_e = bb(10 * u.um, unit="W/(m2 um)") + >>> I_e + + >>> surface.emission(I_e, epsilon, e=e, phi=phi) # doctest: +FLOAT_CMP + -As an alternative to using :math:`(i, e, \phi)`, radiance may be calculated using vectors that define the normal direction, radial vector of the light source, and radial vector of the observer:: +Calculate the bidirectional reflectance for :math:`e=0`, and a range of incident angles:: -.. doctest-requires:: synphot +.. plot:: + :context: reset + :nofigs: - >>> # the following vectors are equivalent to (i, e, phi) = (30, 60, 90) deg + >>> import numpy as np + >>> + >>> F_i = 100 * u.W / u.m**2 / u.um + >>> e = 0 * u.deg + >>> i = np.linspace(-90, 90) * u.deg + >>> phi = np.abs(e - i) # calculate phase angle + >>> r = surface.reflectance(F_i, albedo, i=i, e=e, phi=phi) + >>> r # doctest: +FLOAT_CMP + + +.. plot:: + :include-source: + :context: + + import matplotlib.pyplot as plt + fig, ax = plt.subplots() + ax.plot(i, r) + ax.set(xlabel="$i$ (deg)", ylabel="$F_λ$ (W/m$^2$ μm)") + + +Using Vectors Instead of Angles +------------------------------- + +As an alternative to using :math:`(i, e, \phi)`, results may be calculated using vectors that define the normal direction, radial vector of the light source, and radial vector of the observer:: + + >>> # vectors equivalent to (i, e, phi) = (30, 60, 90) deg >>> n = [1, 0, 0] - >>> rs = [0.866, 0.5, 0] * u.au - >>> ro = [0.5, -0.866, 0] * u.au + >>> r = [0.8660254, 0.5, 0] * u.au + >>> ro = [0.5, -0.8660254, 0] * u.au >>> - >>> surface.radiance_from_vectors(F_i, n, rs, ro) # doctest: +FLOAT_CMP - - >>> surface.scattered_sunlight_from_vectors(wave, n, rs, ro) # doctest: +FLOAT_CMP - - -Notice that heliocentric distance was not needed in the call to ``scattered_sunlight_from_vectors`` because it was already accounted for by the ``rs`` vector. + >>> surface.reflectance_from_vectors(F_i, epsilon, n=n, r=r, ro=ro) # doctest: +FLOAT_CMP + -Building your own surface models --------------------------------- +Build Your Own Surface Models +----------------------------- -Defining your own surface model is typically done by creating a new class based on `~sbpy.surfaces.surface.Surface`, and defining methods for ``absorptance``, ``emittance``, and ``reflectance``. The `~sbpy.surfaces.lambertian.Lambertian` model serves as a good example. For surface scattering problems, most users will combine their class with the `~sbpy.surfaces.scattered.ScatteredLight` or `~sbpy.surfaces.scattered.ScatteredSunlight` classes, which provide the ``radiance`` method to complete the ``Surface`` model. +To define your own surface model create a new class based on `~sbpy.surfaces.surface.Surface`, and define the methods for ``absorption``, ``emission``, and ``reflectance``. The `~sbpy.surfaces.lambertian.LambertianSurface` model serves as a reference. -Here, we define a new surface model with surface scattering proportional to :math:`\cos^2` based on the `~sbpy.surfaces.scattered.ScatteredSunlight` class:: +Here, we define a new surface model with surface properties proportional to :math:`\cos^2`:: -.. doctest-requires:: synphot - - >>> import numpy as np - >>> from sbpy.surfaces import Surface, ScatteredSunlight + >>> # use min_zero_cos(a) to ensure cos(a >= 90 deg) = 0 + >>> from sbpy.surfaces import Surface, min_zero_cos >>> - >>> class Cos2SurfaceScatteredSunlight(ScatteredSunlight): - ... """Absorption and scattering proportional to cos**2.""" + >>> class Cos2Surface(Surface): + ... """Absorption and emission proportional to :math:`\cos^2`.""" ... - ... def absorptance(self, i): - ... return (1 - self.phys["albedo"]) * np.cos(i)**2 + ... def absorption(self, F_i, epsilon, i): + ... return F_i * epsilon * min_zero_cos(i)**2 ... - ... def emittance(self, e, phi): - ... return np.cos(e)**2 + ... def emission(self, I_e, epsilon, e, phi): + ... return I_e * epsilon * min_zero_cos(e)**2 ... - ... def reflectance(self, i, e, phi): - ... return self.phys["albedo"] * np.cos(i)**2 * self.emittance(e, phi) + ... def reflectance(self, F_i, albedo, i, e, phi): + ... return * min_zero_cos(i)**2 * self.emission(e, phi) >>> - >>> surface = Cos2SurfaceScatteredSunlight({"albedo": 0.1}) + >>> surface = Cos2Surface() >>> surface.reflectance(i, e, phi) # doctest: +FLOAT_CMP >>> surface.radiance(F_i, i, e, phi) # doctest: +FLOAT_CMP @@ -116,54 +128,6 @@ Here, we define a new surface model with surface scattering proportional to :mat >>> surface.scattered_sunlight(wave, rh, i, e, phi) # doctest: +FLOAT_CMP -.. for test runs without synphot -.. testsetup:: - - >>> from sbpy.surfaces import Surface, ScatteredSunlight - -However, if a scattering model will be re-used with other classes, e.g., for scattered light and thermal emission modeling, then the most flexible approach is to base the model on ``Surface`` and have derived classes combine the model with scattering or thermal emission classes:: - - >>> class Cos2Surface(Surface): - ... """Abstract base class for absorption and scattering proportional to cos**2. - ... - ... We document this as an "abstract base class" because the ``radiance`` - ... method required by ``Surface`` is not implemented here, and therefore - ... this class cannot be used directly (i.e., it cannot be instantiated). - ... - ... """ - ... - ... def absorptance(self, i): - ... return (1 - self.phys["albedo"]) * np.cos(i)**2 - ... - ... def emittance(self, e, phi): - ... return np.cos(e)**2 - ... - ... def reflectance(self, i, e, phi): - ... return self.phys["albedo"] * np.cos(i)**2 * self.emittance(e, phi) - >>> - >>> class Cos2SurfaceScatteredSunlightAlt(Cos2Surface, ScatteredSunlight): - ... """Absorption and scattering proportional to cos**2. - ... - ... This class combines the ``absorptance``, ``emittance``, and ``reflectance`` - ... methods from ``Cos2Surface`` with the ``radiance`` and ``scattered_sunlight`` - ... methods of ``ScatteredSunlight``. - ... - ... Parameters - ... ... - ... """ - >>> - >>> class Cos2SurfaceThermalEmission(Cos2Surface, ThermalEmission): # doctest: +SKIP - ... """Absorption and thermal emission proportional to cos**2. - ... - ... This class combines the ``absorptance``, ``emittance``, and ``reflectance`` - ... from ``Cos2Surface`` with the ``radiance`` method of a fictitious - ... ``ThermalEmission`` class. - ... - ... Parameters - ... ... - ... """ - -Thanks to their base classes (``Cos2Surface``, ``ScatteredSunlight``, and ``ThermalEmission``), the ``Cos2SurfaceScatteredSunlightAlt`` and ``Cos2SurfaceThermalEmission`` classes are complete and no additional code is needed, just documentation. Reference/API ------------- diff --git a/sbpy/surfaces/__init__.py b/sbpy/surfaces/__init__.py index 2240d50f..cc82ec83 100644 --- a/sbpy/surfaces/__init__.py +++ b/sbpy/surfaces/__init__.py @@ -1,7 +1,2 @@ from .surface import Surface from .lambertian import LambertianSurface -from .scattered import ( - ScatteredLight, - ScatteredSunlight, - LambertianSurfaceScatteredSunlight, -) diff --git a/sbpy/surfaces/lambertian.py b/sbpy/surfaces/lambertian.py index 4f11d035..411b7e4e 100644 --- a/sbpy/surfaces/lambertian.py +++ b/sbpy/surfaces/lambertian.py @@ -1,21 +1,68 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst +from typing import Optional, Union import numpy as np import astropy.units as u from .surface import Surface, min_zero_cos -from ..units.typing import SpectralFluxDensityQuantity +from ..units.typing import SpectralFluxDensityQuantity, SpectralRadianceQuantity class LambertianSurface(Surface): - """Lambertian surface absorption, emission, and reflectance.""" + """Lambertian surface absorption, emission, and reflectance. + + The surface is illuminated at an angle of :math:`i`, and emitted toward an + angle of :math:`e`, both measured from the surface normal direction. + :math:`\phi` is the source-target-observer (phase) angle. Both the source + and the emitted light are assumed to be collimated. + + + Examples + -------- + + Absorption of light: + + >>> import astropy.units as u + >>> from sbpy.surfaces import LambertianSurface + >>> + >>> surface = LambertianSurface() + >>> i, e, phi = [30, 60, 90] * u.deg + >>> epsilon = 0.9 + >>> F_i = 100 * u.W / u.m**2 / u.um # incident light + >>> + >>> surface.absorption(F_i, epsilon, i=i) # doctest: +FLOAT_CMP + + + Thermal emission: + + >>> I_e = 100 * u.W / u.m**2 / u.um / u.sr + >>> surface.emission(I_e, epsilon, e=e, phi=phi) # doctest: +FLOAT_CMP + + + Bidirectional reflectance: + + >>> albedo = 1 - epsilon + >>> surface.reflectance(F_i, albedo, i=i, e=e, phi=phi) # doctest: +FLOAT_CMP + + + Using vector-based arguments: + + >>> n = [1, 0, 0] + >>> r = [0.8660254, 0.5, 0] * u.au + >>> ro = [0.5, -0.8660254, 0] * u.au + >>> surface.reflectance_from_vectors(F_i, epsilon, n=n, r=r, ro=ro) # doctest: +FLOAT_CMP + + + """ @staticmethod @u.quantity_input def absorption( F_i: SpectralFluxDensityQuantity, epsilon: float, + *, i: u.physical.angle, + **kwargs, ) -> u.Quantity: # use min_zero_cos(i) to ensure cos(>= 90 deg) = 0 cos_i = min_zero_cos(i) @@ -24,25 +71,40 @@ def absorption( @staticmethod @u.quantity_input def emission( - X_e: SpectralFluxDensityQuantity, + I_e: SpectralFluxDensityQuantity, epsilon: float, + *, e: u.physical.angle, - phi: u.physical.angle, + **kwargs, ) -> u.Quantity: # use min_zero_cos(e) to ensure cos(>= 90 deg) = 0 cos_e = min_zero_cos(e) - return X_e * epsilon * cos_e / np.pi / u.sr + return I_e * epsilon * cos_e @staticmethod @u.quantity_input def reflectance( F_i: SpectralFluxDensityQuantity, albedo: float, + *, i: u.physical.angle, e: u.physical.angle, - phi: u.physical.angle, + **kwargs, ) -> u.Quantity: # use min_zero_cos(theta) to ensure cos(>= 90 deg) = 0 cos_i = min_zero_cos(i) cos_e = min_zero_cos(e) return F_i * albedo * cos_i * cos_e / np.pi / u.sr + + @u.quantity_input + def emission_from_vectors( + self, + I_e: SpectralRadianceQuantity, + epsilon: float, + *, + n: np.ndarray, + ro: np.ndarray, + **kwargs, + ) -> SpectralRadianceQuantity: + _, e, _ = self._vectors_to_angles(n, n, ro) + return self.emission(I_e, epsilon, e=e) diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index 3fa74049..8ece5ef8 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -1,13 +1,13 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst from abc import ABC, abstractmethod +from typing import Optional import numpy as np from astropy import units as u -from ..data.phys import Phys from ..data.decorators import dataclass_input -from ..units.typing import SpectralFluxDensityQuantity +from ..units.typing import SpectralFluxDensityQuantity, SpectralRadianceQuantity def min_zero_cos(a: u.physical.angle) -> u.Quantity: @@ -29,70 +29,54 @@ class Surface(ABC): @staticmethod @abstractmethod def absorption( - F_i: SpectralFluxDensityQuantity, - epsilon: float, i: u.physical.angle, - ) -> u.Quantity: + ) -> SpectralFluxDensityQuantity: r"""Absorption of directional, incident light. - The surface is illuminated by incident flux density, :math:`F_i`, at an - angle of :math:`i`, measured from the surface normal direction. + The surface is illuminated at an angle of :math:`i`, measured from the + surface normal direction. Parameters ---------- - - F_i : `astropy.units.Quantity` - Incident light (spectral flux density / spectral irradiance). - - epsilon : float - Emissivity of the surface. - i : `~astropy.units.Quantity` Angle from normal of incident light. Returns ------- - F_a : `~astropy.units.Quantity` - Absorbed spectral flux density. + a : `~astropy.units.Quantity` + Fraction of incident light absorbed. """ @staticmethod @abstractmethod def emission( - X_e: SpectralFluxDensityQuantity, - epsilon: float, e: u.physical.angle, phi: u.physical.angle, - ) -> u.Quantity: - r"""Emission of light from a surface, as seen by a distant observer. + ) -> SpectralRadianceQuantity: + r"""Emission of directional light from a surface. The surface is observed at an angle of :math:`e`, measured from the - surface normal direction, and at a solar phase angle of :math:`phi`. + surface normal direction. Anisotropic emission is characterized by the + angle `phi`. Parameters ---------- - X_e : `astropy.units.Quantity` - Emitted spectral radiance. - - epsilon : float - Emissivity of the surface. e : `~astropy.units.Quantity` Observed angle from normal. phi : `~astropy.units.Quantity` - Source-target-observer (phase) angle. + Angle to account for anisotropic emission. Returns ------- F_e : `~astropy.units.Quantity` - Spectral flux density / spectral irradiance received by the - observer. + Spectral radiance / specific intensity received by the observer. """ @@ -101,17 +85,18 @@ def emission( def reflectance( F_i: SpectralFluxDensityQuantity, albedo: float, + *, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle, - ) -> u.Quantity: + ) -> SpectralRadianceQuantity: r"""Bidirectional reflectance. - The surface is illuminated by incident flux density (irradiance), - :math:`F_i`, at an angle of :math:`i`, and emitted toward an angle of - :math:`e`, measured from the surface normal direction. :math:`\phi` is - the source-target-observer (phase) angle. Both the source and the - emitted light are assumed to be collimated. + The surface is illuminated by incident spectral flux density + (irradiance), :math:`F_i`, at an angle of :math:`i`, and emitted toward + an angle of :math:`e`, measured from the surface normal direction. + :math:`\phi` is the source-target-observer (phase) angle. Both the + source and the emitted light are assumed to be collimated. Parameters @@ -134,36 +119,124 @@ def reflectance( Returns ------- - F_r : `~astropy.units.Quantity` - Spectral flux density / spectral irradiance received by the observer. + I_r : `~astropy.units.Quantity` + Spectral radiance / specific intensity received by the observer. """ @staticmethod def _vectors_to_angles( n: np.ndarray, - rs: u.physical.length, - ro: np.ndarray, - ) -> tuple: + r: u.physical.length, + ro: u.physical.length, + ) -> tuple[u.Quantity, u.Quantity, u.Quantity]: n_hat = n / np.linalg.norm(n) - rs_hat = rs / np.linalg.norm(rs) + r_hat = r / np.linalg.norm(r) ro_hat = ro / np.linalg.norm(ro) - i = u.Quantity(np.arccos(np.dot(n_hat, rs_hat)), "rad") + i = u.Quantity(np.arccos(np.dot(n_hat, r_hat)), "rad") e = u.Quantity(np.arccos(np.dot(n_hat, ro_hat)), "rad") - phi = u.Quantity(np.arccos(np.dot(rs_hat, ro_hat)), "rad") + phi = u.Quantity(np.arccos(np.dot(r_hat, ro_hat)), "rad") return i, e, phi + @u.quantity_input + def absorption_from_vectors( + self, + F_i: SpectralFluxDensityQuantity, + epsilon: float, + *, + n: np.ndarray, + r: u.physical.length, + ro: Optional[u.physical.length] = None, + ) -> SpectralFluxDensityQuantity: + """Vector-based alternative to `absorption`. + + Input vectors do not need to be normalized. + + + Parameters + ---------- + F_i : `astropy.units.Quantity` + Incident light (spectral flux density). + + epsilon : float + Emissivity of the surface. + + n : `numpy.ndarray` + Surface normal vector. + + r : `~astropy.units.Quantity` + Radial vector from the surface to the light source. + + ro : `~astropy.units.Quantity`, optional + Radial vector from the surface to the observer (ignored). + + + Returns + ------- + F_a : `~astropy.units.Quantity` + Absorbed spectral flux density. + + """ + + i, _, _ = self._vectors_to_angles(n, r, [1, 0, 0]) + return self.absorption(F_i, epsilon, i) + + @u.quantity_input + def emission_from_vectors( + self, + I_e: SpectralRadianceQuantity, + epsilon: float, + *, + n: np.ndarray, + ro: u.physical.length, + r: Optional[u.physical.length] = None, + ) -> SpectralRadianceQuantity: + r"""Vector-based alternative to `emission`. + + Input vectors do not need to be normalized. + + + Parameters + ---------- + I_e : `astropy.units.Quantity` + Emitted spectral radiance. + + epsilon : float + Emissivity of the surface. + + n : `numpy.ndarray` + Surface normal vector. + + ro : `numpy.ndarray` + Radial vector from the surface to the observer. + + r : `~astropy.units.Quantity`, optional + Radial vector from the surface to the light source (ignored). + + + Returns + ------- + F_e : `~astropy.units.Quantity` + Spectral radiance / specific intensity received by the observer. + + """ + + _, e, phi = self._vectors_to_angles(n, None, ro) + return self.emission(I_e, epsilon, e, phi) + @u.quantity_input def reflectance_from_vectors( self, F_i: SpectralFluxDensityQuantity, + albedo: float, + *, n: np.ndarray, - rs: u.physical.length, - ro: np.ndarray, - ) -> u.Quantity: - """Vector based alternative to reflectance(). + r: u.physical.length, + ro: u.physical.length, + ) -> SpectralRadianceQuantity: + """Vector-based alternative to `reflectance`. Input vectors do not need to be normalized. @@ -173,10 +246,13 @@ def reflectance_from_vectors( F_i : `astropy.units.Quantity` Incident light (spectral flux density). + albedo : float + Surface albedo. + n : `numpy.ndarray` Surface normal vector. - rs : `~astropy.units.Quantity` + r : `~astropy.units.Quantity` Radial vector from the surface to the light source. ro : `numpy.ndarray` @@ -190,4 +266,4 @@ def reflectance_from_vectors( """ - return self.reflectance(F_i, *self._vectors_to_angles(n, rs, ro)) + return self.reflectance(F_i, albedo, *self._vectors_to_angles(n, r, ro)) diff --git a/sbpy/surfaces/tests/test_lambertian.py b/sbpy/surfaces/tests/test_lambertian.py index 323bf476..64545ae8 100644 --- a/sbpy/surfaces/tests/test_lambertian.py +++ b/sbpy/surfaces/tests/test_lambertian.py @@ -1,7 +1,5 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst -import pytest - import numpy as np from astropy.coordinates import Angle from astropy import units as u @@ -14,23 +12,41 @@ def test_absorption(): epsilon = 0.9 i = Angle([0, 30, 45, 60, 90, 100], "deg") expected = 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) * u.Jy - result = LambertianSurface.absorption(F_i, epsilon, i) + + surface = LambertianSurface() + result = surface.absorption(F_i, epsilon, i=i) assert u.allclose(result, expected) def test_emission(): - X_e = 1 * u.Jy + I_e = 1 * u.Jy / u.sr epsilon = 0.9 e = Angle([0, 30, 45, 60, 90, 100], "deg") expected = ( - 0.9 - * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) - / np.pi - * u.Jy - / u.sr + 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) * u.Jy / u.sr ) + + surface = LambertianSurface() + result = surface.emission(I_e, epsilon, e=e) + assert u.allclose(result, expected) + phi = np.random.rand(len(e)) * 360 * u.deg # independent of phi - result = LambertianSurface.emission(X_e, epsilon, e, phi) + result = surface.emission(I_e, epsilon, e=e, phi=phi) + assert u.allclose(result, expected) + + +def test_emission_from_vectors(): + I_e = 1 * u.Jy / u.sr + epsilon = 0.9 + expected = (0.9 * 0.5) * u.Jy / u.sr + + # equivalent to (i, e, phi) = (30, 60, 90) deg + n = [1, 0, 0] + r = [0.8660254, 0.5, 0] * u.au + ro = [0.5, -0.8660254, 0] * u.au + + surface = LambertianSurface() + result = surface.emission_from_vectors(I_e, epsilon, n=n, r=r, ro=ro) assert u.allclose(result, expected) @@ -40,7 +56,9 @@ def test_reflectance(): i = Angle([0, 30, 45, 60, 90, 0, 90, 100] * u.deg) e = Angle([0, 60, 45, 30, 0, 90, 90, 0] * u.deg) phi = np.random.rand(len(i)) * 360 * u.deg # independent of phi - result = LambertianSurface.reflectance(F_i, albedo, i, e, phi) + + surface = LambertianSurface() + result = surface.reflectance(F_i, albedo, i=i, e=e, phi=phi) expected = ( 0.1 * np.array([1, np.sqrt(3) / 4, 1 / 2, np.sqrt(3) / 4, 0, 0, 0, 0]) diff --git a/sbpy/surfaces/tests/test_surface.py b/sbpy/surfaces/tests/test_surface.py index 3a9bff6e..82d97dcb 100644 --- a/sbpy/surfaces/tests/test_surface.py +++ b/sbpy/surfaces/tests/test_surface.py @@ -4,39 +4,18 @@ from astropy.coordinates import Angle from astropy import units as u -from ..surface import Surface - - -class TestingSurface(Surface): - def absorptance(self, i): # pragma: no cover - pass - - def emittance(self, e, phi): # pragma: no cover - pass - - def reflectance(self, i: Angle, e: Angle, phi: Angle) -> u.Quantity: - return np.cos(i) * np.cos(e) * np.sin(phi) - - def radiance(self, F_i, i, e, phi): - return F_i * self.reflectance(i, e, phi) / u.sr +from ..surface import Surface, min_zero_cos def test_min_zero_cos(): a = Angle([-91, -90, 0, 30, 45, 60, 90, 91], "deg") - result = Surface._min_zero_cos(a) + result = min_zero_cos(a) expected = [0, 0, 1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0] assert np.allclose(result, expected) # test scalars for i in range(len(a)): - assert np.isclose(Surface._min_zero_cos(a[i]), expected[i]) - - -def test_radiance(): - surface = TestingSurface({}) - F_i = 1664 * u.W / (u.m**2 * u.um) - result = surface.radiance(F_i, 30 * u.deg, 30 * u.deg, 60 * u.deg) - assert u.isclose(result, F_i * np.sqrt(27) / 8 / u.sr) + assert np.isclose(min_zero_cos(a[i]), expected[i]) def test_vectors_to_angles(): @@ -60,24 +39,3 @@ def test_vectors_to_angles(): i, e, phi = Surface._vectors_to_angles(n, rs, ro) assert u.isclose(e, 60 * u.deg) assert u.isclose(phi, 120 * u.deg) - - -def test_radiance_from_vectors(): - F_i = 1664 * u.W / (u.m**2 * u.um) - surface = TestingSurface({}) - - n = [1, 0, 0] - rs = [1, 1, 0] * u.au - ro = [1, -1, 0] * u.au - result = surface.radiance_from_vectors(F_i, n, rs, ro) - assert u.isclose(result, F_i / 2 / u.sr) - - n = [1, 0, 0] - rs = [1 / np.sqrt(3), 1, 0] * u.au - ro = [np.sqrt(3), -1, 0] * u.au - result = surface.radiance_from_vectors(F_i, n, rs, ro) - assert u.isclose(result, F_i * np.sqrt(3) / 4 / u.sr) - - ro = [1 / np.sqrt(3), -1, 0] * u.au - result = surface.radiance_from_vectors(F_i, n, rs, ro) - assert u.isclose(result, F_i * np.sqrt(3) / 8 / u.sr) diff --git a/sbpy/units/typing.py b/sbpy/units/typing.py index e06b3962..96f20c75 100644 --- a/sbpy/units/typing.py +++ b/sbpy/units/typing.py @@ -7,22 +7,22 @@ import astropy.units as u from astropy import __version__ as astropy_version -UnitLike = Union[str, u.Unit] -SpectralQuantity = Union[ +type UnitLike = Union[str, u.Unit] +type SpectralQuantity = Union[ u.Quantity[u.physical.length], u.Quantity[u.physical.frequency] ] -SpectralFluxDensityQuantity = Union[ +type SpectralFluxDensityQuantity = Union[ u.Quantity[u.physical.spectral_flux_density], u.Quantity[u.physical.spectral_flux_density_wav], ] if Version(astropy_version) < Version("6.1"): - SpectralRadianceQuantity = Union[ + type SpectralRadianceQuantity = Union[ u.Quantity[u.physical.spectral_flux_density / u.sr], u.Quantity[u.physical.spectral_flux_density_wav / u.sr], ] else: - SpectralRadianceQuantity = Union[ + type SpectralRadianceQuantity = Union[ u.Quantity[u.physical.surface_brightness], u.Quantity[u.physical.surface_brightness_wav], ] From bdbf7ed94489922529a93d8ed8e9db40530ddd40 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 15 Dec 2025 10:09:47 -0500 Subject: [PATCH 22/31] Checking in --- sbpy/surfaces/surface.py | 49 +++++++++------------------------------- sbpy/units/physical.py | 8 +++++++ 2 files changed, 19 insertions(+), 38 deletions(-) create mode 100644 sbpy/units/physical.py diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index 8ece5ef8..a182f585 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -8,6 +8,7 @@ from ..data.decorators import dataclass_input from ..units.typing import SpectralFluxDensityQuantity, SpectralRadianceQuantity +from .. import units as sbu def min_zero_cos(a: u.physical.angle) -> u.Quantity: @@ -30,7 +31,7 @@ class Surface(ABC): @abstractmethod def absorption( i: u.physical.angle, - ) -> SpectralFluxDensityQuantity: + ) -> u.physical.dimensionless: r"""Absorption of directional, incident light. The surface is illuminated at an angle of :math:`i`, measured from the @@ -55,7 +56,7 @@ def absorption( def emission( e: u.physical.angle, phi: u.physical.angle, - ) -> SpectralRadianceQuantity: + ) -> u.physical.dimensionless: r"""Emission of directional light from a surface. The surface is observed at an angle of :math:`e`, measured from the @@ -65,7 +66,6 @@ def emission( Parameters ---------- - e : `~astropy.units.Quantity` Observed angle from normal. @@ -75,38 +75,26 @@ def emission( Returns ------- - F_e : `~astropy.units.Quantity` - Spectral radiance / specific intensity received by the observer. + """ @staticmethod @abstractmethod def reflectance( - F_i: SpectralFluxDensityQuantity, - albedo: float, - *, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle, - ) -> SpectralRadianceQuantity: + ) -> sbu.physical.reflectance: r"""Bidirectional reflectance. - The surface is illuminated by incident spectral flux density - (irradiance), :math:`F_i`, at an angle of :math:`i`, and emitted toward - an angle of :math:`e`, measured from the surface normal direction. - :math:`\phi` is the source-target-observer (phase) angle. Both the - source and the emitted light are assumed to be collimated. + The surface is illuminated at an angle of :math:`i`, and observed at an + angle of :math:`e`, measured from the surface normal direction. + :math:`\phi` is the source-target-observer (phase) angle. Parameters ---------- - F_i : `astropy.units.Quantity` - Incident light (spectral flux density). - - albedo : float - Surface albedo. - i : `~astropy.units.Quantity` Angle from normal of incident light. @@ -119,8 +107,8 @@ def reflectance( Returns ------- - I_r : `~astropy.units.Quantity` - Spectral radiance / specific intensity received by the observer. + r : `~astropy.units.Quantity` + Bidirectional reflectance. """ @@ -142,13 +130,7 @@ def _vectors_to_angles( @u.quantity_input def absorption_from_vectors( - self, - F_i: SpectralFluxDensityQuantity, - epsilon: float, - *, - n: np.ndarray, - r: u.physical.length, - ro: Optional[u.physical.length] = None, + self, n: np.ndarray, r: u.physical.length ) -> SpectralFluxDensityQuantity: """Vector-based alternative to `absorption`. @@ -157,21 +139,12 @@ def absorption_from_vectors( Parameters ---------- - F_i : `astropy.units.Quantity` - Incident light (spectral flux density). - - epsilon : float - Emissivity of the surface. - n : `numpy.ndarray` Surface normal vector. r : `~astropy.units.Quantity` Radial vector from the surface to the light source. - ro : `~astropy.units.Quantity`, optional - Radial vector from the surface to the observer (ignored). - Returns ------- diff --git a/sbpy/units/physical.py b/sbpy/units/physical.py new file mode 100644 index 00000000..6370e0e0 --- /dev/null +++ b/sbpy/units/physical.py @@ -0,0 +1,8 @@ +""" +Physical quantity types. See `Astropy Physical Types`_ for more. + +.. _Astropy Physical Types: https://docs.astropy.org/en/latest/units/physical_types.html#physical-types + +""" + +import astropy.units as u From 3bf2b65f54b2fc1e0e6292dcd513307dbedec74b Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 16 Dec 2025 09:14:23 -0500 Subject: [PATCH 23/31] Strip incident and emitted flux densities. --- docs/sbpy/surfaces.rst | 83 ++++++++------------ sbpy/surfaces/lambertian.py | 74 +++++++----------- sbpy/surfaces/surface.py | 104 ++++++++++--------------- sbpy/surfaces/tests/test_lambertian.py | 66 +++++++++------- sbpy/surfaces/tests/test_surface.py | 38 ++++----- sbpy/units/physical.py | 2 + 6 files changed, 158 insertions(+), 209 deletions(-) diff --git a/docs/sbpy/surfaces.rst b/docs/sbpy/surfaces.rst index 72f6d123..6a4ed940 100644 --- a/docs/sbpy/surfaces.rst +++ b/docs/sbpy/surfaces.rst @@ -6,8 +6,7 @@ Introduction .. admonition:: warning - The surface module is being made available on a preview basis. The API is - subject to change. Feedback on the approach is welcome. + The surfaces module is being made available on a preview basis. The API is subject to change. Feedback on the approach is welcome. The ``surfaces`` module describes the interaction of electromagnetic radiation with surfaces. Sbpy uses the :math:`(i, e, \phi)` model (angle of incidence, angle of emittance, and phase angle) to describe how light scatters and emits light. It has a flexible system that can incorporate any surface scattering model that can be described with these three angles. @@ -22,32 +21,20 @@ An implementation of the ``Surface`` model has methods to calculate electromagne Getting Started --------------- -Currently a Lambertian surface model is implemented. A Lambertian surface scatters and emits light uniformly in all directions. +Currently the Lambertian surface model is implemented. A Lambertian surface absorbs and emits light uniformly in all directions. -Create an instance of the ``LambertianSurface`` model, and calculate the absorption and bidirectional reflectance for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)` and an albedo of 0.1 (emissivity of 0.9), given incident spectral flux density of 100 W m:superscript:`-2` μm:superscript:`-1`:: +Create an instance of the ``LambertianSurface`` model, and calculate the absorption and emission scale factors for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)`:: >>> import astropy.units as u >>> from sbpy.surfaces import LambertianSurface >>> >>> surface = LambertianSurface() >>> i, e, phi = [30, 60, 90] * u.deg - >>> albedo = 0.1 - >>> epsilon = 1 - albedo - >>> F_i = 100 * u.W / u.m**2 / u.um >>> - >>> surface.absorption(F_i, epsilon, i=i) # doctest: +FLOAT_CMP - - -Now calculate 10 μm emission from the same surface, assuming it is at a temperature of 200 K:: - - >>> from sbpy.spectroscopy.sources import BlackbodySource - >>> - >>> bb = BlackbodySource(200 * u.K) - >>> I_e = bb(10 * u.um, unit="W/(m2 um)") - >>> I_e - - >>> surface.emission(I_e, epsilon, e=e, phi=phi) # doctest: +FLOAT_CMP - + >>> surface.absorption(i) # doctest: +FLOAT_CMP + + >>> surface.emission(e, phi) # doctest: +FLOAT_CMP + Calculate the bidirectional reflectance for :math:`e=0`, and a range of incident angles:: @@ -57,22 +44,10 @@ Calculate the bidirectional reflectance for :math:`e=0`, and a range of incident >>> import numpy as np >>> - >>> F_i = 100 * u.W / u.m**2 / u.um >>> e = 0 * u.deg >>> i = np.linspace(-90, 90) * u.deg >>> phi = np.abs(e - i) # calculate phase angle - >>> r = surface.reflectance(F_i, albedo, i=i, e=e, phi=phi) - >>> r # doctest: +FLOAT_CMP - + >>> r = surface.reflectance(i, e, phi) .. plot:: :include-source: @@ -81,52 +56,56 @@ Calculate the bidirectional reflectance for :math:`e=0`, and a range of incident import matplotlib.pyplot as plt fig, ax = plt.subplots() ax.plot(i, r) - ax.set(xlabel="$i$ (deg)", ylabel="$F_λ$ (W/m$^2$ μm)") + ax.set(xlabel="$i$ (deg)", ylabel="$r$ (1 / sr)") Using Vectors Instead of Angles ------------------------------- -As an alternative to using :math:`(i, e, \phi)`, results may be calculated using vectors that define the normal direction, radial vector of the light source, and radial vector of the observer:: +As an alternative to using :math:`(i, e, \phi)`, results may be calculated using vectors that define the normal direction, radial vector to the light source, and radial vector to the observer:: >>> # vectors equivalent to (i, e, phi) = (30, 60, 90) deg >>> n = [1, 0, 0] >>> r = [0.8660254, 0.5, 0] * u.au >>> ro = [0.5, -0.8660254, 0] * u.au >>> - >>> surface.reflectance_from_vectors(F_i, epsilon, n=n, r=r, ro=ro) # doctest: +FLOAT_CMP - + >>> surface.reflectance_from_vectors(n, r, ro) # doctest: +FLOAT_CMP + Build Your Own Surface Models ----------------------------- -To define your own surface model create a new class based on `~sbpy.surfaces.surface.Surface`, and define the methods for ``absorption``, ``emission``, and ``reflectance``. The `~sbpy.surfaces.lambertian.LambertianSurface` model serves as a reference. +To define your own surface model create a new class based on `~sbpy.surfaces.surface.Surface`, and define the methods for ``absorption``, ``emission``, and ``reflectance``. The `~sbpy.surfaces.lambertian.LambertianSurface` model may serve as a reference. Here, we define a new surface model with surface properties proportional to :math:`\cos^2`:: >>> # use min_zero_cos(a) to ensure cos(a >= 90 deg) = 0 - >>> from sbpy.surfaces import Surface, min_zero_cos + >>> from sbpy.surfaces.surface import Surface, min_zero_cos >>> >>> class Cos2Surface(Surface): - ... """Absorption and emission proportional to :math:`\cos^2`.""" + ... """Absorption and emission proportional to :math:`\\cos^2`.""" ... - ... def absorption(self, F_i, epsilon, i): - ... return F_i * epsilon * min_zero_cos(i)**2 + ... def absorption(self, i): + ... return min_zero_cos(i)**2 ... - ... def emission(self, I_e, epsilon, e, phi): - ... return I_e * epsilon * min_zero_cos(e)**2 + ... def emission(self, e, phi): + ... return min_zero_cos(e)**2 ... - ... def reflectance(self, F_i, albedo, i, e, phi): - ... return * min_zero_cos(i)**2 * self.emission(e, phi) - >>> + ... def reflectance(self, i, e, phi): + ... return self.absorption(i) * self.emission(e, phi) / np.pi / u.sr + +Create and use an instance of our new model:: + >>> surface = Cos2Surface() + >>> i, e, phi = [30, 60, 90] * u.deg + >>> + >>> surface.absorption(i) # doctest: +FLOAT_CMP + + >>> surface.emission(e, phi) # doctest: +FLOAT_CMP + >>> surface.reflectance(i, e, phi) # doctest: +FLOAT_CMP - - >>> surface.radiance(F_i, i, e, phi) # doctest: +FLOAT_CMP - - >>> surface.scattered_sunlight(wave, rh, i, e, phi) # doctest: +FLOAT_CMP - + Reference/API diff --git a/sbpy/surfaces/lambertian.py b/sbpy/surfaces/lambertian.py index 411b7e4e..0b32c634 100644 --- a/sbpy/surfaces/lambertian.py +++ b/sbpy/surfaces/lambertian.py @@ -1,6 +1,6 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst -from typing import Optional, Union +from typing import Union import numpy as np import astropy.units as u @@ -13,7 +13,7 @@ class LambertianSurface(Surface): The surface is illuminated at an angle of :math:`i`, and emitted toward an angle of :math:`e`, both measured from the surface normal direction. - :math:`\phi` is the source-target-observer (phase) angle. Both the source + :math:`\\phi` is the source-target-observer (phase) angle. Both the source and the emitted light are assumed to be collimated. @@ -27,84 +27,62 @@ class LambertianSurface(Surface): >>> >>> surface = LambertianSurface() >>> i, e, phi = [30, 60, 90] * u.deg - >>> epsilon = 0.9 - >>> F_i = 100 * u.W / u.m**2 / u.um # incident light >>> - >>> surface.absorption(F_i, epsilon, i=i) # doctest: +FLOAT_CMP - + >>> surface.absorption(i) # doctest: +FLOAT_CMP + Thermal emission: - >>> I_e = 100 * u.W / u.m**2 / u.um / u.sr - >>> surface.emission(I_e, epsilon, e=e, phi=phi) # doctest: +FLOAT_CMP - + >>> surface.emission(e, phi) # doctest: +FLOAT_CMP + Bidirectional reflectance: - >>> albedo = 1 - epsilon - >>> surface.reflectance(F_i, albedo, i=i, e=e, phi=phi) # doctest: +FLOAT_CMP - + >>> surface.reflectance(i, e, phi) # doctest: +FLOAT_CMP + Using vector-based arguments: >>> n = [1, 0, 0] >>> r = [0.8660254, 0.5, 0] * u.au >>> ro = [0.5, -0.8660254, 0] * u.au - >>> surface.reflectance_from_vectors(F_i, epsilon, n=n, r=r, ro=ro) # doctest: +FLOAT_CMP - + >>> surface.reflectance_from_vectors(n, r, ro) # doctest: +FLOAT_CMP + """ - @staticmethod @u.quantity_input def absorption( - F_i: SpectralFluxDensityQuantity, - epsilon: float, - *, + self, i: u.physical.angle, - **kwargs, - ) -> u.Quantity: + ) -> u.Quantity[u.dimensionless_unscaled]: # use min_zero_cos(i) to ensure cos(>= 90 deg) = 0 - cos_i = min_zero_cos(i) - return F_i * epsilon * cos_i + return min_zero_cos(i) - @staticmethod @u.quantity_input def emission( - I_e: SpectralFluxDensityQuantity, - epsilon: float, - *, + self, e: u.physical.angle, - **kwargs, - ) -> u.Quantity: + phi: Union[u.physical.angle, None], + ) -> u.Quantity[u.dimensionless_unscaled]: # use min_zero_cos(e) to ensure cos(>= 90 deg) = 0 - cos_e = min_zero_cos(e) - return I_e * epsilon * cos_e + return min_zero_cos(e) - @staticmethod @u.quantity_input def reflectance( - F_i: SpectralFluxDensityQuantity, - albedo: float, - *, + self, i: u.physical.angle, e: u.physical.angle, - **kwargs, - ) -> u.Quantity: - # use min_zero_cos(theta) to ensure cos(>= 90 deg) = 0 - cos_i = min_zero_cos(i) - cos_e = min_zero_cos(e) - return F_i * albedo * cos_i * cos_e / np.pi / u.sr + phi: u.physical.angle, + ) -> u.Quantity[u.sr**-1]: + return self.absorption(i) * self.emission(e, phi) / np.pi / u.sr @u.quantity_input def emission_from_vectors( self, - I_e: SpectralRadianceQuantity, - epsilon: float, - *, n: np.ndarray, - ro: np.ndarray, - **kwargs, - ) -> SpectralRadianceQuantity: - _, e, _ = self._vectors_to_angles(n, n, ro) - return self.emission(I_e, epsilon, e=e) + r: Union[u.physical.length, None], + ro: u.physical.length, + ) -> u.Quantity[u.dimensionless_unscaled]: + e = self._angle(n, ro) + return self.emission(e, None) diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index a182f585..5deb2745 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -1,17 +1,18 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst from abc import ABC, abstractmethod -from typing import Optional +from typing import Union import numpy as np from astropy import units as u -from ..data.decorators import dataclass_input from ..units.typing import SpectralFluxDensityQuantity, SpectralRadianceQuantity -from .. import units as sbu +# for bidirectional_reflectance +from ..units import physical # noqa: F401 -def min_zero_cos(a: u.physical.angle) -> u.Quantity: + +def min_zero_cos(a: u.physical.angle) -> u.Quantity[u.dimensionless_unscaled]: """Use to ensure that cos(>=90 deg) equals 0.""" # handle scalars separately @@ -27,11 +28,8 @@ def min_zero_cos(a: u.physical.angle) -> u.Quantity: class Surface(ABC): """Abstract base class for all small-body surfaces.""" - @staticmethod @abstractmethod - def absorption( - i: u.physical.angle, - ) -> u.physical.dimensionless: + def absorption(self, i: u.physical.angle) -> u.Quantity[u.dimensionless_unscaled]: r"""Absorption of directional, incident light. The surface is illuminated at an angle of :math:`i`, measured from the @@ -51,12 +49,12 @@ def absorption( """ - @staticmethod @abstractmethod def emission( + self, e: u.physical.angle, phi: u.physical.angle, - ) -> u.physical.dimensionless: + ) -> u.Quantity[u.dimensionless_unscaled]: r"""Emission of directional light from a surface. The surface is observed at an angle of :math:`e`, measured from the @@ -75,17 +73,18 @@ def emission( Returns ------- - + em : `~astropy.units.Quantity` + Emission factor. """ - @staticmethod @abstractmethod def reflectance( + self, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle, - ) -> sbu.physical.reflectance: + ) -> u.physical.bidirectional_reflectance: r"""Bidirectional reflectance. The surface is illuminated at an angle of :math:`i`, and observed at an @@ -107,31 +106,24 @@ def reflectance( Returns ------- - r : `~astropy.units.Quantity` + bdr : `~astropy.units.Quantity` Bidirectional reflectance. """ @staticmethod - def _vectors_to_angles( - n: np.ndarray, - r: u.physical.length, - ro: u.physical.length, - ) -> tuple[u.Quantity, u.Quantity, u.Quantity]: - n_hat = n / np.linalg.norm(n) - r_hat = r / np.linalg.norm(r) - ro_hat = ro / np.linalg.norm(ro) - - i = u.Quantity(np.arccos(np.dot(n_hat, r_hat)), "rad") - e = u.Quantity(np.arccos(np.dot(n_hat, ro_hat)), "rad") - phi = u.Quantity(np.arccos(np.dot(r_hat, ro_hat)), "rad") - - return i, e, phi + def _angle(a: u.Quantity, b: u.Quantity) -> u.physical.angle: + a_hat = a / np.linalg.norm(a) + b_hat = b / np.linalg.norm(b) + return u.Quantity(np.arccos(np.dot(a_hat, b_hat)), "rad") @u.quantity_input def absorption_from_vectors( - self, n: np.ndarray, r: u.physical.length - ) -> SpectralFluxDensityQuantity: + self, + n: np.ndarray, + r: u.physical.length, + ro: Union[u.physical.length, None], + ) -> u.Quantity[u.dimensionless_unscaled]: """Vector-based alternative to `absorption`. Input vectors do not need to be normalized. @@ -145,6 +137,10 @@ def absorption_from_vectors( r : `~astropy.units.Quantity` Radial vector from the surface to the light source. + ro : `~astropy.units.Quantity` or ``None`` + Radial vector from the surface to the observer. Parameter is unused + and may be ``None``. + Returns ------- @@ -153,19 +149,16 @@ def absorption_from_vectors( """ - i, _, _ = self._vectors_to_angles(n, r, [1, 0, 0]) - return self.absorption(F_i, epsilon, i) + i = self._angle(n, r) + return self.absorption(i) @u.quantity_input def emission_from_vectors( self, - I_e: SpectralRadianceQuantity, - epsilon: float, - *, n: np.ndarray, + r: u.physical.length, ro: u.physical.length, - r: Optional[u.physical.length] = None, - ) -> SpectralRadianceQuantity: + ) -> u.Quantity[u.dimensionless_unscaled]: r"""Vector-based alternative to `emission`. Input vectors do not need to be normalized. @@ -173,21 +166,15 @@ def emission_from_vectors( Parameters ---------- - I_e : `astropy.units.Quantity` - Emitted spectral radiance. - - epsilon : float - Emissivity of the surface. - n : `numpy.ndarray` Surface normal vector. + r : `~astropy.units.Quantity` + Radial vector from the surface to the light source. + ro : `numpy.ndarray` Radial vector from the surface to the observer. - r : `~astropy.units.Quantity`, optional - Radial vector from the surface to the light source (ignored). - Returns ------- @@ -196,19 +183,17 @@ def emission_from_vectors( """ - _, e, phi = self._vectors_to_angles(n, None, ro) - return self.emission(I_e, epsilon, e, phi) + e = self._angle(n, ro) + phi = self._angle(r, ro) + return self.emission(e, phi) @u.quantity_input def reflectance_from_vectors( self, - F_i: SpectralFluxDensityQuantity, - albedo: float, - *, n: np.ndarray, r: u.physical.length, ro: u.physical.length, - ) -> SpectralRadianceQuantity: + ) -> u.Quantity[u.sr**-1]: """Vector-based alternative to `reflectance`. Input vectors do not need to be normalized. @@ -216,12 +201,6 @@ def reflectance_from_vectors( Parameters ---------- - F_i : `astropy.units.Quantity` - Incident light (spectral flux density). - - albedo : float - Surface albedo. - n : `numpy.ndarray` Surface normal vector. @@ -234,9 +213,12 @@ def reflectance_from_vectors( Returns ------- - F_r : `~astropy.units.Quantity` - Spectral flux density / spectral irradiance received by the observer. + bdr : `~astropy.units.Quantity` + Bidirectional reflectance. """ - return self.reflectance(F_i, albedo, *self._vectors_to_angles(n, r, ro)) + i = self._angle(n, r) + e = self._angle(n, ro) + phi = self._angle(r, ro) + return self.reflectance(i, e, phi) diff --git a/sbpy/surfaces/tests/test_lambertian.py b/sbpy/surfaces/tests/test_lambertian.py index 64545ae8..cbd723a7 100644 --- a/sbpy/surfaces/tests/test_lambertian.py +++ b/sbpy/surfaces/tests/test_lambertian.py @@ -8,62 +8,74 @@ def test_absorption(): - F_i = 1 * u.Jy - epsilon = 0.9 i = Angle([0, 30, 45, 60, 90, 100], "deg") - expected = 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) * u.Jy + expected = np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) surface = LambertianSurface() - result = surface.absorption(F_i, epsilon, i=i) + result = surface.absorption(i) assert u.allclose(result, expected) -def test_emission(): - I_e = 1 * u.Jy / u.sr - epsilon = 0.9 - e = Angle([0, 30, 45, 60, 90, 100], "deg") - expected = ( - 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) * u.Jy / u.sr - ) +def test_absorption_from_vectors(): + # equivalent to (i, e, phi) = (30, 60, 90) deg + n = [1, 0, 0] + r = [0.8660254, 0.5, 0] * u.au + ro = [0.5, -0.8660254, 0] * u.au + + expected = np.sqrt(3) / 2 surface = LambertianSurface() - result = surface.emission(I_e, epsilon, e=e) + result = surface.absorption_from_vectors(n, r, ro) assert u.allclose(result, expected) - phi = np.random.rand(len(e)) * 360 * u.deg # independent of phi - result = surface.emission(I_e, epsilon, e=e, phi=phi) + +def test_emission(): + e = Angle([0, 30, 45, 60, 90, 100], "deg") + expected = np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) + + surface = LambertianSurface() + result = surface.emission(e, None) assert u.allclose(result, expected) def test_emission_from_vectors(): - I_e = 1 * u.Jy / u.sr - epsilon = 0.9 - expected = (0.9 * 0.5) * u.Jy / u.sr - # equivalent to (i, e, phi) = (30, 60, 90) deg n = [1, 0, 0] r = [0.8660254, 0.5, 0] * u.au ro = [0.5, -0.8660254, 0] * u.au + expected = 0.5 + surface = LambertianSurface() - result = surface.emission_from_vectors(I_e, epsilon, n=n, r=r, ro=ro) + result = surface.emission_from_vectors(n, r, ro) + assert u.allclose(result, expected) + + # incident vector is optional + result = surface.emission_from_vectors(n, None, ro) assert u.allclose(result, expected) def test_reflectance(): - F_i = 1 * u.Jy - albedo = 0.1 i = Angle([0, 30, 45, 60, 90, 0, 90, 100] * u.deg) e = Angle([0, 60, 45, 30, 0, 90, 90, 0] * u.deg) phi = np.random.rand(len(i)) * 360 * u.deg # independent of phi surface = LambertianSurface() - result = surface.reflectance(F_i, albedo, i=i, e=e, phi=phi) + result = surface.reflectance(i, e, phi) expected = ( - 0.1 - * np.array([1, np.sqrt(3) / 4, 1 / 2, np.sqrt(3) / 4, 0, 0, 0, 0]) - / np.pi - * u.Jy - / u.sr + np.array([1, np.sqrt(3) / 4, 1 / 2, np.sqrt(3) / 4, 0, 0, 0, 0]) / np.pi / u.sr ) assert u.allclose(result, expected) + + +def test_reflectance_from_vectors(): + # equivalent to (i, e, phi) = (30, 60, 90) deg + n = [1, 0, 0] + r = [0.8660254, 0.5, 0] * u.au + ro = [0.5, -0.8660254, 0] * u.au + + expected = np.sqrt(3) / 4 / np.pi / u.sr + + surface = LambertianSurface() + result = surface.reflectance_from_vectors(n, r, ro) + assert u.allclose(result, expected) diff --git a/sbpy/surfaces/tests/test_surface.py b/sbpy/surfaces/tests/test_surface.py index 82d97dcb..2b50918d 100644 --- a/sbpy/surfaces/tests/test_surface.py +++ b/sbpy/surfaces/tests/test_surface.py @@ -1,5 +1,7 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst +import pytest + import numpy as np from astropy.coordinates import Angle from astropy import units as u @@ -18,24 +20,18 @@ def test_min_zero_cos(): assert np.isclose(min_zero_cos(a[i]), expected[i]) -def test_vectors_to_angles(): - n = [1, 0, 0] - rs = [1, 1, 0] * u.au - ro = [1, -1, 0] * u.au - i, e, phi = Surface._vectors_to_angles(n, rs, ro) - assert u.isclose(i, 45 * u.deg) - assert u.isclose(e, 45 * u.deg) - assert u.isclose(phi, 90 * u.deg) - - n = [1, 0, 0] - rs = [1 / np.sqrt(3), 1, 0] * u.au - ro = [np.sqrt(3), -1, 0] * u.au - i, e, phi = Surface._vectors_to_angles(n, rs, ro) - assert u.isclose(i, 60 * u.deg) - assert u.isclose(e, 30 * u.deg) - assert u.isclose(phi, 90 * u.deg) - - ro = [1 / np.sqrt(3), -1, 0] * u.au - i, e, phi = Surface._vectors_to_angles(n, rs, ro) - assert u.isclose(e, 60 * u.deg) - assert u.isclose(phi, 120 * u.deg) +@pytest.mark.parametrize( + "a,b,theta", + [ + ([1, 0, 0], [1, 1, 0], 45), + ([1, 0, 0], [1, -1, 0], 45), + ([1, 1, 0], [1, -1, 0], 90), + ([1, 0, 0], [1 / np.sqrt(3), 1, 0], 60), + ([1, 0, 0], [np.sqrt(3), -1, 0], 30), + ([1 / np.sqrt(3), 1, 0], [np.sqrt(3), -1, 0], 90), + ([1, 0, 0], [1 / np.sqrt(3), -1, 0], 60), + ([1 / np.sqrt(3), 1, 0], [1 / np.sqrt(3), -1, 0], 120), + ], +) +def test_angle(a, b, theta): + assert u.isclose(Surface._angle(a, b * u.au), theta * u.deg) diff --git a/sbpy/units/physical.py b/sbpy/units/physical.py index 6370e0e0..bdfe4da1 100644 --- a/sbpy/units/physical.py +++ b/sbpy/units/physical.py @@ -6,3 +6,5 @@ """ import astropy.units as u + +u.physical.def_physical_type(u.sr**-1, "bidirectional reflectance") From 688f253b00c45c2f43d54914d0d0883728a115e4 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 16 Dec 2025 10:49:22 -0500 Subject: [PATCH 24/31] Add albedo and emissivity back in. --- docs/sbpy/surfaces.rst | 46 +++++++------ sbpy/surfaces/__init__.py | 1 + sbpy/surfaces/lambertian.py | 35 ++++++---- sbpy/surfaces/scattered.py | 90 ++++++++++++++++++++++++-- sbpy/surfaces/surface.py | 33 ++++++++-- sbpy/surfaces/tests/test_lambertian.py | 53 ++++++++++----- 6 files changed, 200 insertions(+), 58 deletions(-) diff --git a/docs/sbpy/surfaces.rst b/docs/sbpy/surfaces.rst index 6a4ed940..67ab8dac 100644 --- a/docs/sbpy/surfaces.rst +++ b/docs/sbpy/surfaces.rst @@ -23,18 +23,20 @@ Getting Started Currently the Lambertian surface model is implemented. A Lambertian surface absorbs and emits light uniformly in all directions. -Create an instance of the ``LambertianSurface`` model, and calculate the absorption and emission scale factors for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)`:: +Create an instance of the ``LambertianSurface`` model, and calculate the absorption and emission scale factors for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)`. Let the albedo be 0.1 (emissivity = 0.9):: >>> import astropy.units as u >>> from sbpy.surfaces import LambertianSurface >>> - >>> surface = LambertianSurface() + >>> albedo = 0.1 + >>> epsilon = 1 - albedo >>> i, e, phi = [30, 60, 90] * u.deg >>> - >>> surface.absorption(i) # doctest: +FLOAT_CMP - - >>> surface.emission(e, phi) # doctest: +FLOAT_CMP - + >>> surface = LambertianSurface() + >>> surface.absorption(epsilon, i) # doctest: +FLOAT_CMP + + >>> surface.emission(epsilon, e, phi) # doctest: +FLOAT_CMP + Calculate the bidirectional reflectance for :math:`e=0`, and a range of incident angles:: @@ -47,7 +49,7 @@ Calculate the bidirectional reflectance for :math:`e=0`, and a range of incident >>> e = 0 * u.deg >>> i = np.linspace(-90, 90) * u.deg >>> phi = np.abs(e - i) # calculate phase angle - >>> r = surface.reflectance(i, e, phi) + >>> r = surface.reflectance(albedo, i, e, phi) .. plot:: :include-source: @@ -69,8 +71,8 @@ As an alternative to using :math:`(i, e, \phi)`, results may be calculated using >>> r = [0.8660254, 0.5, 0] * u.au >>> ro = [0.5, -0.8660254, 0] * u.au >>> - >>> surface.reflectance_from_vectors(n, r, ro) # doctest: +FLOAT_CMP - + >>> surface.reflectance_from_vectors(albedo, n, r, ro) # doctest: +FLOAT_CMP + Build Your Own Surface Models @@ -86,26 +88,28 @@ Here, we define a new surface model with surface properties proportional to :mat >>> class Cos2Surface(Surface): ... """Absorption and emission proportional to :math:`\\cos^2`.""" ... - ... def absorption(self, i): - ... return min_zero_cos(i)**2 + ... def absorption(self, epsilon, i): + ... return epsilon * min_zero_cos(i)**2 ... - ... def emission(self, e, phi): - ... return min_zero_cos(e)**2 + ... def emission(self, epsilon, e, phi): + ... return epsilon * min_zero_cos(e)**2 ... - ... def reflectance(self, i, e, phi): - ... return self.absorption(i) * self.emission(e, phi) / np.pi / u.sr + ... def reflectance(self, albedo, i, e, phi): + ... return albedo * min_zero_cos(i)**2 * min_zero_cos(e)**2 / np.pi / u.sr Create and use an instance of our new model:: >>> surface = Cos2Surface() + >>> albedo = 0.1 + >>> epsilon = 1 - albedo >>> i, e, phi = [30, 60, 90] * u.deg >>> - >>> surface.absorption(i) # doctest: +FLOAT_CMP - - >>> surface.emission(e, phi) # doctest: +FLOAT_CMP - - >>> surface.reflectance(i, e, phi) # doctest: +FLOAT_CMP - + >>> surface.absorption(epsilon, i) # doctest: +FLOAT_CMP + + >>> surface.emission(epsilon, e, phi) # doctest: +FLOAT_CMP + + >>> surface.reflectance(albedo, i, e, phi) # doctest: +FLOAT_CMP + Reference/API diff --git a/sbpy/surfaces/__init__.py b/sbpy/surfaces/__init__.py index cc82ec83..6c360370 100644 --- a/sbpy/surfaces/__init__.py +++ b/sbpy/surfaces/__init__.py @@ -1,2 +1,3 @@ from .surface import Surface from .lambertian import LambertianSurface +from .scattered import ScatteredLight diff --git a/sbpy/surfaces/lambertian.py b/sbpy/surfaces/lambertian.py index 0b32c634..d45f693b 100644 --- a/sbpy/surfaces/lambertian.py +++ b/sbpy/surfaces/lambertian.py @@ -20,69 +20,76 @@ class LambertianSurface(Surface): Examples -------- - Absorption of light: + Absorption of light for a surface with an albedo of 0.1 and an emissivity of + 0.9: >>> import astropy.units as u >>> from sbpy.surfaces import LambertianSurface >>> >>> surface = LambertianSurface() + >>> albedo = 0.1 + >>> epsilon = 1 - albedo >>> i, e, phi = [30, 60, 90] * u.deg >>> - >>> surface.absorption(i) # doctest: +FLOAT_CMP - + >>> surface.absorption(epsilon, i) # doctest: +FLOAT_CMP + - Thermal emission: + Emission: - >>> surface.emission(e, phi) # doctest: +FLOAT_CMP - + >>> surface.emission(epsilon, e, phi) # doctest: +FLOAT_CMP + Bidirectional reflectance: - >>> surface.reflectance(i, e, phi) # doctest: +FLOAT_CMP - + >>> surface.reflectance(albedo, i, e, phi) # doctest: +FLOAT_CMP + Using vector-based arguments: >>> n = [1, 0, 0] >>> r = [0.8660254, 0.5, 0] * u.au >>> ro = [0.5, -0.8660254, 0] * u.au - >>> surface.reflectance_from_vectors(n, r, ro) # doctest: +FLOAT_CMP - + >>> surface.reflectance_from_vectors(albedo, n, r, ro) # doctest: +FLOAT_CMP + """ @u.quantity_input def absorption( self, + epsilon: float, i: u.physical.angle, ) -> u.Quantity[u.dimensionless_unscaled]: # use min_zero_cos(i) to ensure cos(>= 90 deg) = 0 - return min_zero_cos(i) + return epsilon * min_zero_cos(i) @u.quantity_input def emission( self, + epsilon: float, e: u.physical.angle, phi: Union[u.physical.angle, None], ) -> u.Quantity[u.dimensionless_unscaled]: # use min_zero_cos(e) to ensure cos(>= 90 deg) = 0 - return min_zero_cos(e) + return epsilon * min_zero_cos(e) @u.quantity_input def reflectance( self, + albedo: float, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle, ) -> u.Quantity[u.sr**-1]: - return self.absorption(i) * self.emission(e, phi) / np.pi / u.sr + return albedo * min_zero_cos(i) * min_zero_cos(e) / np.pi / u.sr @u.quantity_input def emission_from_vectors( self, + epsilon: float, n: np.ndarray, r: Union[u.physical.length, None], ro: u.physical.length, ) -> u.Quantity[u.dimensionless_unscaled]: e = self._angle(n, ro) - return self.emission(e, None) + return self.emission(epsilon, e, None) diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py index 14509d9d..afa1d3ed 100644 --- a/sbpy/surfaces/scattered.py +++ b/sbpy/surfaces/scattered.py @@ -1,15 +1,97 @@ -# # Licensed under a 3-clause BSD style license - see LICENSE.rst +# Licensed under a 3-clause BSD style license - see LICENSE.rst # from abc import ABC, abstractmethod # import numpy as np -# import astropy.units as u +from typing import Optional +import astropy.units as u + +from ..spectroscopy.sources import SpectralSource, SinglePointSpectrumError +from .surface import Surface -# from .surface import Surface -# from .lambertian import LambertianSurface # from ..calib import Sun # from ..units.typing import SpectralQuantity, SpectralFluxDensityQuantity, UnitLike +from ..units.typing import SpectralRadianceQuantity, UnitLike + + +class ScatteredLight: + """Scatter light off a surface. + + + Examples + -------- + + Scatter sunlight from a Lambertian surface. + + >>> import astropy.units as u + >>> from sbpy.surfaces import ScatteredLight, LambertianSurface + >>> from sbpy.calib import Sun + >>> + >>> sun = Sun.from_default() + >>> surface = LambertianSurface() + >>> scattered = ScatteredLight(surface, sun) + >>> + >>> wave = 550 * u.nm + >>> albedo = 0.1 + >>> i, e, phi = [30, 60, 90] * u.deg + >>> scattered.radiance(wave, albedo, i, e, phi) # doctest: +FLOAT_CMP + + + """ + + def __init__(self, surface: Surface, source: SpectralSource): + self.surface = surface + self.source = source + + @u.quantity_input + def radiance( + self, + wfb, + albedo: u.physical.dimensionless, + i: u.physical.angle, + e: u.physical.angle, + phi: u.physical.angle, + unit: Optional[UnitLike] = None, + ) -> SpectralRadianceQuantity: + """Spectral radiance scattered from the surface. + + + Parameters + ---------- + wfb : `~astropy.units.Quantity`, `~synphot.SpectralElement`, list + Wavelengths, frequencies, bandpass, or list of + bandpasses of the observation. Bandpasses require + `~synphot`. + + albedo : `~astropy.units.Quantity` + Surface albedo for given `wfb` + + i : `~astropy.units.Quantity` + Angle from normal of incident light. + + e : `~astropy.units.Quantity` + Angle from normal of emitted light. + + phi : `~astropy.units.Quantity` + Source-target-observer (phase) angle. + + unit : str or `~astropy.units.Unit`, optional + Spectral radiance units of the result. + + + Returns + ------- + I : `~astropy.units.Quantity` + + """ + + try: + F_i = self.source.observe(wfb, unit=unit) + except SinglePointSpectrumError: + F_i = self.source(wfb, unit=unit) + + return F_i * self.surface.reflectance(albedo, i, e, phi) # class ScatteredLight(ABC): diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index 5deb2745..6f640319 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -29,7 +29,9 @@ class Surface(ABC): """Abstract base class for all small-body surfaces.""" @abstractmethod - def absorption(self, i: u.physical.angle) -> u.Quantity[u.dimensionless_unscaled]: + def absorption( + self, epsilon: u.physical.dimensionless, i: u.physical.angle + ) -> u.Quantity[u.dimensionless_unscaled]: r"""Absorption of directional, incident light. The surface is illuminated at an angle of :math:`i`, measured from the @@ -38,6 +40,9 @@ def absorption(self, i: u.physical.angle) -> u.Quantity[u.dimensionless_unscaled Parameters ---------- + epsilon : `~astropy.units.Quantity` + Surface emissivity. + i : `~astropy.units.Quantity` Angle from normal of incident light. @@ -52,6 +57,7 @@ def absorption(self, i: u.physical.angle) -> u.Quantity[u.dimensionless_unscaled @abstractmethod def emission( self, + epsilon: u.physical.dimensionless, e: u.physical.angle, phi: u.physical.angle, ) -> u.Quantity[u.dimensionless_unscaled]: @@ -64,6 +70,9 @@ def emission( Parameters ---------- + epsilon : `~astropy.units.Quantity` + Surface emissivity. + e : `~astropy.units.Quantity` Observed angle from normal. @@ -81,6 +90,7 @@ def emission( @abstractmethod def reflectance( self, + albedo: u.physical.dimensionless, i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle, @@ -94,6 +104,9 @@ def reflectance( Parameters ---------- + albedo : `~astropy.units.Quantity` + Surface albedo. + i : `~astropy.units.Quantity` Angle from normal of incident light. @@ -120,6 +133,7 @@ def _angle(a: u.Quantity, b: u.Quantity) -> u.physical.angle: @u.quantity_input def absorption_from_vectors( self, + epsilon: u.physical.dimensionless, n: np.ndarray, r: u.physical.length, ro: Union[u.physical.length, None], @@ -131,6 +145,9 @@ def absorption_from_vectors( Parameters ---------- + epsilon : `~astropy.units.Quantity` + Surface emissivity. + n : `numpy.ndarray` Surface normal vector. @@ -150,11 +167,12 @@ def absorption_from_vectors( """ i = self._angle(n, r) - return self.absorption(i) + return self.absorption(epsilon, i) @u.quantity_input def emission_from_vectors( self, + epsilon: u.physical.dimensionless, n: np.ndarray, r: u.physical.length, ro: u.physical.length, @@ -166,6 +184,9 @@ def emission_from_vectors( Parameters ---------- + epsilon : `~astropy.units.Quantity` + Surface emissivity. + n : `numpy.ndarray` Surface normal vector. @@ -185,11 +206,12 @@ def emission_from_vectors( e = self._angle(n, ro) phi = self._angle(r, ro) - return self.emission(e, phi) + return self.emission(epsilon, e, phi) @u.quantity_input def reflectance_from_vectors( self, + albedo: u.physical.dimensionless, n: np.ndarray, r: u.physical.length, ro: u.physical.length, @@ -201,6 +223,9 @@ def reflectance_from_vectors( Parameters ---------- + albedo : `~astropy.units.Quantity` + Surface albedo. + n : `numpy.ndarray` Surface normal vector. @@ -221,4 +246,4 @@ def reflectance_from_vectors( i = self._angle(n, r) e = self._angle(n, ro) phi = self._angle(r, ro) - return self.reflectance(i, e, phi) + return self.reflectance(albedo, i, e, phi) diff --git a/sbpy/surfaces/tests/test_lambertian.py b/sbpy/surfaces/tests/test_lambertian.py index cbd723a7..fcaff0ba 100644 --- a/sbpy/surfaces/tests/test_lambertian.py +++ b/sbpy/surfaces/tests/test_lambertian.py @@ -8,74 +8,97 @@ def test_absorption(): + epsilon = 0.9 i = Angle([0, 30, 45, 60, 90, 100], "deg") - expected = np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) + expected = 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) surface = LambertianSurface() - result = surface.absorption(i) + result = surface.absorption(epsilon, i) + assert u.allclose(result, expected) def test_absorption_from_vectors(): + epsilon = 0.9 # equivalent to (i, e, phi) = (30, 60, 90) deg n = [1, 0, 0] r = [0.8660254, 0.5, 0] * u.au ro = [0.5, -0.8660254, 0] * u.au - expected = np.sqrt(3) / 2 + expected = 0.9 * np.sqrt(3) / 2 surface = LambertianSurface() - result = surface.absorption_from_vectors(n, r, ro) + result = surface.absorption_from_vectors(epsilon, n, r, ro) + assert u.allclose(result, expected) def test_emission(): + epsilon = 0.9 e = Angle([0, 30, 45, 60, 90, 100], "deg") - expected = np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) + expected = 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) surface = LambertianSurface() - result = surface.emission(e, None) + result = surface.emission(epsilon, e, None) assert u.allclose(result, expected) def test_emission_from_vectors(): + epsilon = 0.9 # equivalent to (i, e, phi) = (30, 60, 90) deg n = [1, 0, 0] r = [0.8660254, 0.5, 0] * u.au ro = [0.5, -0.8660254, 0] * u.au - expected = 0.5 + expected = 0.45 surface = LambertianSurface() - result = surface.emission_from_vectors(n, r, ro) + result = surface.emission_from_vectors(epsilon, n, r, ro) assert u.allclose(result, expected) # incident vector is optional - result = surface.emission_from_vectors(n, None, ro) + result = surface.emission_from_vectors(epsilon, n, None, ro) assert u.allclose(result, expected) def test_reflectance(): + albedo = 0.1 i = Angle([0, 30, 45, 60, 90, 0, 90, 100] * u.deg) e = Angle([0, 60, 45, 30, 0, 90, 90, 0] * u.deg) phi = np.random.rand(len(i)) * 360 * u.deg # independent of phi - - surface = LambertianSurface() - result = surface.reflectance(i, e, phi) expected = ( - np.array([1, np.sqrt(3) / 4, 1 / 2, np.sqrt(3) / 4, 0, 0, 0, 0]) / np.pi / u.sr + np.array( + [ + 0.1, + (0.1 * np.sqrt(3) / 2) / 2, + (0.1 * np.sqrt(2) / 2) * np.sqrt(2) / 2, + (0.1 / 2) * np.sqrt(3) / 2, + 0, + 0, + 0, + 0, + ] + ) + / np.pi + / u.sr ) + + surface = LambertianSurface() + result = surface.reflectance(albedo, i, e, phi) + assert u.allclose(result, expected) def test_reflectance_from_vectors(): + albedo = 0.1 # equivalent to (i, e, phi) = (30, 60, 90) deg n = [1, 0, 0] r = [0.8660254, 0.5, 0] * u.au ro = [0.5, -0.8660254, 0] * u.au - expected = np.sqrt(3) / 4 / np.pi / u.sr + expected = (0.1 * np.sqrt(3) / 2) / 2 / np.pi / u.sr surface = LambertianSurface() - result = surface.reflectance_from_vectors(n, r, ro) + result = surface.reflectance_from_vectors(albedo, n, r, ro) + assert u.allclose(result, expected) From cdf0941d5e705c02c466d7c08508406c746525dc Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 16 Dec 2025 11:03:22 -0500 Subject: [PATCH 25/31] Clean up --- sbpy/surfaces/surface.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index 6f640319..c999dd6e 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -6,11 +6,6 @@ import numpy as np from astropy import units as u -from ..units.typing import SpectralFluxDensityQuantity, SpectralRadianceQuantity - -# for bidirectional_reflectance -from ..units import physical # noqa: F401 - def min_zero_cos(a: u.physical.angle) -> u.Quantity[u.dimensionless_unscaled]: """Use to ensure that cos(>=90 deg) equals 0.""" @@ -94,7 +89,7 @@ def reflectance( i: u.physical.angle, e: u.physical.angle, phi: u.physical.angle, - ) -> u.physical.bidirectional_reflectance: + ) -> u.Quantity[u.sr**-1]: r"""Bidirectional reflectance. The surface is illuminated at an angle of :math:`i`, and observed at an From f63316d30cf2ff2cd6027dd90cf24021e0c03136 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 16 Dec 2025 11:03:28 -0500 Subject: [PATCH 26/31] Finish scattered light module. --- docs/sbpy/surfaces.rst | 2 +- sbpy/surfaces/__init__.py | 2 +- sbpy/surfaces/scattered.py | 238 ++++++++++--------------------------- 3 files changed, 68 insertions(+), 174 deletions(-) diff --git a/docs/sbpy/surfaces.rst b/docs/sbpy/surfaces.rst index 67ab8dac..2191b0a0 100644 --- a/docs/sbpy/surfaces.rst +++ b/docs/sbpy/surfaces.rst @@ -86,7 +86,7 @@ Here, we define a new surface model with surface properties proportional to :mat >>> from sbpy.surfaces.surface import Surface, min_zero_cos >>> >>> class Cos2Surface(Surface): - ... """Absorption and emission proportional to :math:`\\cos^2`.""" + ... """Surface properties proportional to :math:`\\cos^2`.""" ... ... def absorption(self, epsilon, i): ... return epsilon * min_zero_cos(i)**2 diff --git a/sbpy/surfaces/__init__.py b/sbpy/surfaces/__init__.py index 6c360370..bb3bb9d1 100644 --- a/sbpy/surfaces/__init__.py +++ b/sbpy/surfaces/__init__.py @@ -1,3 +1,3 @@ from .surface import Surface from .lambertian import LambertianSurface -from .scattered import ScatteredLight +from .scattered import ScatteredLight, ScatteredSunlight diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py index afa1d3ed..55d072c6 100644 --- a/sbpy/surfaces/scattered.py +++ b/sbpy/surfaces/scattered.py @@ -1,18 +1,16 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst -# from abc import ABC, abstractmethod - -# import numpy as np - from typing import Optional + +import numpy as np import astropy.units as u +from ..calib import Sun +from ..units.typing import SpectralRadianceQuantity, UnitLike from ..spectroscopy.sources import SpectralSource, SinglePointSpectrumError from .surface import Surface -# from ..calib import Sun -# from ..units.typing import SpectralQuantity, SpectralFluxDensityQuantity, UnitLike -from ..units.typing import SpectralRadianceQuantity, UnitLike +__all__ = ["ScatteredLight", "ScatteredSunlight"] class ScatteredLight: @@ -93,189 +91,85 @@ def radiance( return F_i * self.surface.reflectance(albedo, i, e, phi) + def radiance_from_vectors( + self, + wfb, + albedo: u.physical.dimensionless, + n: np.ndarray, + r: u.physical.length, + ro: u.physical.length, + unit: Optional[UnitLike] = None, + ) -> SpectralRadianceQuantity: + """Vector-based alternative to `radiance`. -# class ScatteredLight(ABC): -# """Abstract base class to observe light scattered by a surface.""" - -# @u.quantity_input -# def scattered_light( -# self, -# wave_freq: SpectralQuantity, -# i: u.physical.angle, -# e: u.physical.angle, -# phi: u.physical.angle, -# unit: UnitLike = "W/(m2 sr um)", -# ) -> u.Quantity: -# """Radiance from light scattered by a surface. - - -# Parameters -# ---------- - -# wave_freq : `astropy.units.Quantity` -# Wavelength or frequency at which to evaluate the light source. - -# i : `~astropy.units.Quantity` -# Angle from normal of incident light. - -# e : `~astropy.units.Quantity` -# Angle from normal of emitted light. - -# phi : `~astropy.units.Quantity` -# Source-target-observer (phase) angle. - -# unit : `~astropy.units.Unit`, optional -# Unit of the return value. - - -# Returns -# ------- - -# radiance : `~astropy.units.Quantity` -# Observed radiance. - -# """ - -# @u.quantity_input -# def scattered_light_from_vectors( -# self, -# wave_freq: SpectralQuantity, -# n: np.ndarray, -# rs: u.physical.length, -# ro: u.physical.length, -# unit: UnitLike = "W/(m2 sr um)", -# ) -> u.Quantity: -# """Observed light reflected from a surface. - - -# Parameters -# ---------- - -# wave_freq : `astropy.units.Quantity` -# Wavelength or frequency at which to evaluate the light source. - -# n : `numpy.ndarray` -# Surface normal vector. - -# rs : `~astropy.units.Quantity` -# Radial vector from the surface to the light source. - -# ro : `~astropy.units.Quantity` -# Radial vector from the surface to the observer. - -# unit : `~astropy.units.Unit`, optional -# Unit of the return value. - - -# Returns -# ------- - -# radiance : `~astropy.units.Quantity` -# Observed radiance. - -# """ - - -# class ScatteredSunlight(ScatteredLight): -# """Observe sunlight scattered by a surface.""" - -# @u.quantity_input -# def scattered_light( -# self, -# wave_freq: SpectralQuantity, -# i: u.physical.angle, -# e: u.physical.angle, -# phi: u.physical.angle, -# rh: u.physical.length = 1 * u.au, -# unit: UnitLike = "W/(m2 sr um)", -# ) -> u.Quantity: -# """Radiance from sunlight scattered by a surface. - - -# Parameters -# ---------- - -# wave_freq : `astropy.units.Quantity` -# Wavelength or frequency at which to evaluate the Sun. Arrays are -# evaluated with `sbpy.calib.core.Sun.observe()`. - -# i : `~astropy.units.Quantity` -# Angle from normal of incident light. - -# e : `~astropy.units.Quantity` -# Angle from normal of emitted light. - -# phi : `~astropy.units.Quantity` -# Source-target-observer (phase) angle. - -# rh : `~astropy.units.Quantity` -# Heliocentric distance, default = 1 au. - -# unit : `~astropy.units.Unit`, optional -# Unit of the return value. + Parameters + ---------- + wfb : `~astropy.units.Quantity`, `~synphot.SpectralElement`, list + Wavelengths, frequencies, bandpass, or list of + bandpasses of the observation. Bandpasses require + `~synphot`. -# Returns -# ------- -# radiance : `~astropy.units.Quantity` -# Observed radiance. + albedo : `~astropy.units.Quantity` + Surface albedo for given `wfb` -# """ + n : `numpy.ndarray` + Surface normal vector. -# sun = Sun.from_default() -# flux_density_unit = u.Unit(unit) * u.sr -# if wave_freq.size == 1: -# F_i = sun(wave_freq, unit=flux_density_unit) -# else: -# F_i = sun.observe(wave_freq, unit=flux_density_unit) + r : `~astropy.units.Quantity` + Radial vector from the surface to the light source. -# F_i /= rh.to_value("au") ** 2 -# return self.reflectance(F_i, i, e, phi).to(unit) + ro : `~astropy.units.Quantity` or ``None`` + Radial vector from the surface to the observer. Parameter is unused + and may be ``None``. -# @u.quantity_input -# def scattered_light_from_vectors( -# self, -# wave_freq: SpectralQuantity, -# n: np.ndarray, -# rs: u.physical.length, -# ro: u.physical.length, -# unit: UnitLike = "W/(m2 sr um)", -# ) -> u.Quantity: -# """Observed sunlight reflected from a surface. + unit : str or `~astropy.units.Unit`, optional + Spectral radiance units of the result. -# Parameters -# ---------- + Returns + ------- + I : `~astropy.units.Quantity` -# wave_freq : `astropy.units.Quantity` -# Wavelength or frequency at which to evaluate the Sun. Arrays are -# evaluated with `sbpy.calib.core.Sun.observe()`. + """ -# n : `numpy.ndarray` -# Surface normal vector. + i = self._angle(n, r) + e = self._angle(n, ro) + phi = self._angle(r, ro) + return self.radiance(wfb, albedo, i, e, phi, unit=unit) -# rs : `~astropy.units.Quantity` -# Radial vector from the surface to the light source. -# ro : `~astropy.units.Quantity` -# Radial vector from the surface to the observer. +class ScatteredSunlight(ScatteredLight): + """Scatter sunlight off a surface. -# unit : `~astropy.units.Unit`, optional -# Unit of the return value. + Examples + -------- -# Returns -# ------- + Scatter sunlight from a Lambertian surface. -# radiance : `~astropy.units.Quantity` -# Observed radiance. + >>> import astropy.units as u + >>> from sbpy.surfaces import ScatteredSunlight, LambertianSurface + >>> + >>> surface = LambertianSurface() + >>> scattered = ScatteredSunlight(surface) + >>> + >>> wave = 550 * u.nm + >>> albedo = 0.1 + >>> i, e, phi = [30, 60, 90] * u.deg + >>> scattered.radiance(wave, albedo, i, e, phi) # doctest: +FLOAT_CMP + -# """ + The solar spectrum used is controlled with the ``sbpy.calib`` module:: -# rh = np.linalg.norm(rs).to("au") -# i, e, phi = self._vectors_to_angles(n, rs, ro) -# return self.scattered_sunlight(wave_freq, i, e, phi, rh=rh, unit=unit) + >>> from sbpy.calib import solar_spectrum + >>> with solar_spectrum.set("E490_2014LR"): + ... scattered = ScatteredSunlight(surface) + >>> scattered.radiance(wave, albedo, i, e, phi) # doctest: +FLOAT_CMP + + """ -# class LambertianSurfaceScatteredSunlight(LambertianSurface, ScatteredSunlight): -# """Sunlight scattered from a Lambertian surface.""" + def __init__(self, surface: Surface): + self.surface = surface + self.source = Sun.from_default() From d661f0a50ad06e72b6c72f8c7810993c41b46839 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 16 Dec 2025 12:16:07 -0500 Subject: [PATCH 27/31] Rename methods: absorptance and emittance --- docs/sbpy/surfaces.rst | 24 ++++++++--------- sbpy/surfaces/lambertian.py | 18 ++++++------- sbpy/surfaces/surface.py | 36 +++++++++++++------------- sbpy/surfaces/tests/test_lambertian.py | 18 ++++++------- 4 files changed, 48 insertions(+), 48 deletions(-) diff --git a/docs/sbpy/surfaces.rst b/docs/sbpy/surfaces.rst index 2191b0a0..27850817 100644 --- a/docs/sbpy/surfaces.rst +++ b/docs/sbpy/surfaces.rst @@ -8,14 +8,14 @@ Introduction The surfaces module is being made available on a preview basis. The API is subject to change. Feedback on the approach is welcome. -The ``surfaces`` module describes the interaction of electromagnetic radiation with surfaces. Sbpy uses the :math:`(i, e, \phi)` model (angle of incidence, angle of emittance, and phase angle) to describe how light scatters and emits light. It has a flexible system that can incorporate any surface scattering model that can be described with these three angles. +The ``surfaces`` module describes the interaction of electromagnetic radiation with surfaces. Sbpy uses the :math:`(i, e, \phi)` model (angle of incidence, angle of emittance, and phase angle) to describe how light scatters and emits light. It has a flexible system that can incorporate most surface scattering Models that can be described with these three angles. .. figure:: ../static/scattering-vectors.svg - :alt: Diagram of surface scattering and emission vectors + :alt: Diagram of surface scattering and emittance vectors - Sbpy's geometric basis for surface scattering and emission: :math:`n` is the surface normal vector, :math:`r_s` is the radial vector to the light source, and :math:`r_o` is the radial vector to the observer. The angle of incidence (:math:`i`), angle of emittance (:math:`e`), phase angle (:math:`\phi`) are labeled. + Sbpy's geometric basis for surface scattering and emittance: :math:`n` is the surface normal vector, :math:`r_s` is the radial vector to the light source, and :math:`r_o` is the radial vector to the observer. The angle of incidence (:math:`i`), angle of emittance (:math:`e`), phase angle (:math:`\phi`) are labeled. -An implementation of the ``Surface`` model has methods to calculate electromagnetic absorption, emission, and bidirectional reflectance. +An implementation of the ``Surface`` model has methods to calculate electromagnetic absorptance, emittance, and bidirectional reflectance. Getting Started @@ -23,7 +23,7 @@ Getting Started Currently the Lambertian surface model is implemented. A Lambertian surface absorbs and emits light uniformly in all directions. -Create an instance of the ``LambertianSurface`` model, and calculate the absorption and emission scale factors for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)`. Let the albedo be 0.1 (emissivity = 0.9):: +Create an instance of the ``LambertianSurface`` model, and calculate the absorptance and emittance scale factors for :math:`(i, e, \phi) = (30^\circ, 60^\circ, 90^\circ)`. Let the albedo be 0.1 (emissivity = 0.9):: >>> import astropy.units as u >>> from sbpy.surfaces import LambertianSurface @@ -33,9 +33,9 @@ Create an instance of the ``LambertianSurface`` model, and calculate the absorpt >>> i, e, phi = [30, 60, 90] * u.deg >>> >>> surface = LambertianSurface() - >>> surface.absorption(epsilon, i) # doctest: +FLOAT_CMP + >>> surface.absorptance(epsilon, i) # doctest: +FLOAT_CMP - >>> surface.emission(epsilon, e, phi) # doctest: +FLOAT_CMP + >>> surface.emittance(epsilon, e, phi) # doctest: +FLOAT_CMP Calculate the bidirectional reflectance for :math:`e=0`, and a range of incident angles:: @@ -78,7 +78,7 @@ As an alternative to using :math:`(i, e, \phi)`, results may be calculated using Build Your Own Surface Models ----------------------------- -To define your own surface model create a new class based on `~sbpy.surfaces.surface.Surface`, and define the methods for ``absorption``, ``emission``, and ``reflectance``. The `~sbpy.surfaces.lambertian.LambertianSurface` model may serve as a reference. +To define your own surface model create a new class based on `~sbpy.surfaces.surface.Surface`, and define the methods for ``absorptance``, ``emittance``, and ``reflectance``. The `~sbpy.surfaces.lambertian.LambertianSurface` model may serve as a reference. Here, we define a new surface model with surface properties proportional to :math:`\cos^2`:: @@ -88,10 +88,10 @@ Here, we define a new surface model with surface properties proportional to :mat >>> class Cos2Surface(Surface): ... """Surface properties proportional to :math:`\\cos^2`.""" ... - ... def absorption(self, epsilon, i): + ... def absorptance(self, epsilon, i): ... return epsilon * min_zero_cos(i)**2 ... - ... def emission(self, epsilon, e, phi): + ... def emittance(self, epsilon, e, phi): ... return epsilon * min_zero_cos(e)**2 ... ... def reflectance(self, albedo, i, e, phi): @@ -104,9 +104,9 @@ Create and use an instance of our new model:: >>> epsilon = 1 - albedo >>> i, e, phi = [30, 60, 90] * u.deg >>> - >>> surface.absorption(epsilon, i) # doctest: +FLOAT_CMP + >>> surface.absorptance(epsilon, i) # doctest: +FLOAT_CMP - >>> surface.emission(epsilon, e, phi) # doctest: +FLOAT_CMP + >>> surface.emittance(epsilon, e, phi) # doctest: +FLOAT_CMP >>> surface.reflectance(albedo, i, e, phi) # doctest: +FLOAT_CMP diff --git a/sbpy/surfaces/lambertian.py b/sbpy/surfaces/lambertian.py index d45f693b..a516f056 100644 --- a/sbpy/surfaces/lambertian.py +++ b/sbpy/surfaces/lambertian.py @@ -9,7 +9,7 @@ class LambertianSurface(Surface): - """Lambertian surface absorption, emission, and reflectance. + """Lambertian surface absorptance, emittance, and reflectance. The surface is illuminated at an angle of :math:`i`, and emitted toward an angle of :math:`e`, both measured from the surface normal direction. @@ -20,7 +20,7 @@ class LambertianSurface(Surface): Examples -------- - Absorption of light for a surface with an albedo of 0.1 and an emissivity of + Absorptance of light for a surface with an albedo of 0.1 and an emissivity of 0.9: >>> import astropy.units as u @@ -31,12 +31,12 @@ class LambertianSurface(Surface): >>> epsilon = 1 - albedo >>> i, e, phi = [30, 60, 90] * u.deg >>> - >>> surface.absorption(epsilon, i) # doctest: +FLOAT_CMP + >>> surface.absorptance(epsilon, i) # doctest: +FLOAT_CMP - Emission: + Emittance: - >>> surface.emission(epsilon, e, phi) # doctest: +FLOAT_CMP + >>> surface.emittance(epsilon, e, phi) # doctest: +FLOAT_CMP Bidirectional reflectance: @@ -55,7 +55,7 @@ class LambertianSurface(Surface): """ @u.quantity_input - def absorption( + def absorptance( self, epsilon: float, i: u.physical.angle, @@ -64,7 +64,7 @@ def absorption( return epsilon * min_zero_cos(i) @u.quantity_input - def emission( + def emittance( self, epsilon: float, e: u.physical.angle, @@ -84,7 +84,7 @@ def reflectance( return albedo * min_zero_cos(i) * min_zero_cos(e) / np.pi / u.sr @u.quantity_input - def emission_from_vectors( + def emittance_from_vectors( self, epsilon: float, n: np.ndarray, @@ -92,4 +92,4 @@ def emission_from_vectors( ro: u.physical.length, ) -> u.Quantity[u.dimensionless_unscaled]: e = self._angle(n, ro) - return self.emission(epsilon, e, None) + return self.emittance(epsilon, e, None) diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index c999dd6e..a024fb3d 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -24,10 +24,10 @@ class Surface(ABC): """Abstract base class for all small-body surfaces.""" @abstractmethod - def absorption( + def absorptance( self, epsilon: u.physical.dimensionless, i: u.physical.angle ) -> u.Quantity[u.dimensionless_unscaled]: - r"""Absorption of directional, incident light. + r"""Absorptance of directional, incident light. The surface is illuminated at an angle of :math:`i`, measured from the surface normal direction. @@ -45,21 +45,21 @@ def absorption( Returns ------- a : `~astropy.units.Quantity` - Fraction of incident light absorbed. + Absorptance. """ @abstractmethod - def emission( + def emittance( self, epsilon: u.physical.dimensionless, e: u.physical.angle, phi: u.physical.angle, ) -> u.Quantity[u.dimensionless_unscaled]: - r"""Emission of directional light from a surface. + r"""Emittance of directional light from a surface. The surface is observed at an angle of :math:`e`, measured from the - surface normal direction. Anisotropic emission is characterized by the + surface normal direction. Anisotropic emittance is characterized by the angle `phi`. @@ -72,13 +72,13 @@ def emission( Observed angle from normal. phi : `~astropy.units.Quantity` - Angle to account for anisotropic emission. + Angle to account for anisotropic emittance. Returns ------- em : `~astropy.units.Quantity` - Emission factor. + Emittance. """ @@ -126,14 +126,14 @@ def _angle(a: u.Quantity, b: u.Quantity) -> u.physical.angle: return u.Quantity(np.arccos(np.dot(a_hat, b_hat)), "rad") @u.quantity_input - def absorption_from_vectors( + def absorptance_from_vectors( self, epsilon: u.physical.dimensionless, n: np.ndarray, r: u.physical.length, ro: Union[u.physical.length, None], ) -> u.Quantity[u.dimensionless_unscaled]: - """Vector-based alternative to `absorption`. + """Vector-based alternative to `absorptance`. Input vectors do not need to be normalized. @@ -156,23 +156,23 @@ def absorption_from_vectors( Returns ------- - F_a : `~astropy.units.Quantity` - Absorbed spectral flux density. + a : `~astropy.units.Quantity` + Absorptance. """ i = self._angle(n, r) - return self.absorption(epsilon, i) + return self.absorptance(epsilon, i) @u.quantity_input - def emission_from_vectors( + def emittance_from_vectors( self, epsilon: u.physical.dimensionless, n: np.ndarray, r: u.physical.length, ro: u.physical.length, ) -> u.Quantity[u.dimensionless_unscaled]: - r"""Vector-based alternative to `emission`. + r"""Vector-based alternative to `emittance`. Input vectors do not need to be normalized. @@ -194,14 +194,14 @@ def emission_from_vectors( Returns ------- - F_e : `~astropy.units.Quantity` - Spectral radiance / specific intensity received by the observer. + em : `~astropy.units.Quantity` + Emittance. """ e = self._angle(n, ro) phi = self._angle(r, ro) - return self.emission(epsilon, e, phi) + return self.emittance(epsilon, e, phi) @u.quantity_input def reflectance_from_vectors( diff --git a/sbpy/surfaces/tests/test_lambertian.py b/sbpy/surfaces/tests/test_lambertian.py index fcaff0ba..a542fd8c 100644 --- a/sbpy/surfaces/tests/test_lambertian.py +++ b/sbpy/surfaces/tests/test_lambertian.py @@ -7,18 +7,18 @@ from ..lambertian import LambertianSurface -def test_absorption(): +def test_absorptance(): epsilon = 0.9 i = Angle([0, 30, 45, 60, 90, 100], "deg") expected = 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) surface = LambertianSurface() - result = surface.absorption(epsilon, i) + result = surface.absorptance(epsilon, i) assert u.allclose(result, expected) -def test_absorption_from_vectors(): +def test_absorptance_from_vectors(): epsilon = 0.9 # equivalent to (i, e, phi) = (30, 60, 90) deg n = [1, 0, 0] @@ -28,22 +28,22 @@ def test_absorption_from_vectors(): expected = 0.9 * np.sqrt(3) / 2 surface = LambertianSurface() - result = surface.absorption_from_vectors(epsilon, n, r, ro) + result = surface.absorptance_from_vectors(epsilon, n, r, ro) assert u.allclose(result, expected) -def test_emission(): +def test_emittance(): epsilon = 0.9 e = Angle([0, 30, 45, 60, 90, 100], "deg") expected = 0.9 * np.array([1, np.sqrt(3) / 2, 1 / np.sqrt(2), 0.5, 0, 0]) surface = LambertianSurface() - result = surface.emission(epsilon, e, None) + result = surface.emittance(epsilon, e, None) assert u.allclose(result, expected) -def test_emission_from_vectors(): +def test_emittance_from_vectors(): epsilon = 0.9 # equivalent to (i, e, phi) = (30, 60, 90) deg n = [1, 0, 0] @@ -53,11 +53,11 @@ def test_emission_from_vectors(): expected = 0.45 surface = LambertianSurface() - result = surface.emission_from_vectors(epsilon, n, r, ro) + result = surface.emittance_from_vectors(epsilon, n, r, ro) assert u.allclose(result, expected) # incident vector is optional - result = surface.emission_from_vectors(epsilon, n, None, ro) + result = surface.emittance_from_vectors(epsilon, n, None, ro) assert u.allclose(result, expected) From 76dbe9fd2dfd98f781bf6d9281da4738890ee69e Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 5 Jan 2026 12:03:47 -0500 Subject: [PATCH 28/31] Start defining Sphere --- sbpy/shape/core.py | 152 +++++++++++++++++++++++++++++++-------------- 1 file changed, 106 insertions(+), 46 deletions(-) diff --git a/sbpy/shape/core.py b/sbpy/shape/core.py index 95a8c6be..cd98b472 100644 --- a/sbpy/shape/core.py +++ b/sbpy/shape/core.py @@ -1,72 +1,132 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst """ -sbpy Shape Module - -created on June 26, 2017 +Asteroid and comet nucleus shape models """ -__all__ = ['ModelClass', 'Kaasalainen', 'Lightcurve'] +from typing import Callable +import astropy.units as u + +from ..data.ephem import Ephem +from ..data.decorators import dataclass_input + + +class Shape: + """Model asteroid or comet nucleus shape.""" + + +class Sphere(Shape): + """A spherical object. + + + Parameters + ---------- + + + """ + + def __init__(self, radius: u.physical.length): + self.radius = radius + + def to_faceted_model(self): + pass + + @dataclass_input + def integrate_over_surface( + self, func: Callable, eph: Ephem, *args, **kwargs + ) -> u.Quantity: + """Integrate the function over the surface. + + The integration is over the observed area:: + + + + + Parameters + ---------- + func : callable + The function to integrate: ``func(*args, i, e, phi, **kwargs)``, + where :math:`i` is the angle of incidence, :math:`e` is the the + angle of emittance, and :math:`phi` is the phase angle. + + eph : `Ephem` + The observing geometry as an ephemeris (or equivalent) object: + - rh + - delta + - phase + + *args, **kwargs + Additional arguments passed to the function. + + + Returns + ------- + total : `~sbpy.units.Quantity` + + """ + + +# __all__ = ["ModelClass", "Kaasalainen", "Lightcurve"] -class ModelClass(): +# class ModelClass: - def __init__(self): - self.shape = None +# def __init__(self): +# self.shape = None - def load_obj(self, filename): - """Load .OBJ shape model""" +# def load_obj(self, filename): +# """Load .OBJ shape model""" - def get_facet(self, facetid): - """Retrieve information for one specific surface facet""" +# def get_facet(self, facetid): +# """Retrieve information for one specific surface facet""" - def iof(self, facetid, eph): - """Derive I/F for one specific surface facet""" +# def iof(self, facetid, eph): +# """Derive I/F for one specific surface facet""" - def integrated_flux(self, eph): - """Derive surface-integrated flux""" +# def integrated_flux(self, eph): +# """Derive surface-integrated flux""" - def lightcurve(self, eph, time): - """Derive lightcurve""" +# def lightcurve(self, eph, time): +# """Derive lightcurve""" -class Kaasalainen(ModelClass): +# class Kaasalainen(ModelClass): - def __init__(self): - self.properties = None - # setup model properties +# def __init__(self): +# self.properties = None +# # setup model properties - def convexinv(self, eph): - """Obtain shape model through lightcurve inversion, optimizing all - parameters and uses spherical harmonics functions for shape - representation""" +# def convexinv(self, eph): +# """Obtain shape model through lightcurve inversion, optimizing all +# parameters and uses spherical harmonics functions for shape +# representation""" - def conjgradinv(self, eph): - """Obtain shape model through lightcurve inversion, optimizing only - shape and uses directly facet areas as parameters""" +# def conjgradinv(self, eph): +# """Obtain shape model through lightcurve inversion, optimizing only +# shape and uses directly facet areas as parameters""" -class Lightcurve(): +# class Lightcurve: - def __init__(self, eph): - self.eph = eph - self.fouriercoeff = None - self.period = None - self.pole = (0, 90) # ecliptic coordinates +# def __init__(self, eph): +# self.eph = eph +# self.fouriercoeff = None +# self.period = None +# self.pole = (0, 90) # ecliptic coordinates - def axis_ratio(self): - """Derive axis ratio from lightcurve amplitude""" +# def axis_ratio(self): +# """Derive axis ratio from lightcurve amplitude""" - def derive_period(self, method='lomb-scargle'): - """Derive lightcurve period using different methods""" +# def derive_period(self, method="lomb-scargle"): +# """Derive lightcurve period using different methods""" - def fit_fourier(self, order): - """Fit Fourier coefficients to lightcurve""" +# def fit_fourier(self, order): +# """Fit Fourier coefficients to lightcurve""" - def fit_pole(self): - """Fit pole orientation""" +# def fit_pole(self): +# """Fit pole orientation""" - def fit(self): - """Fit period, pole orientation, and Fourier coefficients at the same time""" +# def fit(self): +# """Fit period, pole orientation, and Fourier coefficients at the same time""" - def simulate(self): - """Simulate a lightcurve from period, Fourier coefficients, pole orientation""" +# def simulate(self): +# """Simulate a lightcurve from period, Fourier coefficients, pole orientation""" From 8ce0874e5fcdbda6e6629fa1dd432e591d65bb72 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 23 Dec 2025 12:34:05 -0500 Subject: [PATCH 29/31] Remove "type" --- sbpy/units/typing.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/sbpy/units/typing.py b/sbpy/units/typing.py index 96f20c75..e06b3962 100644 --- a/sbpy/units/typing.py +++ b/sbpy/units/typing.py @@ -7,22 +7,22 @@ import astropy.units as u from astropy import __version__ as astropy_version -type UnitLike = Union[str, u.Unit] -type SpectralQuantity = Union[ +UnitLike = Union[str, u.Unit] +SpectralQuantity = Union[ u.Quantity[u.physical.length], u.Quantity[u.physical.frequency] ] -type SpectralFluxDensityQuantity = Union[ +SpectralFluxDensityQuantity = Union[ u.Quantity[u.physical.spectral_flux_density], u.Quantity[u.physical.spectral_flux_density_wav], ] if Version(astropy_version) < Version("6.1"): - type SpectralRadianceQuantity = Union[ + SpectralRadianceQuantity = Union[ u.Quantity[u.physical.spectral_flux_density / u.sr], u.Quantity[u.physical.spectral_flux_density_wav / u.sr], ] else: - type SpectralRadianceQuantity = Union[ + SpectralRadianceQuantity = Union[ u.Quantity[u.physical.surface_brightness], u.Quantity[u.physical.surface_brightness_wav], ] From 7515e3a70ba98d9e091c6303dbc177268bb57d8d Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 23 Dec 2025 12:34:20 -0500 Subject: [PATCH 30/31] Use a substitutions file for docs so that astropy docstrings can be used --- docs/conf.py | 96 ++++++++++++++++++++++-------------------- docs/substitutions.txt | 34 +++++++++++++++ 2 files changed, 84 insertions(+), 46 deletions(-) create mode 100644 docs/substitutions.txt diff --git a/docs/conf.py b/docs/conf.py index a71dc65e..85b64215 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -30,26 +30,30 @@ import sys import matplotlib -matplotlib.use('agg') + +matplotlib.use("agg") try: from sphinx_astropy.conf.v1 import * except ImportError: - print('ERROR: the documentation requires the sphinx-astropy package to be' - ' installed') + print( + "ERROR: the documentation requires the sphinx-astropy package to be" + " installed" + ) sys.exit(1) # Get configuration information from setup.cfg from configparser import ConfigParser # noqa: E402 + conf = ConfigParser() -conf.read([os.path.join(os.path.dirname(__file__), '..', 'setup.cfg')]) -setup_cfg = dict(conf.items('metadata')) +conf.read([os.path.join(os.path.dirname(__file__), "..", "setup.cfg")]) +setup_cfg = dict(conf.items("metadata")) # -- General configuration ---------------------------------------------------- # By default, highlight as Python 3. -highlight_language = 'python3' +highlight_language = "python3" # If your documentation needs a minimal Sphinx version, state it here. # needs_sphinx = '1.2' @@ -60,16 +64,18 @@ # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. -exclude_patterns.append('_templates') +exclude_patterns.append("_templates") # This is added to the end of RST files - a good place to put substitutions to # be used globally. -rst_epilog += """ -""" +with open("substitutions.txt") as inf: + rst_epilog += inf.read() -extensions += ['sphinx.ext.intersphinx', - 'sphinx_automodapi.smart_resolver', - 'sphinx.ext.autosectionlabel'] +extensions += [ + "sphinx.ext.intersphinx", + "sphinx_automodapi.smart_resolver", + "sphinx.ext.autosectionlabel", +] # For example, index:Introduction for a section called Introduction that # appears in document index.rst. @@ -83,20 +89,19 @@ # -- Project information ------------------------------------------------------ # This does not *have* to match the package name, but typically does -project = setup_cfg['package_name'] -author = setup_cfg['author'] -copyright = '{0}, {1}'.format( - datetime.datetime.now().year, setup_cfg['author']) +project = setup_cfg["package_name"] +author = setup_cfg["author"] +copyright = "{0}, {1}".format(datetime.datetime.now().year, setup_cfg["author"]) # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the # built documents. -__import__(setup_cfg['package_name']) -package = sys.modules[setup_cfg['package_name']] +__import__(setup_cfg["package_name"]) +package = sys.modules[setup_cfg["package_name"]] # The short X.Y version. -version = package.__version__.split('-', 1)[0] +version = package.__version__.split("-", 1)[0] # The full version, including alpha/beta/rc tags. release = package.__version__ @@ -113,68 +118,68 @@ # Add any paths that contain custom themes here, relative to this directory. # To use a different custom theme, add the directory containing the theme. -#html_theme_path = [] +# html_theme_path = [] # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. To override the custom theme, set this to the # name of a builtin theme or the name of a custom theme in html_theme_path. -#html_theme = 'sphinx_rtd_theme' +# html_theme = 'sphinx_rtd_theme' # Please update these texts to match the name of your package. html_theme_options = { - 'logotext1': 'sb', # white, semi-bold - 'logotext2': 'py', # orange, light - 'logotext3': ':docs' # white, light + "logotext1": "sb", # white, semi-bold + "logotext2": "py", # orange, light + "logotext3": ":docs", # white, light } # Custom sidebar templates, maps document names to template names. -#html_sidebars = {} +# html_sidebars = {} # The name of an image file (relative to this directory) to place at the top # of the sidebar. -#html_logo = '' +# html_logo = '' # The name of an image file (within the static path) to use as favicon of the # docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 # pixels large. -#html_favicon = '' +# html_favicon = '' # If not '', a 'Last updated on:' timestamp is inserted at every page bottom, # using the given strftime format. -#html_last_updated_fmt = '' +# html_last_updated_fmt = '' # The name for this set of Sphinx documents. If None, it defaults to # " v documentation". -html_title = '{0} v{1}'.format(project, release) +html_title = "{0} v{1}".format(project, release) # Output file base name for HTML help builder. -htmlhelp_basename = project + 'doc' +htmlhelp_basename = project + "doc" # -- Options for LaTeX output ------------------------------------------------- # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, author, documentclass [howto/manual]). -latex_documents = [('index', project + '.tex', project + u' Documentation', - author, 'manual')] +latex_documents = [ + ("index", project + ".tex", project + " Documentation", author, "manual") +] # -- Options for manual page output ------------------------------------------- # One entry per manual page. List of tuples # (source start file, name, description, authors, manual section). -man_pages = [('index', project.lower(), project + u' Documentation', - [author], 1)] +man_pages = [("index", project.lower(), project + " Documentation", [author], 1)] # -- Options for the edit_on_github extension --------------------------------- -if eval(setup_cfg.get('edit_on_github')): - extensions += ['sphinx_astropy.ext.edit_on_github'] +if eval(setup_cfg.get("edit_on_github")): + extensions += ["sphinx_astropy.ext.edit_on_github"] - versionmod = __import__(setup_cfg['package_name'] + '.version') - edit_on_github_project = setup_cfg['github_project'] + versionmod = __import__(setup_cfg["package_name"] + ".version") + edit_on_github_project = setup_cfg["github_project"] if versionmod.version.release: edit_on_github_branch = "v" + versionmod.version.version else: @@ -184,18 +189,17 @@ edit_on_github_doc_root = "docs" # -- Resolving issue number to links in changelog ----------------------------- -github_issues_url = 'https://github.com/{0}/issues/'.format( - setup_cfg['github_project']) +github_issues_url = "https://github.com/{0}/issues/".format(setup_cfg["github_project"]) # -- compile list of field names # import compile_fieldnames # --- intersphinx setup -intersphinx_mapping['astroquery'] = ( - 'https://astroquery.readthedocs.io/en/stable/', None) +intersphinx_mapping["astroquery"] = ( + "https://astroquery.readthedocs.io/en/stable/", + None, +) -intersphinx_mapping['synphot'] = ( - 'https://synphot.readthedocs.io/en/stable/', None) +intersphinx_mapping["synphot"] = ("https://synphot.readthedocs.io/en/stable/", None) -intersphinx_mapping['astropy'] = ( - 'https://docs.astropy.org/en/stable/', None) +intersphinx_mapping["astropy"] = ("https://docs.astropy.org/en/stable/", None) diff --git a/docs/substitutions.txt b/docs/substitutions.txt new file mode 100644 index 00000000..0becf7d8 --- /dev/null +++ b/docs/substitutions.txt @@ -0,0 +1,34 @@ +.. ReST substitutions and links for the docs and docstrings. + +.. This file is included in the conf.py rst_epilog variable. + +.. Keep common items aligned with the astropy substitutions +.. (astropy/docs/common_links.rst) to avoid issues with inheriting +.. substitutions from astropy classes + +.. NumPy +.. |ndarray| replace:: :class:`numpy.ndarray` + +.. Astropy +.. Coordinates +.. |Angle| replace:: `~astropy.coordinates.Angle` +.. |Latitude| replace:: `~astropy.coordinates.Latitude` +.. |Longitude| replace:: :class:`~astropy.coordinates.Longitude` +.. |BaseFrame| replace:: `~astropy.coordinates.BaseCoordinateFrame` +.. |SkyCoord| replace:: :class:`~astropy.coordinates.SkyCoord` + +.. Table +.. |Column| replace:: :class:`~astropy.table.Column` +.. |MaskedColumn| replace:: :class:`~astropy.table.MaskedColumn` +.. |TableColumns| replace:: :class:`~astropy.table.TableColumns` +.. |Row| replace:: :class:`~astropy.table.Row` +.. |Table| replace:: :class:`~astropy.table.Table` +.. |QTable| replace:: :class:`~astropy.table.QTable` + +.. Time +.. |Time| replace:: :class:`~astropy.time.Time` + +.. Units +.. |PhysicalType| replace:: :class:`~astropy.units.PhysicalType` +.. |Quantity| replace:: :class:`~astropy.units.Quantity` +.. |Unit| replace:: :class:`~astropy.units.UnitBase` From 3adda27cc04918ec39f63ec52ae0c0ef0293f158 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 23 Dec 2025 12:34:50 -0500 Subject: [PATCH 31/31] Unused imports --- sbpy/surfaces/surface.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py index a024fb3d..2a60cf4b 100644 --- a/sbpy/surfaces/surface.py +++ b/sbpy/surfaces/surface.py @@ -6,6 +6,9 @@ import numpy as np from astropy import units as u +# for bidirectional_reflectance +from ..units import physical # noqa: F401 + def min_zero_cos(a: u.physical.angle) -> u.Quantity[u.dimensionless_unscaled]: """Use to ensure that cos(>=90 deg) equals 0."""