From 30edfc954fbd717efbcc737914527786f0ae320c Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Sat, 9 Nov 2024 08:33:08 -0500 Subject: [PATCH 01/20] 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 6ae456a8131e1564ffe3a9712eaaef21ee433733 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 08:40:54 -0500 Subject: [PATCH 02/20] 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 050c1b01..2d02d42c 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 826bda5653d9d27838474e0d114e5fabcaf94891 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 09:44:30 -0500 Subject: [PATCH 03/20] 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 ef171c771a2f2b1154c8ee7d9bfdf97de9f9606c Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 12:58:48 -0500 Subject: [PATCH 04/20] 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 abece2bdefc0adcef296d80f447e86346e4a69d0 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 13:11:31 -0500 Subject: [PATCH 05/20] 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 bc77d7fad846b0e8b730e4aaa5529de6dce0a116 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 13:16:10 -0500 Subject: [PATCH 06/20] 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 c0a910f74a35afbb4fcaeac4bc1996163ab9dabf Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 13:19:40 -0500 Subject: [PATCH 07/20] 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 22a77fd415c795b3688dd8205b1cd46042a4488d Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Mon, 11 Nov 2024 13:24:15 -0500 Subject: [PATCH 08/20] 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 abd03e8d915b9f0b331d620a5d174132a4034d40 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 10:46:51 -0500 Subject: [PATCH 09/20] 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 ec78c1a2545bed635f3af7432d36807456ec10d2 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:38:31 -0500 Subject: [PATCH 10/20] 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 8dc08ebb96e1e6c4e3985f42632b9634f42a2b40 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:39:15 -0500 Subject: [PATCH 11/20] 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 8185f25a32442f26d2a716c9d31852b92e64925b Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:39:31 -0500 Subject: [PATCH 12/20] 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 7749336fec6e9d42f3ac336353630b0a91d8f2a8 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:40:03 -0500 Subject: [PATCH 13/20] 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 47679668ffa0ae0e4ccde102eea1e4e9859d5489 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:55:08 -0500 Subject: [PATCH 14/20] 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 2689e8e8cd01cdeed57076c8f951ffbf77a3a6ed Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 12:59:38 -0500 Subject: [PATCH 15/20] 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 3166910253379c14bfd1add2033b2157d82a0dce Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 13:30:46 -0500 Subject: [PATCH 16/20] 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 191e53cc0941985b13f310ce6ab7b1ac28c37e46 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 12 Nov 2024 16:23:31 -0500 Subject: [PATCH 17/20] 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 15ef0810c385a9f49cc5ddec5e5caf8f1e56dd9a Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 26 Nov 2024 22:49:03 -0500 Subject: [PATCH 18/20] 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 8461d791f3fe54264702c9994794189afeb0dfad Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 10 Dec 2024 09:50:03 -0500 Subject: [PATCH 19/20] 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 7b8b265a159ed31fe382c437001c587ec9353b80 Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Fri, 5 Sep 2025 09:20:13 -0400 Subject: [PATCH 20/20] 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)