Surfaces and shapes

[mkelley]

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.sr

Surface 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 T

At 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 sphere

Solar 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"]**2

Asteroid 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.

[mkelley]

If this looks good, I would proceed in a few steps:

• [ ] Surface, LambertianSurface, SurfaceReflectance
• [ ] SurfaceThermalEmission, InstantaneousThermalEquilibrium, RapidRotatorThermalEquilibrium
• [ ] Shape and Sphere
• [ ] SolarSystemObjectModel, NEATM, RapidRotatorThermalModel

[jianyangli]

I think this is a reasonable approach. My feeling is that there might be some redundancy in the architecture. For example, I still need to make sense of the combination of Shape, SolarSystemObjectModel, and Sphere, and maybe the organization can be simplified. But I don't have a good idea how yet. Let's move on this way and I'll keep thinking and get back if I have something.

[mkelley]

Thanks. The concept behind the SolarSystemObjectModel (I'm open to better names) is to require the total_fluxd method. This is a method that requires a surface and a shape. I didn't think it made sense to add total_fluxd to any of the other classes, hence SSOM was invented.