diff --git a/docs/conf.py b/docs/conf.py index a71dc65eb..85b64215a 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/index.rst b/docs/index.rst index 33c62534b..129213d40 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 000000000..27850817f --- /dev/null +++ b/docs/sbpy/surfaces.rst @@ -0,0 +1,119 @@ +Surfaces Module (`sbpy.surfaces`) +================================= + +Introduction +------------ + +.. admonition:: warning + + 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 most surface scattering Models that can be described with these three angles. + +.. figure:: ../static/scattering-vectors.svg + :alt: Diagram of surface scattering and emittance vectors + + 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 absorptance, emittance, and bidirectional reflectance. + + +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 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 + >>> + >>> albedo = 0.1 + >>> epsilon = 1 - albedo + >>> i, e, phi = [30, 60, 90] * u.deg + >>> + >>> surface = LambertianSurface() + >>> surface.absorptance(epsilon, i) # 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:: + +.. plot:: + :context: reset + :nofigs: + + >>> import numpy as np + >>> + >>> e = 0 * u.deg + >>> i = np.linspace(-90, 90) * u.deg + >>> phi = np.abs(e - i) # calculate phase angle + >>> r = surface.reflectance(albedo, i, e, phi) + +.. plot:: + :include-source: + :context: + + import matplotlib.pyplot as plt + fig, ax = plt.subplots() + ax.plot(i, r) + 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 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(albedo, 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 ``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`:: + + >>> # use min_zero_cos(a) to ensure cos(a >= 90 deg) = 0 + >>> from sbpy.surfaces.surface import Surface, min_zero_cos + >>> + >>> class Cos2Surface(Surface): + ... """Surface properties proportional to :math:`\\cos^2`.""" + ... + ... def absorptance(self, epsilon, i): + ... return epsilon * min_zero_cos(i)**2 + ... + ... def emittance(self, epsilon, e, phi): + ... return epsilon * min_zero_cos(e)**2 + ... + ... 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.absorptance(epsilon, i) # doctest: +FLOAT_CMP + + >>> surface.emittance(epsilon, e, phi) # doctest: +FLOAT_CMP + + >>> surface.reflectance(albedo, i, e, phi) # doctest: +FLOAT_CMP + + + +Reference/API +------------- +.. automodapi:: sbpy.surfaces + :no-heading: + :inherited-members: diff --git a/docs/static/scattering-vectors.svg b/docs/static/scattering-vectors.svg new file mode 100644 index 000000000..0e437bea2 --- /dev/null +++ b/docs/static/scattering-vectors.svg @@ -0,0 +1,127 @@ + + + + + + + + + + + + + + + + + + i + e + φ + n + rs + ro + + diff --git a/docs/substitutions.txt b/docs/substitutions.txt new file mode 100644 index 000000000..0becf7d85 --- /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` diff --git a/sbpy/shape/core.py b/sbpy/shape/core.py index 95a8c6be6..cd98b4729 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""" diff --git a/sbpy/surfaces/__init__.py b/sbpy/surfaces/__init__.py new file mode 100644 index 000000000..bb3bb9d19 --- /dev/null +++ b/sbpy/surfaces/__init__.py @@ -0,0 +1,3 @@ +from .surface import Surface +from .lambertian import LambertianSurface +from .scattered import ScatteredLight, ScatteredSunlight diff --git a/sbpy/surfaces/lambertian.py b/sbpy/surfaces/lambertian.py new file mode 100644 index 000000000..a516f0561 --- /dev/null +++ b/sbpy/surfaces/lambertian.py @@ -0,0 +1,95 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +from typing import Union +import numpy as np +import astropy.units as u + +from .surface import Surface, min_zero_cos +from ..units.typing import SpectralFluxDensityQuantity, SpectralRadianceQuantity + + +class LambertianSurface(Surface): + """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. + :math:`\\phi` is the source-target-observer (phase) angle. Both the source + and the emitted light are assumed to be collimated. + + + Examples + -------- + + Absorptance 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.absorptance(epsilon, i) # doctest: +FLOAT_CMP + + + Emittance: + + >>> surface.emittance(epsilon, e, phi) # doctest: +FLOAT_CMP + + + Bidirectional reflectance: + + >>> 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(albedo, n, r, ro) # doctest: +FLOAT_CMP + + + """ + + @u.quantity_input + def absorptance( + self, + epsilon: float, + i: u.physical.angle, + ) -> u.Quantity[u.dimensionless_unscaled]: + # use min_zero_cos(i) to ensure cos(>= 90 deg) = 0 + return epsilon * min_zero_cos(i) + + @u.quantity_input + def emittance( + 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 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 albedo * min_zero_cos(i) * min_zero_cos(e) / np.pi / u.sr + + @u.quantity_input + def emittance_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.emittance(epsilon, e, None) diff --git a/sbpy/surfaces/scattered.py b/sbpy/surfaces/scattered.py new file mode 100644 index 000000000..55d072c6a --- /dev/null +++ b/sbpy/surfaces/scattered.py @@ -0,0 +1,175 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +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 + +__all__ = ["ScatteredLight", "ScatteredSunlight"] + + +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) + + 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`. + + + 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` + + n : `numpy.ndarray` + Surface normal vector. + + 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``. + + unit : str or `~astropy.units.Unit`, optional + Spectral radiance units of the result. + + + Returns + ------- + I : `~astropy.units.Quantity` + + """ + + 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) + + +class ScatteredSunlight(ScatteredLight): + """Scatter sunlight off a surface. + + + Examples + -------- + + Scatter sunlight from a Lambertian surface. + + >>> 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:: + + >>> 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 + + + """ + + def __init__(self, surface: Surface): + self.surface = surface + self.source = Sun.from_default() diff --git a/sbpy/surfaces/surface.py b/sbpy/surfaces/surface.py new file mode 100644 index 000000000..2a60cf4b7 --- /dev/null +++ b/sbpy/surfaces/surface.py @@ -0,0 +1,247 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +from abc import ABC, abstractmethod +from typing import Union + +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.""" + + # handle scalars separately + if a.ndim == 0 and u.isclose(np.abs(a), 90 * u.deg): + return u.Quantity(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.""" + + @abstractmethod + def absorptance( + self, epsilon: u.physical.dimensionless, i: u.physical.angle + ) -> u.Quantity[u.dimensionless_unscaled]: + r"""Absorptance of directional, incident light. + + The surface is illuminated at an angle of :math:`i`, measured from the + surface normal direction. + + + Parameters + ---------- + epsilon : `~astropy.units.Quantity` + Surface emissivity. + + i : `~astropy.units.Quantity` + Angle from normal of incident light. + + + Returns + ------- + a : `~astropy.units.Quantity` + Absorptance. + + """ + + @abstractmethod + def emittance( + self, + epsilon: u.physical.dimensionless, + e: u.physical.angle, + phi: u.physical.angle, + ) -> u.Quantity[u.dimensionless_unscaled]: + 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 emittance is characterized by the + angle `phi`. + + + Parameters + ---------- + epsilon : `~astropy.units.Quantity` + Surface emissivity. + + e : `~astropy.units.Quantity` + Observed angle from normal. + + phi : `~astropy.units.Quantity` + Angle to account for anisotropic emittance. + + + Returns + ------- + em : `~astropy.units.Quantity` + Emittance. + + """ + + @abstractmethod + def reflectance( + self, + albedo: u.physical.dimensionless, + i: u.physical.angle, + e: u.physical.angle, + phi: u.physical.angle, + ) -> u.Quantity[u.sr**-1]: + r"""Bidirectional reflectance. + + 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 + ---------- + albedo : `~astropy.units.Quantity` + Surface albedo. + + 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 + ------- + bdr : `~astropy.units.Quantity` + Bidirectional reflectance. + + """ + + @staticmethod + 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 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 `absorptance`. + + Input vectors do not need to be normalized. + + + Parameters + ---------- + epsilon : `~astropy.units.Quantity` + Surface emissivity. + + n : `numpy.ndarray` + Surface normal vector. + + 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 + ------- + a : `~astropy.units.Quantity` + Absorptance. + + """ + + i = self._angle(n, r) + return self.absorptance(epsilon, i) + + @u.quantity_input + 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 `emittance`. + + Input vectors do not need to be normalized. + + + Parameters + ---------- + epsilon : `~astropy.units.Quantity` + Surface emissivity. + + 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. + + + Returns + ------- + em : `~astropy.units.Quantity` + Emittance. + + """ + + e = self._angle(n, ro) + phi = self._angle(r, ro) + return self.emittance(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, + ) -> u.Quantity[u.sr**-1]: + """Vector-based alternative to `reflectance`. + + Input vectors do not need to be normalized. + + + Parameters + ---------- + albedo : `~astropy.units.Quantity` + Surface albedo. + + 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. + + + Returns + ------- + bdr : `~astropy.units.Quantity` + Bidirectional reflectance. + + """ + + i = self._angle(n, r) + e = self._angle(n, ro) + phi = self._angle(r, ro) + return self.reflectance(albedo, i, e, phi) diff --git a/sbpy/surfaces/tests/__init__.py b/sbpy/surfaces/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/sbpy/surfaces/tests/test_lambertian.py b/sbpy/surfaces/tests/test_lambertian.py new file mode 100644 index 000000000..a542fd8c3 --- /dev/null +++ b/sbpy/surfaces/tests/test_lambertian.py @@ -0,0 +1,104 @@ +# 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 ..lambertian import LambertianSurface + + +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.absorptance(epsilon, i) + + assert u.allclose(result, expected) + + +def test_absorptance_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.9 * np.sqrt(3) / 2 + + surface = LambertianSurface() + result = surface.absorptance_from_vectors(epsilon, n, r, ro) + + assert u.allclose(result, expected) + + +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.emittance(epsilon, e, None) + assert u.allclose(result, expected) + + +def test_emittance_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.45 + + surface = LambertianSurface() + result = surface.emittance_from_vectors(epsilon, n, r, ro) + assert u.allclose(result, expected) + + # incident vector is optional + result = surface.emittance_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 + expected = ( + 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 = (0.1 * np.sqrt(3) / 2) / 2 / np.pi / u.sr + + surface = LambertianSurface() + result = surface.reflectance_from_vectors(albedo, n, r, ro) + + 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 000000000..7824b4653 --- /dev/null +++ b/sbpy/surfaces/tests/test_scattered.py @@ -0,0 +1,50 @@ +# # Licensed under a 3-clause BSD style license - see LICENSE.rst + +# import pytest +# import numpy as np +# from astropy import units as u + +# from ...calib import Sun, solar_spectrum +# from ..scattered import LambertianSurfaceScatteredSunlight + + +# 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(): +# pytest.importorskip("synphot") + +# 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.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) diff --git a/sbpy/surfaces/tests/test_surface.py b/sbpy/surfaces/tests/test_surface.py new file mode 100644 index 000000000..2b50918d2 --- /dev/null +++ b/sbpy/surfaces/tests/test_surface.py @@ -0,0 +1,37 @@ +# 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 ..surface import Surface, min_zero_cos + + +def test_min_zero_cos(): + a = Angle([-91, -90, 0, 30, 45, 60, 90, 91], "deg") + 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(min_zero_cos(a[i]), expected[i]) + + +@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 new file mode 100644 index 000000000..bdfe4da1b --- /dev/null +++ b/sbpy/units/physical.py @@ -0,0 +1,10 @@ +""" +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 + +u.physical.def_physical_type(u.sr**-1, "bidirectional reflectance") diff --git a/sbpy/units/typing.py b/sbpy/units/typing.py new file mode 100644 index 000000000..e06b39620 --- /dev/null +++ b/sbpy/units/typing.py @@ -0,0 +1,28 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +"""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[u.physical.length], u.Quantity[u.physical.frequency] +] +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], + ]