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/calib.rst b/docs/sbpy/calib.rst index 12ba9d86f..46db1e45e 100644 --- a/docs/sbpy/calib.rst +++ b/docs/sbpy/calib.rst @@ -41,6 +41,8 @@ The names of the built-in sources are stored as an internal array. They can be >>> print(sun) +.. _default-spectra: + Controlling the default spectra ------------------------------- @@ -164,6 +166,7 @@ Get the default solar spectrum, observe it through the Johnson V-band filter (di >>> print(fluxd) # doctest: +FLOAT_CMP -26.744715028702647 mag(JM) +.. _calib-binning-vs-interpolation: Binning versus interpolation with ``observe()`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 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..c4b1a5b8a --- /dev/null +++ b/docs/substitutions.txt @@ -0,0 +1,41 @@ +.. 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` + +.. sbpy +.. DataClass +.. |Ephem| replace:: :class:`~sbpy.data.ephem.Ephem` +.. |Orbit| replace:: :class:`~sbpy.data.orbit.Orbit` +.. |Obs| replace:: :class:`~sbpy.data.obs.Obs` +.. |Phys| replace:: :class:`~sbpy.data.phys.Phys` diff --git a/sbpy/shape/__init__.py b/sbpy/shape/__init__.py index 75634340e..b8bdcaf7e 100644 --- a/sbpy/shape/__init__.py +++ b/sbpy/shape/__init__.py @@ -2,4 +2,5 @@ SBPy Module for small body shapes """ -from .core import * +from .core import Shape +from .sphere import Sphere diff --git a/sbpy/shape/core.py b/sbpy/shape/core.py index 95a8c6be6..cfdff71dd 100644 --- a/sbpy/shape/core.py +++ b/sbpy/shape/core.py @@ -1,72 +1,101 @@ # 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'] +__all__ = ["Shape"] + +import astropy.units as u + +from ..calib import Sun +from ..data.ephem import Ephem +from ..spectroscopy.sources import SinglePointSpectrumError + + +class Shape: + """Model asteroid or comet nucleus shape.""" + + @staticmethod + def _incident_sunlight( + wfb, + rh: u.Quantity, + unit: u.Unit, + interpolate: bool, + ) -> u.Quantity: + """Calculate incident sunlight for a normal surface.""" + + sun = Sun.from_default() + try: + S = sun.observe(wfb, unit=unit, interpolate=interpolate) + except SinglePointSpectrumError: + S = sun(wfb, unit=unit) + S /= rh.to_value("au") ** 2 + + return S + + +# __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/shape/sphere.py b/sbpy/shape/sphere.py new file mode 100644 index 000000000..f6e7b08a6 --- /dev/null +++ b/sbpy/shape/sphere.py @@ -0,0 +1,355 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +Asteroid and comet nucleus shape models +""" + +__all__ = ["Sphere"] + +from typing import Callable + +import numpy as np +from numpy import pi +import astropy.units as u +from scipy.integrate import dblquad + +from ..data.ephem import Ephem +from ..data.decorators import dataclass_input +from ..surfaces.surface import Surface +from .core import Shape +from .transformations import twovec + + +class Sphere(Shape): + """A spherical object. + + + Parameters + ---------- + radius : |Quantity| + The radius of the sphere. + + pole : tuple of |Quantity|, optional + The direction of the pole as longitude and latitude. + + """ + + @u.quantity_input + @dataclass_input + def __init__( + self, + radius: u.Quantity["length"], + pole: tuple[u.Quantity["angle"]] | None = None, + ): + self.radius: u.Quantity = radius + self.pole = None if pole is None else (pole[0], pole[1]) + + @dataclass_input + def integrate( + self, + func: Callable, + **kwargs, + ) -> tuple[u.Quantity, u.Quantity]: + """Integrate arbitrary function over the surface of the sphere. + + + Parameters + ---------- + func : callable + The function to integrate: ``func(phi, theta, *args)``. + + kwargs : dict, optional + Keyword arguments passed to the integrator, + `~scipy.integrate.dblquad`. + + + Returns + ------- + total : |Quantity| + The result. + + err : |Quantity| + Estimated integration error on ``total``. + + """ + + result, err = dblquad(func, 0, pi, 0, 2 * pi, **kwargs) * self.radius**2 + + return result, err + + def integrate_i_e_phi( + self, + func: Callable, + phase: u.physical.angle, + **kwargs, + ) -> tuple[u.Quantity, u.Quantity]: + """Integrate a function of $(i, e, \phi)$ over the sphere. + + + Parameters + ---------- + func : callable + The function to integrate: ``func(i, e, phase, *args)``, + where :math:`i` is the angle of incidence, :math:`e` is the the + angle of emittance, and :math:`phase` is the phase angle. + + phase : |Quantity| + Sun-target-observer (phase) angle. + + kwargs : dict, optional + Additional keyword arguments are passed to the integrator, + `~scipy.integrate.dblquad`. + + + Returns + ------- + total : |Quantity| + The result. + + err : |Quantity| + Estimated integration error on ``total``. + + """ + + # Let target-observer vector be +z + r_obs = np.array([0, 0, 1]) + + # Define target-Sun vector using the phase angle + r_sun = np.array([np.sin(phase), 0, np.cos(phase)]) + + # Define a matrix that transforms vectors in a reference frame with + # target-Sun along +z + if np.isclose(phase, 0): + # if r_obs == r_sun, then use the identity matrix + M = np.eye(3) + else: + M = twovec(r_sun, 2, r_obs, 0) + + def f(phi, theta, *args): + """ + z-axis is to the observer + theta is the angle from the z-axis + phi is the angle from the x-axis + """ + + r = np.array( + [ + np.cos(phi) * np.sin(theta), + np.sin(phi) * np.sin(theta), + np.cos(theta), + ] + ) + + # angle of emittance + e = theta * u.rad + + # transform to Sun-centered frame to find angle of incidence + rprime = np.dot(M, r) + i = np.arctan2(np.hypot(rprime[0], rprime[1]), rprime[2]) * u.rad + + a = func(i, e, phase, *args) * np.sin(theta) + return a + + return self.integrate(f, **kwargs) + + @u.quantity_input + @dataclass_input + def absorption( + self, + I: u.Quantity, + epsilon: u.Quantity["dimensionless"], + surface: Surface, + **kwargs, + ) -> tuple[u.Quantity, u.Quantity]: + """Surface total absorbed light. + + + Parameters + ---------- + I : |Quantity| + Incident light. + + epsilon : |Quantity| or float + Surface emissivity. + + surface : `~sbpy.surface.surface.Surface` + The surface being illuminated. + + kwargs : dict, optional + Additional keyword arguments are passed to the integrator, + `~scipy.integrate.dblquad`. + + + Returns + ------- + A : |Quantity| + Surface total absorbed energy. + + err : |Quantity| + Estimated surface integration error on ``A``. + + """ + + def f(i, e, phi): + return surface.absorptance(1, i) + + # set phase to zero: absorption does not depend on an observer + A, err = self.integrate_i_e_phi(f, 0 * u.deg, **kwargs) + return (I * epsilon * A), (I * epsilon * err) + + @dataclass_input + def absorption_of_sunlight( + self, + wfb, + epsilon: u.Quantity["dimensionless"] | float, + surface: Surface, + eph: Ephem, + interpolate: bool = False, + unit: u.Unit = "W / um", + **kwargs, + ) -> tuple[u.Quantity, u.Quantity]: + """Surface total absorbed sunlight. + + Uses the current default ``Sun``. See :ref:`_sbpy-calib` and + :ref:`_default-spectra` for more information. + + + Parameters + ---------- + wfb : |Quantity|, |SpectralElement|, str, or list + Wavelengths, frequencies, or bandpasses. + + epsilon : |Quantity| or float + Surface emissivity. + + surface : `~sbpy.surface.surface.Surface` + The surface being illuminated. + + eph : |Ephem| + Distance to the sun. + + interpolate : bool, optional + Set to ``True`` to interpolate the solar spectrum for ``wfb``. See + :ref:`_calib-binning-vs-interpolation` for guidance. + + unit : |Unit| or str, optional + Spectral density units for the return value (e.g., W/μm). + + kwargs : dict, optional + Additional keyword arguments are passed to the integrator, + `~scipy.integrate.dblquad`. + + + Returns + ------- + A : |Quantity| + Surface total absorbed sunlight at each ``wfb``. + + err : |Quantity| + Estimated surface integration error on ``A``. + + """ + + unit = u.Unit(unit) + S = Shape._incident_sunlight(wfb, eph["rh"][0], unit / u.m**2, interpolate) + A, err = self.absorption(S, epsilon, surface, **kwargs) + return A.to(unit), err.to(unit) + + @dataclass_input + def reflected_light( + self, + I: u.Quantity, + albedo: u.Quantity["dimensionless"], + surface: Surface, + eph: Ephem, + **kwargs, + ) -> tuple[u.Quantity, u.Quantity]: + """Surface total reflected light. + + + Parameters + ---------- + I : |Quantity| + Incident light. + + albedo : |Quantity| or float + Surface albedo. + + surface : `~sbpy.surface.surface.Surface` + The surface being illuminated. + + eph : |Ephem| + Target-observer distance and Sun-target-observer (phase) angle. + + kwargs : dict, optional + Additional keyword arguments are passed to the integrator, + `~scipy.integrate.dblquad`. + + + Returns + ------- + R : |Quantity| + Reflected light at observer. + + err : |Quantity| + Estimated surface integration error on ``R``. + + """ + + delta = eph["delta"][0] + phase = eph["phase"][0] + + def f(i, e, phi): + return surface.reflectance(1, i, e, phi) + + R, err = self.integrate_i_e_phi(f, phase, **kwargs) + return ( + (I * albedo * R / delta**2).to(I.unit), + (I * albedo * err / delta**2).to(I.unit), + ) + + @dataclass_input + def reflected_sunlight( + self, + wfb, + albedo: u.Quantity["dimensionless"], + surface: Surface, + eph: Ephem, + unit: u.Unit | str = "W / (m2 um)", + interpolate: bool = False, + **kwargs, + ) -> tuple[u.Quantity, u.Quantity]: + """Surface total reflected sunlight. + + Uses the current default ``Sun``. See :ref:`_sbpy-calib` and + :ref:`_default-spectra` for more information. + + + Parameters + ---------- + wfb : |Quantity|, |SpectralElement|, str, or list + Wavelengths, frequencies, or bandpasses. + + albedo : |Quantity| or float + Surface albedo. + + surface : `~sbpy.surface.surface.Surface` + The surface being illuminated. + + eph : |Ephem| + Heliocentric distance, target-observer distance and + Sun-target-observer (phase) angle. + + interpolate : bool, optional + Set to ``True`` to interpolate the solar spectrum for ``wfb``. See + :ref:`_calib-binning-vs-interpolation` for guidance. + + kwargs : dict, optional + Additional keyword arguments are passed to the integrator, + `~scipy.integrate.dblquad`. + + """ + + unit = u.Unit(unit) + S = Shape._incident_sunlight(wfb, eph["rh"][0], unit, interpolate) + R, err = self.reflected_light(S, albedo, surface, eph, **kwargs) + return R.to(unit), err.to(unit) diff --git a/sbpy/shape/tests/__init__.py b/sbpy/shape/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/sbpy/shape/tests/test_core.py b/sbpy/shape/tests/test_core.py new file mode 100644 index 000000000..40d0a50e7 --- /dev/null +++ b/sbpy/shape/tests/test_core.py @@ -0,0 +1,20 @@ +# 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 ..core import Shape + + +def test_Shape_incident_sunlight(): + # Test branching on wfb input, scaling by rh + wave = np.linspace(0.4, 0.9) * u.um + unit = u.Jy + rh = 2 * u.au + with solar_spectrum.set(Sun.from_array(wave, 1000 * unit)): + S = Shape._incident_sunlight(wave, rh, unit, True) + assert u.allclose(S, 250 * u.Jy) + + S = Shape._incident_sunlight(550 * u.nm, rh, unit, False) + assert u.isclose(S, 250 * u.Jy, atol=0.02 * unit) diff --git a/sbpy/shape/tests/test_sphere.py b/sbpy/shape/tests/test_sphere.py new file mode 100644 index 000000000..c9911de8e --- /dev/null +++ b/sbpy/shape/tests/test_sphere.py @@ -0,0 +1,84 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import numpy as np +from numpy import pi +from astropy import units as u + +from ..sphere import Sphere +from ...calib import Sun, solar_spectrum +from ...surfaces.surface import Surface, min_zero_cos + + +class SimpleSurface(Surface): + def absorptance(self, epsilon, i): + return epsilon * min_zero_cos(i) + + def emittance(self, albedo, e, phi): + pass # pragma: no cover + + def reflectance(self, albedo, i, e, phi): + return albedo * min_zero_cos(i) * min_zero_cos(e) + + +class TestSphere: + def test_integrate_i_e_phi(self): + s = Sphere(1 * u.km) + + def func(i, e, phase): + return 1 + + # test phase = 0 special case + phase = 0 * u.deg + area, _ = s.integrate_i_e_phi(func, phase) + assert u.isclose(area, 4 * pi * u.km**2) + + # test phase > 0 + phase = 30 * u.deg + area, _ = s.integrate_i_e_phi(func, phase) + assert u.isclose(area, 4 * pi * u.km**2) + + # test hemisphere + def func(i, e, phase): + return 1.0 * (e.to_value(u.rad) < pi / 2) + + area, _ = s.integrate_i_e_phi(func, phase) + assert u.isclose(area, 2 * pi * u.km**2) + + # test cos(i < 90) + def func(i, e, phase): + return min_zero_cos(i) + + area, _ = s.integrate_i_e_phi(func, 0 * u.deg) + assert u.isclose(area, pi * u.km**2) + + def test_absorption_of_sunlight(self): + surface = SimpleSurface() + s = Sphere(1 * u.km) + eph = dict(rh=2 * u.au, delta=2 * u.au, phase=0 * u.deg) + wave = np.logspace(-0.5, 0.5) * u.um + + sun = Sun.from_array(wave, 1000 * u.W / u.m**2 / u.um) + + with solar_spectrum.set(sun): + a, _ = s.absorption_of_sunlight(wave, 0.95, surface, eph, interpolate=True) + + expected = 950 / 4 * sun.fluxd.unit * pi * s.radius**2 + assert u.allclose(a, expected) + + def test_reflected_sunlight(self): + surface = SimpleSurface() + s = Sphere(1 * u.km) + eph = dict(rh=2 * u.au, delta=2 * u.au, phase=0 * u.deg) + wave = np.logspace(-0.5, 0.5) * u.um + + sun = Sun.from_array(wave, 1000 * u.W / u.m**2 / u.um) + + with solar_spectrum.set(sun): + r, _ = s.reflected_sunlight(wave, 0.1, surface, eph, interpolate=True) + + # (%i2) integrate (2 * %pi * cos(x)**2 * sin(x), x, 0, %pi / 2); + # 2 %pi + # (%o2) ----- + # 3 + expected = 100 / 4 * sun.fluxd.unit * 2 / 3 * pi * s.radius**2 / (2 * u.au) ** 2 + assert u.allclose(r, expected) diff --git a/sbpy/shape/tests/test_transformations.py b/sbpy/shape/tests/test_transformations.py new file mode 100644 index 000000000..8c77ec18a --- /dev/null +++ b/sbpy/shape/tests/test_transformations.py @@ -0,0 +1,57 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import pytest +import numpy as np +from astropy import units as u + +from ..transformations import twovec, xyz2sph, sph2xyz + + +def test_twovec(): + x = [1, 0, 0] * u.au + y = [0, 1, 0] * u.au + z = [0, 0, 1] * u.au + + # x -> y' + # y -> z' + # z -> x' + M = twovec(z, 0, y, 2) + assert u.allclose(np.dot(M, x), y) + assert u.allclose(np.dot(M, y), z) + assert u.allclose(np.dot(M, z), x) + + # x -> y' + # y -> -x' + # z -> z' + M = twovec(x, 1, z, 2) + assert u.allclose(np.dot(M, x), y) + assert u.allclose(np.dot(M, y), -x) + assert u.allclose(np.dot(M, z), z) + + with pytest.raises(ValueError): + M = twovec(z, 0, 2 * z, 1) + + +def test_xyz2sph(): + x = [1, 0, 0] * u.au + y = [0, 1, 0] * u.au + z = [0, 0, 1] * u.au + + assert u.allclose(xyz2sph(*x), [0, 0] * u.deg) + assert u.allclose(xyz2sph(*y), [90, 0] * u.deg) + assert u.allclose(xyz2sph(*z), [0, 90] * u.deg) + + +def test_sph2xyz(): + x = [1, 0, 0] * u.au + y = [0, 1, 0] * u.au + z = [0, 0, 1] * u.au + + assert u.allclose(sph2xyz(0 * u.deg, 0 * u.deg), x.value) + assert u.allclose(sph2xyz(90 * u.deg, 0 * u.deg), y.value, atol=1e-16) + assert u.allclose(sph2xyz(0 * u.deg, 90 * u.deg), z.value, atol=1e-16) + + r = 1 * u.au + assert u.allclose(sph2xyz(0 * u.deg, 0 * u.deg, r), x) + assert u.allclose(sph2xyz(90 * u.deg, 0 * u.deg, r), y, atol=1e-16 * u.au) + assert u.allclose(sph2xyz(0 * u.deg, 90 * u.deg, r), z, atol=1e-16 * u.au) diff --git a/sbpy/shape/transformations.py b/sbpy/shape/transformations.py new file mode 100644 index 000000000..58ee68627 --- /dev/null +++ b/sbpy/shape/transformations.py @@ -0,0 +1,118 @@ +import numpy as np +import astropy.units as u + + +def twovec( + axdef: np.ndarray, indexa: int, plndef: np.ndarray, indexp: int +) -> np.ndarray: + """Transformation matrix to a new coordinate defined by two input vectors. + + + Parameters + ---------- + axdef : array-like float containing 3 elements + The vector (x, y, z) that defines one of the principal axes of the new + coordinate frame. + + indexa : int 0, 1, or 2 + Specify which of the three coordinate axes is defined by ``axdef``. 0 + for x-axis, 1 for y-axis, and 2 for z-axis + + plndef : array-like float containing 3 elements + The vector (x, y, z) that defines (with ``axdef``) a principal plane of + the new coordinate frame. + + indexp : int 0, 1, or 2 + Specify the second axis of the principal frame determined by ``axdef`` + and ``plndef`` + + + Returns + ------- + M : numpy array of shape 3x3 + The transformation matrix that convert a vector from the old coordinate + to the coordinate frame defined by the input vectors via a dot product. + + + Raises + ------ + ValueError + When ``axdef`` and ``plndef`` are linearly dependent. + + + Notes + ----- + This routine is directly translated form SPICE lib routine twovec.f (cf. + `SPICE manual + `_). + + The indexing of array elements are different in FORTRAN (that SPICE is + originally based) from Python. Here 0-based index is used. + + Note that the original twovec.f in SPICE toolkit returns matrix that + converts a vector in the new frame to the original frame, opposite to what + is implemented here. + + """ + + axdef = np.asarray(axdef).flatten() + plndef = np.asarray(plndef).flatten() + + if np.linalg.norm(np.cross(axdef, plndef)) == 0: + raise ValueError( + "The input vectors AXDEF and PLNDEF are linearly" + " correlated and can't define a coordinate frame." + ) + + M = np.eye(3) + i1 = indexa % 3 + i2 = (indexa + 1) % 3 + i3 = (indexa + 2) % 3 + + M[i1, :] = axdef / np.linalg.norm(axdef) + if indexp == i2: + xv = np.cross(axdef, plndef) + M[i3, :] = xv / np.linalg.norm(xv) + xv = np.cross(xv, axdef) + M[i2, :] = xv / np.linalg.norm(xv) + else: + xv = np.cross(plndef, axdef) + M[i2, :] = xv / np.linalg.norm(xv) + xv = np.cross(axdef, xv) + M[i3, :] = xv / np.linalg.norm(xv) + + return M + + +def xyz2sph( + x: np.ndarray | u.Quantity, + y: np.ndarray | u.Quantity, + z: np.ndarray | u.Quantity, +) -> np.ndarray | u.Quantity: + """Convert (x, y, z) to (lon, lat).""" + x_ = np.asanyarray(x) + y_ = np.asanyarray(y) + z_ = np.asanyarray(z) + lon = np.arctan2(y_, x_) + complete_angle = ( + u.Quantity(2 * np.pi, u.rad) if isinstance(lon, u.Quantity) else 2 * np.pi + ) + lon = (lon + complete_angle) % complete_angle + lat = np.arctan2(z_, np.sqrt(x_ * x_ + y_ * y_)) + return np.stack([lon, lat]) + + +def sph2xyz( + lon: np.ndarray | u.Quantity, + lat: np.ndarray | u.Quantity, + r: float | u.Quantity | None = None, +) -> np.ndarray | u.Quantity: + """Convert (lon, lat) to (x, y, z), with a default length of unity.""" + if r is None: + r = 1.0 * u.dimensionless_unscaled if isinstance(lon, u.Quantity) else 1.0 + lon_ = np.asanyarray(lon) + lat_ = np.asanyarray(lat) + x = r * np.cos(lat_) * np.cos(lon_) + y = r * np.cos(lat_) * np.sin(lon_) + z = r * np.sin(lat_) + return np.stack([x, y, z]) 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..a024fb3de --- /dev/null +++ b/sbpy/surfaces/surface.py @@ -0,0 +1,244 @@ +# 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 + + +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/surfaces/tests/test_thermal.py b/sbpy/surfaces/tests/test_thermal.py new file mode 100644 index 000000000..189bbd209 --- /dev/null +++ b/sbpy/surfaces/tests/test_thermal.py @@ -0,0 +1,69 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import pytest + +import numpy as np +from numpy import pi +from astropy.coordinates import Angle +from astropy import units as u +import astropy.constants as const +from astropy.modeling.models import BlackBody + +from ..thermal import InstantaneousEquilibrium, ThermalEmission +from ..lambertian import LambertianSurface + + +class ConstantTemperature(ThermalEmission): + def __init__(self, T0): + super().__init__(0 * u.W / u.m**2, 0, 1, 1, LambertianSurface()) + self._T0 = T0 + + @property + def T0(self): + return self._T0 + + def T(self, i): + return self._T0 + + +class TestThermalEmission: + def test_emission(self): + bb200 = ConstantTemperature(200 * u.K) + wave = 10 * u.um + L = bb200.emission(wave, 0 * u.deg, 0 * u.deg, 0 * u.deg) + L200 = BlackBody(200 * u.K)(wave).to(L.unit, u.spectral_density(wave)) + assert u.isclose(L, L200) + + L306060 = bb200.emission(wave, 30 * u.deg, 60 * u.deg, 60 * u.deg) + assert u.isclose(L306060, L200 / 2) + + bb0 = ConstantTemperature(0 * u.K) + wave = 10 * u.um + L = bb0.emission(wave, 0 * u.deg, 0 * u.deg, 0 * u.deg) + assert u.isclose(L, 0 * L.unit) + + +class TestInstantaneousEquilibrium: + def test_T0(self): + S = 1361.16646541 * u.W / u.m**2 + thermal = InstantaneousEquilibrium(S, 0, 1, 1, LambertianSurface()) + T = (S.value / (pi * 5.6703744191844314e-08)) ** (1 / 4) * u.K + + assert u.isclose(thermal.T0, T) + + # only albedo and eta affect temperature (emissivity is for emission) + thermal.albedo = 0.1 + T = ((1 - 0.1) * S.value / (pi * 5.6703744191844314e-08)) ** (1 / 4) * u.K + assert u.isclose(thermal.T0, T) + + thermal.eta = 1.5 + T = ((1 - 0.1) * S.value / (pi * 1.5 * 5.6703744191844314e-08)) ** (1 / 4) * u.K + assert u.isclose(thermal.T0, T) + + def test_temperature(self): + S = 1361.16646541 * u.W / u.m**2 + thermal = InstantaneousEquilibrium(S, 0, 1, 1, LambertianSurface()) + assert thermal.T(0 * u.deg) == thermal.T0 + + T60 = (S.value / (2 * pi * 5.6703744191844314e-08)) ** (1 / 4) * u.K + assert u.isclose(thermal.T(60 * u.deg), T60) diff --git a/sbpy/surfaces/thermal.py b/sbpy/surfaces/thermal.py new file mode 100644 index 000000000..324285873 --- /dev/null +++ b/sbpy/surfaces/thermal.py @@ -0,0 +1,118 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import abc +import numpy as np +from numpy import pi +import astropy.units as u +import astropy.constants as const +from astropy.modeling.models import BlackBody +from .surface import Surface + + +class ThermalEmission(abc.ABC): + """Abstract base class for calculating surface thermal emission.""" + + def __init__(self, I, albedo, epsilon, eta, surface: Surface): + self.I = I + self.albedo = albedo + self.epsilon = epsilon + self.eta = eta + self.surface: Surface = surface + self._update_T0 = True + + @property + def I(self) -> u.Quantity["W/m2"]: + return self._I + + @I.setter + @u.quantity_input + def I(self, i: u.Quantity["W/m2"]): + self._I = i + self._update_T0 = True + + @property + def albedo(self) -> u.Quantity[""]: + return self._albedo + + @albedo.setter + @u.quantity_input + def albedo(self, a: u.Quantity[""]): + self._albedo = a + self._update_T0 = True + + @property + def epsilon(self) -> u.Quantity[""]: + return self._epsilon + + @epsilon.setter + @u.quantity_input + def epsilon(self, e: u.Quantity[""]): + self._epsilon = e + self._update_T0 = True + + @property + def eta(self) -> u.Quantity[""]: + return self._eta + + @eta.setter + @u.quantity_input + def eta(self, e: u.Quantity[""]): + self._eta = e + self._update_T0 = True + + @property + def surface(self) -> Surface: + return self._surface + + @surface.setter + @u.quantity_input + def surface(self, s: Surface): + self._surface = s + self._update_T0 = True + + @abc.abstractmethod + def T0(self) -> u.Quantity[u.K]: + """Temperature of a surface normal to the incident radiation.""" + pass + + @property + @abc.abstractmethod + def T(self, i: u.Quantity["angle"]) -> u.Quantity[u.K]: + """Temperature of a surface at angle ``i`` with incident radiation.""" + pass + + def emission(self, wave_freq, i, e, phi, unit=None): + """Observed spectral radiance of a surface element.""" + + if unit is None: + unit = u.W / u.m**2 / u.sr / wave_freq.unit + + T = self.T(i) + + if u.isclose(T, 0 * u.K): + return 0 * unit + + B = BlackBody(T)(wave_freq) + # emissivity = 1; it was already accounted for in the temperature + I_e = B * self.surface.emittance(1, e, phi) + + return I_e.to(unit, u.spectral_density(wave_freq)) + + +class InstantaneousEquilibrium(ThermalEmission): + @property + def T0(self) -> u.Quantity[u.K]: + """Temperature of a surface normal to the incident radiation.""" + if self._update_T0: + I_a = self.I * self.surface.absorptance(1 - self.albedo, 0 * u.deg) + self._T0 = ( + (I_a / (pi * self.epsilon * self.eta * const.sigma_sb)) ** (1 / 4) + ).decompose() + self._update_T0 = False + + return self._T0 + + @u.quantity_input + def T(self, i: u.Quantity["angle"]) -> u.Quantity[u.K]: + # emissivity = 1; it was already accounted when T0 was calculated + return self.T0 * self.surface.absorptance(1, i) ** (1 / 4) diff --git a/sbpy/thermal/core.py b/sbpy/thermal/core.py index 1807b75d6..d2223911a 100644 --- a/sbpy/thermal/core.py +++ b/sbpy/thermal/core.py @@ -5,70 +5,253 @@ created on June 27, 2017 """ -__all__ = ['ThermalClass', 'STM', 'FRM', 'NEATM'] +__all__ = ["ThermalClass", "NonRotThermalModel", "FastRotThermalModel"] -__doctest_requires__ = { - "ThermalClass.flux": "astroquery" -} +import abc +import numpy as np +from numpy import pi +from numpy.linalg import norm +from scipy.integrate import dblquad +import astropy.units as u +import astropy.constants as const +from astropy.modeling.models import BlackBody +from ..data import Phys, Obs, Ephem, dataclass_input, quantity_to_dataclass +from ..shape import Sphere +from ..surfaces import LambertianSurface +from ..surfaces.thermal import InstantaneousEquilibrium -class ThermalClass(): +def neatm(wave_freq, diameter, albedo, epsilon, eta, eph, unit="W m-2 um-1"): + """Near-Earth Asteroid Thermal Model.""" - def flux(phys, eph, lam): - """Model flux density for a given wavelength `lam`, or a list/array thereof + shape = Sphere(diameter / 2) + surface = LambertianSurface() + S = const.L_sun / (4 * pi * eph["rh"][0] ** 2) + thermal = InstantaneousEquilibrium(S, albedo, epsilon, eta, surface) - Parameters - ---------- - phys : `sbpy.data.Phys` instance, mandatory - provide physical properties - eph : `sbpy.data.Ephem` instance, mandatory - provide object ephemerides - lam : `astropy.units` quantity or list-like, mandatory - wavelength or list thereof + def f(i, e, phi): + return thermal.emission(wave_freq, i, e, phi, unit=unit * u.m**2) - Examples - -------- - >>> from astropy.time import Time - >>> from astropy import units as u - >>> from sbpy.thermal import STM - >>> from sbpy.data import Ephem, Phys - >>> epoch = Time('2019-03-12 12:30:00', scale='utc') - >>> eph = Ephem.from_horizons('2015 HW', location='568', epochs=epoch) # doctest: +REMOTE_DATA - >>> phys = PhysProp('diam'=0.3*u.km, 'pv'=0.3) # doctest: +SKIP - >>> lam = np.arange(1, 20, 5)*u.micron # doctest: +SKIP - >>> flux = STM.flux(phys, eph, lam) # doctest: +SKIP + fluxd = shape.integrate_i_e_phi(f, eph["phase"][0]) / eph["delta"][0] ** 2 + return fluxd.to(unit) - not yet implemented - """ +# @u.quantity_input(rh=u.km, R=u.km) +# def __init__(self, rh, R, albedo=0.1, emissivity=1.0, beaming=1.0): +# """ +# Parameters +# ---------- +# rh : u.Quantity +# Heliocentric distance +# R : u.Quantity +# Radius of asteroid +# albedo : float, u.Quantity +# Bolometric Bond albedo +# emissivity : float, u.Quantity +# Emissivity of surface +# beaming : float, u.Quantity +# Beaming parameter +# """ +# self.rh = rh +# self.R = R +# self.albedo = albedo +# self.emissivity = emissivity +# # self.beaming = beaming +# class ThermalEmission(abc.ABC): +# """Thermal emission from an asteroid or comet nucleus.""" - def fit(self, eph): - """Fit thermal model to observations stored in `sbpy.data.Ephem` instance +# def __init__(self, shape: Shape, surface: Surface): +# self.shape: Shape = shape +# self.surface: Surface = surface - Parameters - ---------- - eph : `sbpy.data.Ephem` instance, mandatory - provide object ephemerides and flux measurements +# @abc.abstractmethod +# def T( +# self, lon: u.Quantity["angle"], lat: u.Quantity["angle"] +# ) -> u.Quantity["temperature"]: +# """Temperature on the surface of an object. - Examples - -------- - >>> from sbpy.thermal import STM - >>> stmfit = STM.fit(eph) # doctest: +SKIP - not yet implemented +# Parameters +# ---------- +# lon : |Quantity| +# Longitude - """ +# lat : |Quantity| +# Latitude +# """ -class STM(ThermalClass): - pass +# # Needs to be overridden in subclasses. This function needs to be able +# # to return a valid quantity for the full range of lon and lat, i.e., +# # include the night side of an object. +# pass -class FRM(ThermalClass): - pass +# # @u.quantity_input(wave_freq=u.m, equivalencies=u.spectral()) +# # def _int_func(self, lon, lat, m, unit, wave_freq): +# # """Integral function for `fluxd`. +# # Parameters +# # ---------- +# # lon : float +# # Longitude in radiance +# # lat : float +# # Latitude in fradiance +# # m : numpy array of shape (3, 3) +# # Transformation matrix to convert a vector in the frame to perform +# # integration to body-fixed frame. This matrix can be calculated +# # with private method `_transfer_to_bodyframe`. The integration +# # is performed in a frame where the sub-observer point is defined at +# # lon = 0, lat = 0. +# # unit : str or astropy.units.Unit +# # Unit of the integral function. +# # wave_freq : u.Quantity +# # Wavelength or frequency of calculation -class NEATM(ThermalClass): - def __init__(self): - from .. import bib - bib.register('sbpy.thermal.NEATM', {'method': '1998Icar..131..291H'}) +# # Returns +# # ------- +# # float : Integral function to calculate total flux density. +# # """ +# # _, lon1, lat1 = xyz2sph(m.dot(sph2xyz(lon, lat))) +# # T = self.T(lon1 * u.rad, lat1 * u.rad) +# # if np.isclose(T, 0 * u.K): +# # return 0.0 +# # else: +# # # the integral term needs to include a correction for latitude +# # # with cos(lat), and a Lambertian emission term cos(lat) + cos(lon) +# # coslat = np.cos(lat) +# # coslon = np.cos(lon) +# # f = BlackBody(T)(wave_freq) * coslat * coslat * coslon +# # return f.to_value(unit, u.spectral_density(wave_freq)) + +# # @staticmethod +# # @u.quantity_input(sublon=u.deg, sublat=u.deg) +# # def _transfer_to_bodyframe(sublon, sublat): +# # """Calculate transformation matrix. + +# # The numerical integration to calculate total flux density is performed +# # in a reference frame where the sub-observer point is at +# # (lon, lat) = (0, 0). This matrix supports the transformation from +# # this frame to the body-fixed frame to facilitate the calculation of +# # surface temperature. +# # """ +# # coslat = np.cos(sublat).value +# # if abs(coslat) < np.finfo(type(coslat)).resolution: +# # if sublat.value > 0: +# # m = np.array([[0, 0, -1], [0, 1, 0], [1, 0, 0]]) +# # else: +# # m = np.array([[0, 0, 1], [0, 1, 0], [-1, 0, 0]]) +# # else: +# # m = twovec( +# # [sublon.to_value("deg"), sublat.to_value("deg")], 0, [0, 90], 2 +# # ).T +# # return m + +# @u.quantity_input( +# wave_freq=u.m, delta=u.m, lon=u.deg, lat=u.deg, equivalencies=u.spectral() +# ) +# def fluxd( +# self, +# wave_freq: u.Quantity["length"] | u.Quantity["frequency"], +# eph: Ephem, +# unit="W m-2 um-1", +# **kwargs +# ): +# """Model thermal flux density of an object. + +# Parameters +# ---------- +# eph : `sbpy.data.Ephem` instance, mandatory +# provide object ephemerides and flux measurements + +# Examples +# -------- +# >>> from sbpy.thermal import STM +# >>> stmfit = STM.fit(eph) # doctest: +SKIP + +# not yet implemented + +# """ + + +# class NonRotThermalModel(ThermalEmission): +# """Non-rotating / instantaneous temperature distribution.""" + +# # @property +# # def Tss(self): +# # f_sun = const.L_sun / (4 * np.pi * self.rh**2) +# # return ( +# # ( +# # (1 - self.albedo) +# # * f_sun +# # / (self.beaming * self.emissivity * const.sigma_sb) +# # ) +# # ** 0.25 +# # ).decompose() + +# u.quantity_input + +# def T( +# self, lon: u.Quantity["angle"], lat: u.Quantity["angle"] +# ) -> u.Quantity["temperature"]: +# """Temperature on the surface of an object. + + +# Parameters +# ---------- +# lon : |Quantity| +# Longitude + +# lat : |Quantity| +# Latitude + +# """ + +# # @u.quantity_input(lon=u.deg, lat=u.deg) +# # def T(self, lon, lat): +# # """Surface temperature at specific (lat, lon) + +# # lon : u.Quantity in units equivalent to deg +# # Longitude +# # lat : u.Quantity in units equivalent to deg +# # Latitude + +# # Returns +# # ------- +# # u.Quantity : Surface temperature. +# # """ +# # coslon = np.cos(lon) +# # coslat = np.cos(lat) +# # prec = np.finfo(coslat.value).resolution +# # if (abs(coslon) < prec) or (abs(coslat) < prec) or (coslon < 0): +# # return 0 * u.K +# # else: +# # return self.Tss * (coslon * coslat) ** 0.25 + + +# class FastRotThermalModel(ThermalClass): +# """Fast-rotating object temperature distribution, i.e., FRM""" + +# @property +# def Tss(self): +# f_sun = const.L_sun / (4 * np.pi * self.rh**2) +# return ( +# ((1 - self.albedo) * f_sun / (np.pi * self.emissivity * const.sigma_sb)) +# ** 0.25 +# ).decompose() + +# @u.quantity_input(lon=u.deg, lat=u.deg) +# def T(self, lon, lat): +# """Surface temperature at specific (lat, lon) + +# lon : u.Quantity in units equivalent to deg +# Longitude +# lat : u.Quantity in units equivalent to deg +# Latitude + +# Returns +# ------- +# u.Quantity : Surface temperature. +# """ +# coslat = np.cos(lat) +# return self.Tss * coslat**0.25 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], + ]