Added `named_arrays.AbstractCartesian3dVectorArray.solid_angle_cell()…

…` method. (#60)

docs/conf.py

@@ -36,6 +36,7 @@
     'sphinx.ext.autosummary',
     'sphinx.ext.intersphinx',
     'sphinx.ext.inheritance_diagram',
+    'sphinxcontrib.bibtex',
     'jupyter_sphinx',
     'sphinx_favicon'
 ]
@@ -45,11 +46,9 @@
 # autoclass_content = 'both'
 autodoc_typehints = "description"
 
-# autosummary_filename_map = {
-#     'kgpy.optics.Surface': 'kgpy.optics.Surface_cls',
-# }
-
-# typehints_fully_qualified = True
+suppress_warnings = [
+    'autosummary.import_cycle',
+]
 
 graphviz_output_format = 'png'
 inheritance_graph_attrs = dict(rankdir='TB')
@@ -106,6 +105,10 @@
 # https://github.com/readthedocs/readthedocs.org/issues/2569
 master_doc = 'index'
 
+bibtex_bibfiles = ['refs.bib']
+bibtex_default_style = 'plain'
+bibtex_reference_style = 'author_year'
+
 intersphinx_mapping = {
     'python': ('https://docs.python.org/3', None),
     'numpy': ('https://numpy.org/doc/stable/', None),

docs/index.rst

@@ -11,6 +11,13 @@ named_arrays is an implementation of named tensors that includes support for lin
     named_arrays
 
 
+References
+==========
+
+.. bibliography::
+
+|
+
 
 Indices and tables
 ==================

docs/refs.bib

@@ -0,0 +1,15 @@
+@article{Eriksson1990,
+    author = {Folke Eriksson},
+    title = {On the Measure of Solid Angles},
+    journal = {Mathematics Magazine},
+    volume = {63},
+    number = {3},
+    pages = {184--187},
+    year = {1990},
+    publisher = {Taylor \& Francis},
+    doi = {10.1080/0025570X.1990.11977515},
+    URL = {https://doi.org/10.1080/0025570X.1990.11977515},
+    eprint = {https://doi.org/10.1080/0025570X.1990.11977515}
+}
+
+

named_arrays/_vectors/cartesian/tests/test_vectors_cartesian_3d.py

@@ -134,6 +134,21 @@ class AbstractTestAbstractCartesian3dVectorArray(
     def test_xy(self, array: na.AbstractCartesian3dVectorArray):
         assert isinstance(array.xy, na.Cartesian2dVectorArray)
 
+    def test_solid_area_cell(
+        self,
+        array: na.AbstractCartesian3dVectorArray,
+    ):
+        for c in array.components:
+            component = na.as_named_array(array.components[c])
+            if not isinstance(component, na.AbstractScalar):
+                return
+        if array.ndim != 2:
+            return
+        result = array.solid_angle_cell()
+        assert isinstance(result, na.AbstractScalar)
+        assert np.all(result >= 0)
+        assert result.unit.is_equivalent(u.sr)
+
     @pytest.mark.parametrize("array_2", _cartesian3d_arrays_2())
     def test_cross(
             self,

named_arrays/_vectors/cartesian/vectors_cartesian_3d.py

@@ -3,6 +3,8 @@
 from typing_extensions import Self
 import abc
 import dataclasses
+import numpy as np
+import astropy.units as u
 import named_arrays as na
 
 __all__ = [
@@ -108,6 +110,94 @@ def volume_cell(
 
         return result
 
+    @classmethod
+    def _sold_angle(
+        cls,
+        a: na.AbstractCartesian3dVectorArray,
+        b: na.AbstractCartesian3dVectorArray,
+        c: na.AbstractCartesian3dVectorArray,
+    ) -> na.ScalarLike:
+
+        numerator = a @ b.cross(c)
+
+        a_ = a.length
+        b_ = b.length
+        c_ = c.length
+
+        d0 = a_ * b_ * c_
+        d1 = (a @ b) * c_
+        d2 = (a @ c) * b_
+        d3 = (b @ c) * a_
+        denomerator = d0 + d1 + d2 + d3
+
+        unit = numerator.unit
+
+        if unit is not None:
+            numerator = numerator.to(unit).value
+            denomerator = denomerator.to(unit).value
+
+        angle = 2 * np.arctan2(numerator, denomerator)
+
+        return angle << u.sr
+
+    def solid_angle_cell(
+        self,
+        axis: None | tuple[str, str] = None,
+    ) -> na.AbstractScalar:
+        r"""
+        Compute the solid angle of each cell formed by interpreting this
+        array as a logically-rectangular 2D grid of vertices.
+
+        Note that this method is usually only used for sorted arrays
+
+        Parameters
+        ----------
+        axis
+            The two axes defining the logically-rectangular 2D grid.
+            If :obj:`None` (the default), :attr:`axes` is used and must have
+            only two elements.
+
+        Notes
+        -----
+        The solid angle :math:`\Omega` of a triangle formed by the vertices
+        :math:`\vec{a}`, :math:`\vec{b}`, and :math:`\vec{c}` is given by
+        :cite:t:`Eriksson1990` as
+
+        .. math::
+
+            \tan \left( \frac{1}{2} \Omega \right)
+                = \frac{\vec{a} \cdot (\vec{b} \times \vec{c})}
+                    {a b c + (\vec{a} \cdot \vec{b}) c + (\vec{a} \cdot \vec{c}) b + (\vec{b} \cdot \vec{c}) a}.
+
+        Each rectangular cell is decomposed into two triangles and then the
+        solid angle of each triangle is computed.
+        """
+
+        shape = self.shape
+
+        if axis is None:
+            axis = tuple(shape)
+        else:
+            if not set(axis).issubset(shape):   # pragma: nocover
+                raise ValueError(
+                    f"{axis=} should be a subset of {self.shape=}."
+                )
+
+        ax, ay = axis
+
+        s0 = slice(None, ~0)
+        s1 = slice(+1, None)
+
+        a = self[{ax: s0, ay: s0}]
+        b = self[{ax: s1, ay: s0}]
+        c = self[{ax: s1, ay: s1}]
+        d = self[{ax: s0, ay: s1}]
+
+        angle_1 = self._sold_angle(a, b, c)
+        angle_2 = self._sold_angle(c, d, a)
+
+        return angle_1 + angle_2
+
     def cross(
             self,
             other: AbstractCartesian3dVectorArray,

pyproject.toml

@@ -33,6 +33,7 @@ doc = [
     "pytest",
     "graphviz",
     "sphinx-autodoc-typehints",
+    "sphinxcontrib-bibtex",
     "pydata-sphinx-theme",
     "ipykernel",
     "jupyter-sphinx",