Surfaces and shapes #415
Description
This issue is for designing new classes around the concepts of surfaces and shapes, with the goal of being able to model surface reflectance, thermal emission, and observations of arbitrary shape models.
I may have some misunderstanding of the concepts and terminology, so comments and improvements are welcome! Also, the code below is incomplete, but there should be enough sketched out to see how it would work.
Surface
The Surface class accounts for all surface properties, e.g., albedo, emissivity, and how light is reflected off the surface. It has four methods. Three methods are used to describe reflectance, absorbance, and emittance. A fourth calculates the surface brightness (spectral radiance):
class Surface:
"""Abstract base class for all surface characteristics.
Parameters
----------
phys : Phys
Surface physical parameters, e.g., albedo.
"""
@abstractmethod
def reflectance(self, i, e, phi):
r"""Reflectance.
The surface is illuminated by incident flux density, :math:`F_i`, at an
angle of :math:`i`, and emitted toward an angle of :math:`e`, measured
from the surface normal direction. :math:`\phi` is the
Sun-target-observer (phase) angle.
"""
@abstractmethod
def absorptance(self, i) -> float:
r"""Absorption of incident light.
The surface is illuminated by incident flux density, :math:`F_i`, at an
angle of :math:`i`, measured from the surface normal direction.
"""
@abstractmethod
def emittance(self, e, phi):
r"""Emission of light.
The surface is observed at an angle of :math:`e`, measured from the
surface normal direction, and at a solar phase angle of :math:`phi`.
"""
@abstractmethod
def radiance(self, F_i, n, rs, ro):
"""Observed radiance from a surface.
Parameters
----------
F_i : `astropy.units.Quantity`
Incident light, spectral flux density.
n : `numpy.ndarray`
Surface normal vector.
rs : `astropy.units.Quantity`
Radial vector from the surface to the light source.
ro : `astropy.units.Quantity`
Radial vector from the surface to the observer.
"""A Lambertian surface could be implemented like so:
class LambertianSurface(Surface):
"""Lambertian surface."""
def reflectance(self, i, e, phi):
return self.phys["albedo"] * np.cos(i) * self.emittance(e, phi)
def absorptance(self, i):
return (1 - self.phys["albedo"]) * np.cos(i)
def emittance(self, e, phi)
# independent of phase angle
return np.cos(e)Surface reflectance or thermal emission classes would implement the radiance() method. Since we are a solar system focused project, we assume the surface is illuminated by sunlight:
class SurfaceReflectance(Surface):
"""Reflectance from a surface illuminated by sunlight."""
def radiance(self, wave_freq, n, rs, ro):
"""Surface radiance.
Parameters
----------
wave_freq: `astropy.units.Quantity`
Wavelength or frequency.
n : `numpy.ndarray`
Surface normal vector.
rs : `astropy.units.Quantity`
Radial vector from the surface to the light source.
ro : `astropy.units.Quantity`
Radial vector from the surface to the observer.
"""
sun = Sun.from_default()
F_i = sun.observe(wave_freq)
n_hat = np.linalg.norm(n.value)
rs_hat = np.linalg.norm(rs.value)
ro_hat = np.linalg.norm(ro.value)
i = np.arccos(np.dot(n_hat, rs_hat))
e = np.arccos(np.dot(n_hat, ro_hat))
phi = np.arccos(np.dot(rs_hat, ro_hat))
return F_i * self.reflectance(i, e, phi) / u.srSurface thermal emission needs a temperature to produce a blackbody spectrum. The temperature method will be defined by sub-classes. Below we implement an instantaneous LTE class:
class SurfaceThermalEmission(Surface):
"""Thermal emission from a surface illuminated by sunlight.
Parameters
----------
phys : Phys
Surface physical parameters:
- albedo
- emissivity
- beaming parameter
"""
@abstractmethod
def T(self, i, rh):
"""Surface temperature given incidence angle and heliocentric distance."""
def radiance(self, wave_freq, n, rs, ro):
n_hat = np.linalg.norm(n.value)
rs_hat = np.linalg.norm(rs.value)
ro_hat = np.linalg.norm(ro.value)
i = np.arccos(np.dot(n_hat, rs_hat))
e = np.arccos(np.dot(n_hat, ro_hat))
phi = np.arccos(np.dot(rs_hat, ro_hat))
rh = np.sqrt(np.dot(rs, rs))
T = self.T(i, rh)
epsilon = self.phys["emissivity"]
# Note the phase dependency of emittance doesn't make sense here. Should thermal models override the emittance method and drop phi?
return epsilon * BlackBody1D(temperature=T) * self.emittance(e, phi)
class InstantaneousThermalEquilibrium(SurfaceThermalEquilibrium):
"""A surface in instantaneous thermal equilibrium with absorbed sunlight."""
def T(self, i, rh):
sun = const.L / (4 * pi * rh**2)
epsilon = self.phys["emissivity"]
eta = self.phys["eta"]
T = (self.absorptance(i) * sun / epsilon / eta / const.sigma_sb) ** (1 / 4)
return TAt this point concrete classes can be made and used:
class ModelReflectance(SurfaceReflectance, LambertianSurface):
"""Simple surface scattering."""
class ModelThermalEmission(InstantaneousThermalEquilibrium, LambertianSurface):
"""Simple thermal emission."""Shape
Shape classes would account for the scattering and observing geometry for an object. Arbitrary shapes should be possible, but I think just three classes could be provided by sbpy: Sphere, Ellipsoid, and TriangularMesh.
Since the problem at hand is focused on observed brightness of an object, I think we just need one method would integrate a function (e.g., radiance) over the observed surface:
class Shape:
@abstractmethod
def integrate_over_surface(self, eph, func):
"""Integrate a function over the observed surface.
Parameters
----------
eph : `Ephem`
The observing geometry as an ephemeris object:
- rh
- delta
- phase
func : callable
The function to integrate.
"""Derived classes would take into account the object's actual geometry. For example, a Sphere class would run a double integral in spherical coordinates:
class Sphere(Shape):
def integrate_over_surface(self, func, eph):
# integrate func over the observed sphereSolar system object models
Before we can implement an asteroid thermal model, we need a mechanism that combines the concepts of surfaces and shapes. Specifically, we want something that can integrate the surface radiance over the object.
class SolarSystemObjectModel(Surface, Shape):
class total_fluxd(self, wave_freq, eph):
"""Total observed spectral flux density.
Parameters
----------
wave_freq: `astropy.units.Quantity`
Wavelength or frequency.
"""
def func(*args):
self.radiance(wave_freq, *args)
return self.integrate_over_surface(eph, func) / eph["delta"]**2Asteroid thermal models
Now we have everything we need to implement a model like the Near-Earth Asteroid Thermal Model. The NEATM assumes a Lambertian surface illuminated by the Sun in instantaneous thermal equilibrium:
class NEATM(Sphere, InstantaneousThermalEquilibrium, LambertianSurface, SolarSystemObjectModel):
"""The Near-Earth Asteroid Thermal Model of Harris (1998)."""A chart that illustrates the inheritance:
Similarly, a rapid rotator model could be derived using a RapidRotatorThermalEquilibrium class (not shown):
class RapidRotatorThermalModel(Sphere, RapidRotatorThermalEquilibrium, LambertianSurface, SolarSystemObjectModel):
"""Rapid (fast) rotator model."""We also want to implement a Standard Thermal Model. Since it doesn't explicitly use the shape model, it may be simpler to implement as a stand-alone class.