Added optika.rulings.HolographicRulingSpacing class to represent ru…

…lings created by interfering two beams. (#77)

README.md

@@ -20,7 +20,7 @@ pip install optika
 ## Features
 - Sequential raytrace modeling
 - Spherical, conical and toroidal surface sag profiles
-- Ruled surfaces, constant and variable line spacing
+- Ruled surfaces, constant, variable, and holographic line spacing
 - Circular, rectangular, and polygonal apertures
 - multilayer reflectivity and transmissivity
 - n-dimensional configurations of the optical system using [named-arrays](https://github.com/sun-data/named-arrays)

docs/refs.bib

@@ -191,5 +191,17 @@ @INPROCEEDINGS{Stern2004
        adsurl = {https://ui.adsabs.harvard.edu/abs/2004SPIE.5171...77S},
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
-
+@article{Welford1975,
+    title = {A vector raytracing equation for hologram lenses of arbitrary shape},
+    journal = {Optics Communications},
+    volume = {14},
+    number = {3},
+    pages = {322-323},
+    year = {1975},
+    issn = {0030-4018},
+    doi = {https://doi.org/10.1016/0030-4018(75)90327-2},
+    url = {https://www.sciencedirect.com/science/article/pii/0030401875903272},
+    author = {W.T. Welford},
+    abstract = {A single vector equation is given, including the effect of wavelength change, for tracing rays in the geometrical optics approximation.}
+}
 

optika/materials/_snells_law.py

@@ -334,7 +334,11 @@ def snells_law(
     if np.any(m != 0):
         d = spacing_rulings
         g = normal_rulings
-        a = a - np.sign(-a @ normal) * (m * wavelength_new * g) / (n1 * d)
+        a = np.where(
+            condition=np.isfinite(d),
+            x=a - np.sign(-a @ normal) * (m * wavelength_new * g) / (n1 * d),
+            y=a,
+        )
 
     c = -a @ normal
 

optika/rulings/__init__.py

@@ -7,6 +7,7 @@
     AbstractRulingSpacing,
     ConstantRulingSpacing,
     Polynomial1dRulingSpacing,
+    HolographicRulingSpacing,
 )
 from ._rulings import (
     AbstractRulings,
@@ -23,6 +24,7 @@
     "AbstractRulingSpacing",
     "ConstantRulingSpacing",
     "Polynomial1dRulingSpacing",
+    "HolographicRulingSpacing",
     "AbstractRulings",
     "Rulings",
     "MeasuredRulings",

optika/rulings/_rulings.py

@@ -294,7 +294,7 @@ def efficiency(
         light, and :math:`\theta` is the angle of incidence inside the medium.
         """
 
-        normal_rulings = self.spacing_(rays.position).normalized
+        normal_rulings = self.spacing_(rays.position, normal).normalized
 
         parallel_rulings = normal.cross(normal_rulings).normalized
 
@@ -441,7 +441,7 @@ def efficiency(
         light, and :math:`\theta` is the angle of incidence inside the medium.
         """
 
-        normal_rulings = self.spacing_(rays.position).normalized
+        normal_rulings = self.spacing_(rays.position, normal).normalized
 
         parallel_rulings = normal.cross(normal_rulings).normalized
 
@@ -594,7 +594,7 @@ def efficiency(
         light, and :math:`\theta` is the angle of incidence inside the medium.
         """
 
-        normal_rulings = self.spacing_(rays.position).normalized
+        normal_rulings = self.spacing_(rays.position, normal).normalized
 
         parallel_rulings = normal.cross(normal_rulings).normalized
 
@@ -741,7 +741,7 @@ def efficiency(
         light, and :math:`\theta` is the angle of incidence inside the medium.
         """
 
-        normal_rulings = self.spacing_(rays.position).normalized
+        normal_rulings = self.spacing_(rays.position, normal).normalized
 
         parallel_rulings = normal.cross(normal_rulings).normalized
 
@@ -901,7 +901,7 @@ def efficiency(
         light, and :math:`\theta` is the angle of incidence inside the medium.
         """
 
-        normal_rulings = self.spacing_(rays.position).normalized
+        normal_rulings = self.spacing_(rays.position, normal).normalized
 
         parallel_rulings = normal.cross(normal_rulings).normalized
 

optika/rulings/_spacing.py

@@ -7,6 +7,7 @@
     "AbstractRulingSpacing",
     "ConstantRulingSpacing",
     "Polynomial1dRulingSpacing",
+    "HolographicRulingSpacing",
 ]
 
 
@@ -26,6 +27,7 @@ class AbstractRulingSpacing(
     def __call__(
         self,
         position: na.AbstractCartesian3dVectorArray,
+        normal: na.AbstractCartesian3dVectorArray,
     ) -> na.AbstractCartesian3dVectorArray:
         """
         The local ruling vector at the given position
@@ -34,6 +36,8 @@ def __call__(
         ----------
         position
             The location on the grating at which to evaluate the ruling spacing.
+        normal
+            The unit vector perpendicular to the surface at the given position.
         """
 
 
@@ -65,6 +69,7 @@ def transformation(self) -> None:
     def __call__(
         self,
         position: na.AbstractCartesian3dVectorArray,
+        normal: na.AbstractCartesian3dVectorArray,
     ) -> na.Cartesian3dVectorArray:
         return self.constant * self.normal
 
@@ -104,19 +109,208 @@ def shape(self) -> dict[str, int]:
     def __call__(
         self,
         position: na.AbstractCartesian3dVectorArray,
+        normal: na.AbstractCartesian3dVectorArray,
     ) -> na.Cartesian3dVectorArray:
 
         coefficients = self.coefficients
-        normal = self.normal
+        normal_rulings = self.normal
         transformation = self.transformation
 
         if transformation is not None:
             position = transformation(position)
 
-        x = position @ normal
+        x = position @ normal_rulings
 
         result = 0 * u.mm
         for power, coefficient in coefficients.items():
             result = result + coefficient * (x**power)
 
-        return result * normal
+        return result * normal_rulings
+
+
+@dataclasses.dataclass(eq=False, repr=False)
+class HolographicRulingSpacing(
+    AbstractRulingSpacing,
+):
+    r"""
+    Rulings created by interfering two beams.
+
+    Examples
+    --------
+
+    Create some holographic rulings from two source points,
+    launch rays from the first source point and confirm they are refocused
+    onto the second source point.
+
+    .. jupyter-execute::
+
+        import matplotlib.pyplot as plt
+        import astropy.units as u
+        import astropy.visualization
+        import named_arrays as na
+        import optika
+
+        # Define the origins of the two recording beams
+        x1 = na.Cartesian3dVectorArray(5, 0, -10) * u.mm
+        x2 = na.Cartesian3dVectorArray(10, 0, -10) * u.mm
+
+        # Define the wavelength of the recording beams
+        wavelength = 500 * u.nm
+
+        # Define the surface normal
+        normal = na.Cartesian3dVectorArray(0, 0, -1)
+
+        # Define input rays emanating from the origin
+        # of the first recording beam
+        position = na.Cartesian3dVectorArray(
+            x=na.linspace(-5, +5, axis="x", num=5) * u.mm,
+            y=0 * u.mm,
+            z=0 * u.mm,
+        )
+        direction_input = position - x1
+
+        # Initialize the holographic ruling spacing
+        # representation
+        rulings = optika.rulings.HolographicRulingSpacing(
+            x1=x1,
+            x2=x2,
+            wavelength=wavelength,
+        )
+
+        # Evaluate the ruling spacing where
+        # the rays strike the surface
+        d = rulings(position, normal)
+
+        # Compute the new direction of some diffracted
+        # rays
+        direction_output = optika.materials.snells_law(
+            wavelength=wavelength,
+            direction=direction_input.normalized,
+            index_refraction=1,
+            index_refraction_new=1,
+            normal=normal,
+            is_mirror=True,
+            diffraction_order=1,
+            spacing_rulings=d.length,
+            normal_rulings=d.normalized,
+        )
+        direction_output = direction_output * 20 * u.mm
+
+        # Plot the results
+        with astropy.visualization.quantity_support():
+            fig, ax = plt.subplots()
+            na.plt.plot(
+                na.stack([x1, position], axis="t"),
+                components=("z", "x"),
+                axis="t",
+                color="tab:blue",
+            )
+            na.plt.plot(
+                na.stack([position, position + direction_output], axis="t"),
+                components=("z", "x"),
+                axis="t",
+                color="tab:orange",
+            )
+            ax.axvline(0)
+            ax.scatter(x1.z, x1.x)
+            ax.scatter(x2.z, x2.x)
+
+    Notes
+    -----
+
+    From :cite:t:`Welford1975`, the ruling spacing is given by
+
+    .. math::
+
+        \mathbf{d} = \frac{\lambda}{a} \mathbf{q} \times \mathbf{n}
+
+    where :math:`\lambda` is the wavelength of the recording beams,
+    :math:`\mathbf{n}` is a unit vector perpendicular to the surface,
+
+    .. math::
+
+        a \mathbf{q} = \mathbf{n} \times (\pm \mathbf{r}_1 \mp \mathbf{r}_2),
+
+    :math:`\mathbf{r}_1` is a unit vector in the direction of the first
+    recording beam,
+    and :math:`\mathbf{r}_2` is a unit vector in the direction of the second
+    recording beam.
+    If rays are diverging from the origin of the recording beams,
+    the top branch is used, otherwise the bottom branch is used.
+    """
+
+    x1: na.AbstractCartesian3dVectorArray
+    """
+    The origin of the first recording beam in local coordinates.
+    """
+
+    x2: na.AbstractCartesian3dVectorArray
+    """
+    The origin of the second recording beam in local coordinates.
+    """
+
+    wavelength: u.Quantity | na.AbstractScalar
+    """
+    The wavelength of the recording beams.
+    """
+
+    is_diverging_1: bool | na.AbstractScalar = True
+    """
+    A boolean flag indicating if rays are diverging from the origin of the
+    first beam.
+    """
+
+    is_diverging_2: bool | na.AbstractScalar = False
+    """
+    A boolean flag indicating if rays are diverging from the origin of the
+    second beam.
+    """
+
+    transformation: None | na.transformations.AbstractTransformation = None
+    """
+    A transformation from surface coordinates to ruling coordinates.
+    """
+
+    @property
+    def shape(self) -> dict[str, int]:
+        return na.broadcast_shapes(
+            optika.shape(self.x1),
+            optika.shape(self.x2),
+            optika.shape(self.is_diverging_1),
+            optika.shape(self.is_diverging_2),
+        )
+
+    def __call__(
+        self,
+        position: na.AbstractCartesian3dVectorArray,
+        normal: na.AbstractCartesian3dVectorArray,
+    ) -> na.Cartesian3dVectorArray:
+
+        x1 = self.x1
+        x2 = self.x2
+        wavelength = self.wavelength
+        d1 = self.is_diverging_1
+        d2 = self.is_diverging_2
+        n = normal
+
+        d1 = 2 * d1 - 1
+        d2 = 2 * d2 - 1
+
+        r1 = position - x1
+        r2 = position - x2
+
+        r1 = d1 * r1.normalized
+        r2 = d2 * r2.normalized
+
+        dr = r1 - r2
+
+        aq = n.cross(dr)
+
+        a = aq.length
+        q = aq / a
+
+        spacing = wavelength / a
+
+        result = spacing * q.cross(n)
+
+        return result

optika/rulings/_spacing_test.py

@@ -22,12 +22,19 @@ class AbstractTestAbstractRulingSpacing(
             ),
         ],
     )
+    @pytest.mark.parametrize(
+        argnames="normal",
+        argvalues=[
+            na.Cartesian3dVectorArray(0, 0, -1),
+        ],
+    )
     def test__call__(
         self,
         a: optika.rulings.AbstractRulingSpacing,
         position: na.AbstractCartesian3dVectorArray,
+        normal: na.Cartesian3dVectorArray,
     ):
-        result = a(position)
+        result = a(position, normal)
         assert isinstance(result, na.AbstractCartesian3dVectorArray)
         assert np.all(result.length > (0 * u.mm))
 
@@ -65,3 +72,19 @@ class TestPolynomial1dRulingSpacing(
     AbstractTestAbstractRulingSpacing,
 ):
     pass
+
+
+@pytest.mark.parametrize(
+    argnames="a",
+    argvalues=[
+        optika.rulings.HolographicRulingSpacing(
+            x1=na.Cartesian3dVectorArray(0, 1, 2) * u.mm,
+            x2=na.Cartesian3dVectorArray(1, 2, 3) * u.mm,
+            wavelength=500 * u.nm,
+        ),
+    ],
+)
+class TestHolographicRulingSpacing(
+    AbstractTestAbstractRulingSpacing,
+):
+    pass