Merge pull request #208 from mkelley/projected-size-equiv

Projected size equivalencies

docs/sbpy/activity/index.rst

@@ -40,7 +40,7 @@ Apertures may be converted between linear and angular units using :func:`~sbpy.a
 
   >>> ap = sba.CircularAperture(1 * u.arcsec)
   >>> ap.as_length(1 * u.au)  # doctest: +FLOAT_CMP
-  <CircularAperture: radius [725270.94380784] m>
+  <CircularAperture: radius [725.27094381] km>
 
 Ideal comae (constant production rate, free-expansion, infinite lifetime) have *1/ρ* surface brightness distributions.  With :func:`~sbpy.activity.Aperture.coma_equivalent_radius`, we may convert the aperture into a circular aperture that would contain the same total flux at the telescope:
 

docs/sbpy/units.rst

@@ -4,9 +4,7 @@ Units Module (`sbpy.units`)
 Introduction
 ------------
 
-`~sbpy.units` provides common planetary astronomy units and Vega-based
-magnitude conversions.  ``sbpy`` units may be added to the top-level
-``astropy.units`` namespace via:
+`~sbpy.units` provides common planetary astronomy units and unit conversions, including Vega-based magnitudes.  ``sbpy`` units may be added to the top-level ``astropy.units`` namespace via:
 
   >>> import sbpy.units
   >>> sbpy.units.enable()    # doctest: +IGNORE_OUTPUT
@@ -124,6 +122,21 @@ phase angle can be calculated:
   [0.0021763  0.00201223 0.0022041  0.00269637 0.00292785] 1 / sr
 
 
+Projected Sizes
+---------------
+
+With the `~sbpy.units.projected_size` equivalencies, one may convert between angles and lengths, for a given distance:
+
+  >>> import astropy.units as u
+  >>> import sbpy.units as sbu
+  >>> (1 * u.arcsec).to('km', sbu.projected_size(1 * u.au))
+  ... # doctest: +FLOAT_CMP
+  <Quantity [725.27094381] km>
+  >>> (725.27 * u.km).to('arcsec', sbu.projected_size(1 * u.au))
+  ... # doctest: +FLOAT_CMP
+  <Quantity [1.00] arcsec>
+
+
 Reference/API
 -------------
 .. automodapi:: sbpy.units

sbpy/activity/core.py

@@ -10,8 +10,6 @@
 
 
 __all__ = [
-    'rho_as_angle',
-    'rho_as_length',
     'Aperture',
     'CircularAperture',
     'AnnularAperture',
@@ -26,64 +24,7 @@
 import astropy.units as u
 
 from .. import data as sbd
-
-
-def rho_as_angle(rho, eph):
-    """Projected linear distance to angular distance.
-
-
-    Parameters
-    ----------
-    rho : `~astropy.units.Quantity`
-        Projected distance in units of length.
-
-    eph : dictionary-like, `~sbpy.data.Ephem`
-        Ephemerides; requires geocentric distance as `delta`.
-
-
-    Returns
-    -------
-    rho_l : `~astropy.units.Quantity`
-
-    """
-
-    if rho.unit.is_equivalent(u.m):
-        rho_a = np.arctan(rho / eph['delta'].to(u.m))
-    elif rho.unit.is_equivalent(u.rad):
-        rho_a = rho
-    else:
-        raise u.UnitConversionError('rho must have units of length or angle')
-
-    return rho_a
-
-
-def rho_as_length(rho, eph):
-    """Angular distance to projected linear distance.
-
-
-    Parameters
-    ----------
-    rho : `~astropy.units.Quantity`
-        Projected distance in units of angle.
-
-    eph : dictionary-like, `~sbpy.data.Ephem`
-        Ephemerides; requires geocentric distance as `delta`.
-
-
-    Returns
-    -------
-    rho_l : `~astropy.units.Quantity`
-
-    """
-
-    if rho.unit.is_equivalent(u.rad):
-        rho_l = eph['delta'].to(u.m) * np.tan(rho)
-    elif rho.unit.is_equivalent(u.m):
-        rho_l = rho
-    else:
-        raise u.UnitConversionError('rho must have units of length or angle.')
-
-    return rho_l
+from .. import units as sbu
 
 
 class Aperture(ABC):
@@ -131,7 +72,7 @@ def as_angle(self, eph):
 
         """
 
-        dim = rho_as_angle(self.dim, eph)
+        dim = self.dim.to('arcsec', sbu.projected_size(eph))
         return type(self)(dim)
 
     @sbd.dataclass_input(eph=sbd.Ephem)
@@ -152,21 +93,9 @@ def as_length(self, eph):
 
         """
 
-        dim = rho_as_length(self.dim, eph)
+        dim = self.dim.to('km', sbu.projected_size(eph))
         return type(self)(dim)
 
-    def _convert_unit(self, rho, eph):
-        """Make units match those of self."""
-        if not self.dim.unit.is_equivalent(rho.unit):
-            if rho.unit.is_equivalent(u.m):
-                x = rho_as_angle(rho, eph)
-            else:
-                x = rho_as_length(rho, eph)
-        else:
-            x = rho
-
-        return x
-
     @abstractmethod
     def coma_equivalent_radius(self):
         """Circular aperture radius that yields same flux for a 1/ρ coma.
@@ -303,6 +232,7 @@ def coma_equivalent_radius(self):
 class GaussianAperture(Aperture):
     """Gaussian-shaped aperture, e.g., for radio observations.
 
+    The aperture is normalized to 1.0 at the center.
 
     Parameters
     ----------
@@ -343,22 +273,28 @@ def fwhm(self):
         """Beam full-width at half-maximum."""
         return self.dim * 2.3548200450309493
 
-    def __call__(self, rho, eph=None):
+    @sbd.dataclass_input
+    @sbd.quantity_to_dataclass(eph=(sbd.Ephem, 'delta'))
+    def __call__(self, rho, eph: sbd.Ephem=None):
         """Evaluate the aperture.
 
+
         Parameters
         ----------
         rho : `~astropy.units.Quantity`
             Position to evaluate, in units of length or angle.
 
-        eph : dictionary-like, `~sbpy.data.Ephem`, optional
-            Ephemerides at epoch; requires geocentric distance as
-            `delta` keyword if the aperture's units and `rho`'s units
-            do not match.
+        eph : dictionary-like, `~sbpy.data.Ephem`, or `~astropy.units.Quantity`, optional
+            The observer-target distance (``delta``).  Use ``eph`` to
+            convert between angles and lengths, as needed.
 
         """
 
-        x = self._convert_unit(rho, eph)
+        if eph is not None:
+            equiv = sbu.projected_size(eph)
+        else:
+            equiv = []
+        x = rho.to(self.dim.unit, equiv)
 
         # normalize to 1.0 at the center
         return np.exp(-x**2 / self.sigma**2 / 2)

sbpy/activity/dust.py

@@ -26,8 +26,9 @@
 from ..spectroscopy import BlackbodySource
 from .. import data as sbd
 from .. import exceptions as sbe
+from .. import units as sbu
 from ..spectroscopy.sources import SinglePointSpectrumError
-from .core import Aperture, rho_as_length
+from .core import Aperture
 
 
 @bib.cite({
@@ -254,7 +255,7 @@ def to_fluxd(self, wfb, aper, eph, unit=None, **kwargs):
         else:
             rho = aper
             ndim = np.ndim(rho)
-        rho = rho_as_length(rho, eph)
+        rho = rho.to('km', sbu.projected_size(eph))
 
         ndim = max(ndim, np.ndim(self))
 

sbpy/activity/gas/core.py

@@ -30,7 +30,6 @@
 from ...exceptions import RequiredPackageUnavailable
 from .. core import (Aperture, RectangularAperture, GaussianAperture,
                      AnnularAperture, CircularAperture)
-from .. core import rho_as_length
 
 
 def photo_lengthscale(species, source=None):
@@ -293,9 +292,9 @@ def column_density(self, rho, eph=None):
             Projected distance to the region of interest on the plane
             of the sky in units of length or angle.
 
-        eph : dictionary-like, `~sbpy.data.Ephem`, `~astropy.units.Quantity`
-            Target-observer distance, or ephemeris with `delta`.
-            Required if the aperture has angular units.
+        eph : dictionary-like, `~sbpy.data.Ephem`, `~astropy.units.Quantity`. optional
+            Target-observer distance, or ephemeris with ``'delta'``
+            field.  Required to convert rho to a projected size.
 
 
         Returns
@@ -306,8 +305,12 @@ def column_density(self, rho, eph=None):
 
         """
 
-        rho_m = rho_as_length(rho, eph=eph).to('m').value
-        return self._column_density(rho_m) / u.m**2
+        equiv = []
+        if eph is not None:
+            equiv = sbu.projected_size(eph)
+
+        rho = rho.to('m', equiv).value
+        return self._column_density(rho) / u.m**2
 
     @sbd.dataclass_input(eph=sbd.Ephem)
     @sbd.quantity_to_dataclass(eph=(sbd.Ephem, 'delta'))
@@ -333,7 +336,9 @@ def total_number(self, aper, eph=None):
 
         """
 
-        return self._integrate_column_density(aper.as_length(eph))[0]
+        if eph is not None:
+            aper = aper.as_length(eph)
+        return self._integrate_column_density(aper)[0]
 
     @abstractmethod
     def _volume_density(self, r):
@@ -601,7 +606,8 @@ def total_number(self, aper, eph=None):
         if isinstance(aper, u.Quantity):
             aper = CircularAperture(aper)
 
-        aper = aper.as_length(eph)
+        if eph is not None:
+            aper = aper.as_length(eph)
 
         # Inspect aper and handle as appropriate
         if isinstance(aper, (RectangularAperture, GaussianAperture)):

sbpy/activity/tests/test_core.py

@@ -3,46 +3,10 @@
 import pytest
 import numpy as np
 import astropy.units as u
+from ... import units as sbu
 from ..core import *
 
 
-def test_rho_as_angle_length():
-    # arctan(100 km, 1 au) * 206264.806 = 0.13787950659645942
-    rho = rho_as_angle(100 * u.km, {'delta': 1 * u.au})
-    assert np.isclose(rho.to(u.arcsec).value, 0.13787950659645942)
-
-
-def test_rho_as_angle_angle():
-    assert rho_as_angle(1 * u.rad, {'delta': 1 * u.au}).value == 1
-
-
-def test_rho_as_angle_error():
-    with pytest.raises(u.UnitConversionError):
-        rho_as_angle(1 * u.s, {'delta': 1 * u.au})
-
-
-def test_rho_as_length_angle():
-    # 1 au * tan(1") = 725.2709438078363
-    rho = rho_as_length(1 * u.arcsec, {'delta': 1 * u.au})
-    assert np.isclose(rho.to(u.km).value, 725.2709438078363)
-
-
-def test_rho_as_length_length():
-    assert rho_as_length(1 * u.km, {'delta': 1 * u.au}).value == 1
-
-
-def test_rho_as_length_error():
-    with pytest.raises(u.UnitConversionError):
-        rho_as_length(1 * u.s, {'delta': 1 * u.au})
-
-
-def test_rho_roundtrip():
-    a = 10 * u.arcsec
-    eph = {'delta': 1 * u.au}
-    b = rho_as_angle(rho_as_length(a, eph), eph)
-    assert np.isclose(a.value, b.to(u.arcsec).value)
-
-
 class TestCircularAperture:
     def test_init_unit(self):
         with pytest.raises(u.UnitTypeError):
@@ -61,13 +25,15 @@ def test_as_length(self):
         r = 1 * u.arcsec
         aper = CircularAperture(r)
         eph = {'delta': 1 * u.au}
-        assert aper.as_length(eph).dim == rho_as_length(r, eph)
+        length = r.to('km', sbu.projected_size(eph))
+        assert np.isclose(aper.as_length(eph).dim.value, length.value)
 
     def test_as_angle(self):
         r = 100 * u.km
         aper = CircularAperture(r)
         eph = {'delta': 1 * u.au}
-        assert aper.as_angle(eph).dim == rho_as_angle(r, eph)
+        angle = r.to('arcsec', sbu.projected_size(eph))
+        assert np.isclose(aper.as_angle(eph).dim.value, angle.value)
 
 
 class TestAnnularAperture:
@@ -85,18 +51,6 @@ def test_coma_equivalent_radius(self):
         aper = AnnularAperture(shape)
         assert aper.coma_equivalent_radius() == 1 * u.arcsec
 
-    def test_as_length(self):
-        shape = [1, 2] * u.arcsec
-        aper = AnnularAperture(shape)
-        eph = {'delta': 1 * u.au}
-        assert all(aper.as_length(eph).dim == rho_as_length(shape, eph))
-
-    def test_as_angle(self):
-        shape = [100, 200] * u.km
-        aper = AnnularAperture(shape)
-        eph = {'delta': 1 * u.au}
-        assert all(aper.as_angle(eph).dim == rho_as_angle(shape, eph))
-
     def test_shape_error(self):
         with pytest.raises(ValueError):
             AnnularAperture([1, 2, 3] * u.km)
@@ -113,18 +67,6 @@ def test_coma_equivalent_radius(self):
         r = aper.coma_equivalent_radius()
         assert np.isclose(r.value, 0.66776816346357259)
 
-    def test_as_length(self):
-        shape = [1, 2] * u.arcsec
-        aper = RectangularAperture(shape)
-        eph = {'delta': 1 * u.au}
-        assert all(aper.as_length(eph).dim == rho_as_length(shape, eph))
-
-    def test_as_angle(self):
-        shape = [100, 200] * u.km
-        aper = RectangularAperture(shape)
-        eph = {'delta': 1 * u.au}
-        assert all(aper.as_angle(eph).dim == rho_as_angle(shape, eph))
-
     def test_shape_error(self):
         with pytest.raises(ValueError):
             RectangularAperture([1, 2, 3] * u.km)
@@ -149,18 +91,6 @@ def test_coma_equivalent_radius(self):
         aper = GaussianAperture(sigma)
         assert aper.coma_equivalent_radius() == np.sqrt(np.pi / 2) * u.arcsec
 
-    def test_as_length(self):
-        sig = 1 * u.arcsec
-        aper = GaussianAperture(sig)
-        eph = {'delta': 1 * u.au}
-        assert aper.as_length(eph).dim == rho_as_length(sig, eph)
-
-    def test_as_angle(self):
-        sig = 100 * u.km
-        aper = GaussianAperture(sig)
-        eph = {'delta': 1 * u.au}
-        assert aper.as_angle(eph).dim == rho_as_angle(sig, eph)
-
     def test_call_angle(self):
         aper = GaussianAperture(1 * u.arcsec)
         assert np.isclose(aper(1 * u.arcsec), 0.6065306597126334)