Added optika.direction() and optika.angles() functions. (#62)

optika/__init__.py

@@ -3,7 +3,7 @@
 """
 
 from . import mixins
-from ._util import shape
+from ._util import shape, direction, angles
 from . import vectors
 from . import targets
 from . import rays
@@ -21,6 +21,8 @@
 __all__ = [
     "mixins",
     "shape",
+    "direction",
+    "angles",
     "vectors",
     "targets",
     "rays",

optika/_tests/test_systems.py

@@ -68,20 +68,36 @@ def test_rayfunction_stops(self, a: optika.systems.AbstractSequentialSystem):
     def test_field_min(self, a: optika.systems.AbstractSequentialSystem):
         result = a.field_min
         assert isinstance(result, na.AbstractCartesian2dVectorArray)
+        if a.object_is_at_infinity:
+            assert na.unit(result).is_equivalent(u.deg)
+        else:
+            assert na.unit(result).is_equivalent(u.m)
 
     def test_field_max(self, a: optika.systems.AbstractSequentialSystem):
         result = a.field_max
         assert isinstance(result, na.AbstractCartesian2dVectorArray)
         assert np.all(result > a.field_min)
+        if a.object_is_at_infinity:
+            assert na.unit(result).is_equivalent(u.deg)
+        else:
+            assert na.unit(result).is_equivalent(u.m)
 
     def test_pupil_min(self, a: optika.systems.AbstractSequentialSystem):
         result = a.pupil_min
         assert isinstance(result, na.AbstractCartesian2dVectorArray)
+        if a.object_is_at_infinity:
+            assert na.unit(result).is_equivalent(u.m)
+        else:
+            assert na.unit(result).is_equivalent(u.deg)
 
     def test_pupil_max(self, a: optika.systems.AbstractSequentialSystem):
         result = a.pupil_max
         assert isinstance(result, na.AbstractCartesian2dVectorArray)
         assert np.all(result > a.pupil_min)
+        if a.object_is_at_infinity:
+            assert na.unit(result).is_equivalent(u.m)
+        else:
+            assert na.unit(result).is_equivalent(u.deg)
 
     @pytest.mark.parametrize(
         argnames="wavelength,field,pupil",

optika/_util.py

@@ -1,9 +1,13 @@
 from typing import Any
+import numpy as np
+import astropy.units as u
 import named_arrays as na
 import optika
 
 __all__ = [
     "shape",
+    "direction",
+    "angles",
 ]
 
 
@@ -32,3 +36,62 @@ def shape(a: Any) -> dict[str, int]:
         return a.shape
     else:
         return dict()
+
+
+def direction(
+    angles: na.AbstractCartesian2dVectorArray,
+) -> na.Cartesian3dVectorArray:
+    r"""
+    Given a 2D vector of azimuth and elevation angles, convert to a
+    3D vector of direction cosines.
+
+    Parameters
+    ----------
+    angles
+        A vector of azimuth and elevation angles.
+
+    Notes
+    -----
+    If the azimuth and elevation angles are taken to :math:`\phi_x` and
+    :math:`\phi_y`, the direction cosines :math:`\vec{d}` can be found by
+
+    .. math::
+
+        \vec{d} = R_y(\phi_x) R_x(\phi_y) \hat{z}
+
+    where :math:`R_x(\theta)` and :math:`R_y(\theta)` are the rotation matrices
+    about the :math:`x` and :math:`y` axes respectively.
+
+    See Also
+    --------
+    :func:`angles` : Inverse of this function
+    """
+    return na.Cartesian3dVectorArray(
+        x=-np.cos(angles.y) * np.sin(angles.x),
+        y=-np.sin(angles.y),
+        z=+np.cos(angles.y) * np.cos(angles.x),
+    )
+
+
+def angles(
+    direction: na.AbstractCartesian3dVectorArray,
+) -> na.Cartesian2dVectorArray:
+    """
+    Convert a 3D vector of direction cosines to a 2D vector of azimuth and
+    elevation angles.
+
+    Parameters
+    ----------
+    direction
+        A vector of direction cosines.
+
+    See Also
+    --------
+    :func:`direction` : Inverse of this function
+    """
+    if na.unit(direction) is None:
+        direction = direction << u.dimensionless_unscaled
+    return na.Cartesian2dVectorArray(
+        x=-np.arctan2(direction.x, direction.z).to(u.deg),
+        y=-np.arcsin(direction.y / direction.length).to(u.deg),
+    )

optika/_util_test.py

@@ -0,0 +1,48 @@
+import pytest
+import numpy as np
+import astropy.units as u
+import named_arrays as na
+import optika
+
+
+@pytest.mark.parametrize(
+    argnames="angles",
+    argvalues=[
+        na.Cartesian2dVectorArray(1, 2) * u.deg,
+        na.Cartesian2dVectorLinearSpace(
+            start=0,
+            stop=1,
+            axis=na.Cartesian2dVectorArray("x", "y"),
+            num=5,
+        ),
+    ],
+)
+def test_direction(angles: na.AbstractCartesian2dVectorArray):
+    result = optika.direction(angles)
+
+    rotation_x = na.Cartesian3dYRotationMatrixArray(-angles.x)
+    rotation_y = na.Cartesian3dXRotationMatrixArray(+angles.y)
+    result_expected = rotation_x @ rotation_y @ na.Cartesian3dVectorArray(z=1)
+    # result_expected = rotation_y @ rotation_x @ na.Cartesian3dVectorArray(z=1)
+
+    print(f"{result=}")
+    print(f"{result_expected=}")
+
+    assert isinstance(result, na.AbstractCartesian3dVectorArray)
+    assert np.allclose(result, result_expected)
+
+
+@pytest.mark.parametrize(
+    argnames="direction",
+    argvalues=[
+        na.Cartesian3dVectorArray(1, 2, 5).normalized,
+    ],
+)
+def test_angles(direction: na.AbstractCartesian3dVectorArray):
+    result = optika.angles(direction)
+
+    print(f"{direction=}")
+    print(f"{optika.direction(result)=}")
+
+    assert isinstance(result, na.AbstractCartesian2dVectorArray)
+    assert np.allclose(direction, optika.direction(result))

optika/systems.py

@@ -401,7 +401,8 @@ def field_min(self) -> na.AbstractCartesian2dVectorArray:
         """
         axis = (self._axis_field_stop, self._axis_pupil_stop)
         if self.object_is_at_infinity:
-            return self.rayfunction_stops.outputs.direction.xy.min(axis)
+            angles = optika.angles(self.rayfunction_stops.outputs.direction)
+            return angles.min(axis)
         else:
             return self.rayfunction_stops.outputs.position.xy.min(axis)
 
@@ -412,7 +413,8 @@ def field_max(self) -> na.AbstractCartesian2dVectorArray:
         """
         axis = (self._axis_field_stop, self._axis_pupil_stop)
         if self.object_is_at_infinity:
-            return self.rayfunction_stops.outputs.direction.xy.max(axis)
+            angles = optika.angles(self.rayfunction_stops.outputs.direction)
+            return angles.max(axis)
         else:
             return self.rayfunction_stops.outputs.position.xy.max(axis)
 
@@ -426,7 +428,8 @@ def pupil_min(self) -> na.AbstractCartesian2dVectorArray:
         if self.object_is_at_infinity:
             return self.rayfunction_stops.outputs.position.xy.min(axis)
         else:
-            return self.rayfunction_stops.outputs.direction.xy.min(axis)
+            angles = optika.angles(self.rayfunction_stops.outputs.direction)
+            return angles.min(axis)
 
     @property
     def pupil_max(self):
@@ -438,7 +441,8 @@ def pupil_max(self):
         if self.object_is_at_infinity:
             return self.rayfunction_stops.outputs.position.xy.max(axis)
         else:
-            return self.rayfunction_stops.outputs.direction.xy.max(axis)
+            angles = optika.angles(self.rayfunction_stops.outputs.direction)
+            return angles.max(axis)
 
     def _calc_rayfunction_input(
         self,