Added optika.chemical.Chemical.n()

optika/chemicals/_chemicals.py

@@ -81,14 +81,21 @@ def file_nk(self) -> na.ScalarArray:
 
         return result
 
-    @functools.cached_property
-    def _wavelength_n_k(self) -> tuple[na.AbstractScalar, ...]:
+    def n(
+        self,
+        wavelength: u.Quantity | na.AbstractScalar,
+    ) -> na.AbstractScalar:
         """
-        Return arrays of wavelength, :math:`n`, and :math:`k` from :attr:`file_nk`.
+        The complex index of refaction of this chemical for a given wavelength
         """
         file_nk = self.file_nk
         shape_base = file_nk.shape
 
+        shape = na.broadcast_shapes(shape_base, wavelength.shape)
+
+        wavelength = na.broadcast_to(wavelength, shape)
+        result = np.empty_like(wavelength.value, dtype=complex)
+
         for i, index in enumerate(file_nk.ndindex()):
             file_nk_index = file_nk[index].ndarray
 
@@ -100,9 +107,9 @@ def _wavelength_n_k(self) -> tuple[na.AbstractScalar, ...]:
                     else:
                         break
 
-            wavelength_index, n_index, k_index = np.genfromtxt(
+            wavelength_index, n_index, k_index = np.loadtxt(
                 fname=file_nk_index,
-                skip_header=skip_header,
+                skiprows=skip_header,
                 unpack=True,
             )
 
@@ -111,35 +118,31 @@ def _wavelength_n_k(self) -> tuple[na.AbstractScalar, ...]:
             n_index = na.ScalarArray(n_index, axes="wavelength")
             k_index = na.ScalarArray(k_index, axes="wavelength")
 
-            shape_index = na.shape_broadcasted(wavelength_index, n_index, k_index)
-            shape = na.broadcast_shapes(shape_base, shape_index)
-
-            if i == 0:
-                wavelength = na.ScalarArray.empty(shape) << u.AA
-                n = na.ScalarArray.empty(shape)
-                k = na.ScalarArray.empty(shape)
-
-            wavelength[index] = wavelength_index
-            n[index] = n_index
-            k[index] = k_index
+            result[index] = na.interp(
+                x=wavelength,
+                xp=wavelength_index,
+                fp=n_index + 1j * k_index
+            )
 
-        return wavelength, n, k
+        return result
 
-    @property
-    def index_refraction(self) -> na.FunctionArray[na.ScalarArray, na.ScalarArray]:
+    def index_refraction(
+        self,
+        wavelength: u.Quantity | na.AbstractScalar,
+    ) -> na.AbstractScalar:
         """
-        The index of refraction of this material as a function of wavelength.
+        The index of refraction of this chemical for the given wavelength.
         """
-        wavelength, n, k = self._wavelength_n_k
-        return na.FunctionArray(inputs=wavelength, outputs=n)
+        return np.real(self.n(wavelength))
 
-    @property
-    def wavenumber(self) -> na.FunctionArray[na.ScalarArray, na.ScalarArray]:
+    def wavenumber(
+        self,
+        wavelength: u.Quantity | na.AbstractScalar,
+    ) -> na.AbstractScalar:
         """
-        The wavenumber of this material as a function of wavelength.
+        The wavenumber of this chemical for the given wavelength.
         """
-        wavelength, n, k = self._wavelength_n_k
-        return na.FunctionArray(inputs=wavelength, outputs=k)
+        return np.imag(self.n(wavelength))
 
 
 @dataclasses.dataclass(eq=False, repr=False)
@@ -159,35 +162,38 @@ class Chemical(
     .. jupyter-execute::
 
         import matplotlib.pyplot as plt
+        import astropy.units as u
         import named_arrays as na
         import optika
 
         si = optika.chemicals.Chemical("Si")
         sio2 = optika.chemicals.Chemical("SiO2")
 
-        n_si = si.index_refraction
-        n_sio2 = sio2.index_refraction
+        wavelength = na.geomspace(10, 10000, axis="wavelength", num=1001) * u.AA
+
+        n_si = si.index_refraction(wavelength)
+        n_sio2 = sio2.index_refraction(wavelength)
 
         fig, ax = plt.subplots(constrained_layout=True)
-        na.plt.plot(n_si.inputs, n_si.outputs, label="silicon");
-        na.plt.plot(n_sio2.inputs, n_sio2.outputs, label="silicon dioxide");
+        na.plt.plot(wavelength, n_si, label="silicon");
+        na.plt.plot(wavelength, n_sio2, label="silicon dioxide");
         ax.set_xscale("log");
-        ax.set_xlabel(f"wavelength ({n_si.inputs.unit:latex_inline})");
+        ax.set_xlabel(f"wavelength ({wavelength.unit:latex_inline})");
         ax.set_ylabel("index of refraction");
         ax.legend();
 
     Plot the wavenumbers of silicon and silicon dioxide
 
     .. jupyter-execute::
 
-        k_si = si.wavenumber
-        k_sio2 = sio2.wavenumber
+        k_si = si.wavenumber(wavelength)
+        k_sio2 = sio2.wavenumber(wavelength)
 
         fig, ax = plt.subplots(constrained_layout=True)
-        na.plt.plot(k_si.inputs, k_si.outputs, label="silicon");
-        na.plt.plot(k_sio2.inputs, k_sio2.outputs, label="silicon dioxide");
+        na.plt.plot(wavelength, k_si, label="silicon");
+        na.plt.plot(wavelength, k_sio2, label="silicon dioxide");
         ax.set_xscale("log");
-        ax.set_xlabel(f"wavelength ({k_si.inputs.unit:latex_inline})");
+        ax.set_xlabel(f"wavelength ({wavelength.unit:latex_inline})");
         ax.set_ylabel("wavenumber");
         ax.legend();
     """

optika/chemicals/_tests/test_chemicals.py

@@ -2,10 +2,18 @@
 import abc
 import pathlib
 import numpy as np
+import astropy.units as u
 import named_arrays as na
 import optika
 
 
+_wavelength = [
+    500 * u.nm,
+    na.geomspace(10, 10000, axis="wavelength", num=101) * u.AA,
+    na.NormalUncertainScalarArray(500 * u.nm, width=10 * u.nm),
+]
+
+
 class AbstractTestAbstractChemical(
     abc.ABC,
 ):
@@ -28,13 +36,25 @@ def test_file_nk(self, a: optika.chemicals.AbstractChemical):
             assert isinstance(result[index].ndarray, pathlib.Path)
             assert result[index].ndarray.exists()
 
-    def test_index_refraction(self, a: optika.chemicals.AbstractChemical):
-        result = a.index_refraction
-        assert isinstance(result, na.AbstractFunctionArray)
+    @pytest.mark.parametrize("wavelength", _wavelength)
+    def test_index_refraction(
+        self,
+        a: optika.chemicals.AbstractChemical,
+        wavelength: u.Quantity | na.AbstractScalar,
+    ):
+        result = a.index_refraction(wavelength)
+        assert isinstance(result, na.AbstractScalar)
+        assert np.all(result > 0)
 
-    def test_wavenumber(self, a: optika.chemicals.AbstractChemical):
-        result = a.wavenumber
-        assert isinstance(result, na.AbstractFunctionArray)
+    @pytest.mark.parametrize("wavelength", _wavelength)
+    def test_wavenumber(
+        self,
+        a: optika.chemicals.AbstractChemical,
+        wavelength: u.Quantity | na.AbstractScalar,
+    ):
+        result = a.wavenumber(wavelength)
+        assert isinstance(result, na.AbstractScalar)
+        assert np.all(result >= 0)
 
 
 @pytest.mark.parametrize(

optika/materials/_multilayers.py

@@ -195,17 +195,7 @@ def multilayer_efficiency(
 
         # Compute the complex index of refraction for the silicon substrate
         silicon = optika.chemicals.Chemical("Si")
-        index_refraction_substrate = na.interp(
-            x=wavelength,
-            xp=silicon.index_refraction.inputs,
-            fp=silicon.index_refraction.outputs,
-        )
-        wavenumber_substrate = na.interp(
-            x=wavelength,
-            xp=silicon.wavenumber.inputs,
-            fp=silicon.wavenumber.outputs,
-        )
-        n_substrate = index_refraction_substrate + wavenumber_substrate * 1j
+        n_substrate = silicon.n(wavelength)
 
         # Compute the reflectivity and transmissivity of this multilayer stack
         reflectivity, transmissivity = optika.materials.multilayer_efficiency(
@@ -409,23 +399,7 @@ def multilayer_efficiency(
                     formula=formula_j,
                 )
 
-                index_refaction_j = chemical_j.index_refraction
-                index_refaction_j = na.interp(
-                    x=wavelength,
-                    xp=index_refaction_j.inputs,
-                    fp=index_refaction_j.outputs,
-                    axis="wavelength",
-                )
-
-                wavenumber_j = chemical_j.wavenumber
-                wavenumber_j = na.interp(
-                    x=wavelength,
-                    xp=wavenumber_j.inputs,
-                    fp=wavenumber_j.outputs,
-                    axis="wavelength",
-                )
-
-                n_j = index_refaction_j + wavenumber_j * 1j
+                n_j = chemical_j.n(wavelength)
                 n_cache[formula_j] = n_j
 
         direction_j = snells_law(

optika/materials/_tests/test_multilayers.py

@@ -160,17 +160,7 @@ def test_multilayer_transmissivity_vs_file(
     transmissivity_file = na.ScalarArray(transmissivity_file, axes=axis_layers)
 
     substrate = optika.chemicals.Chemical(material_substrate)
-    index_refraction_substrate = na.interp(
-        x=wavelength_ambient,
-        xp=substrate.index_refraction.inputs,
-        fp=substrate.index_refraction.outputs,
-    )
-    wavenumber_substrate = na.interp(
-        x=wavelength_ambient,
-        xp=substrate.wavenumber.inputs,
-        fp=substrate.wavenumber.outputs,
-    )
-    n_substrate = index_refraction_substrate + wavenumber_substrate * 1j
+    n_substrate = substrate.n(wavelength_ambient)
 
     reflectivity, transmissivity = optika.materials.multilayer_efficiency(
         material_layers=material_layers,

optika/sensors/_materials/_materials.py

@@ -340,17 +340,7 @@ def quantum_efficiency_effective(
 
     if n_substrate is None:
         substrate = optika.chemicals.Chemical("Si")
-        index_refraction_substrate = na.interp(
-            x=wavelength,
-            xp=substrate.index_refraction.inputs,
-            fp=substrate.index_refraction.outputs,
-        )
-        wavenumber_substrate = na.interp(
-            x=wavelength,
-            xp=substrate.wavenumber.inputs,
-            fp=substrate.wavenumber.outputs,
-        )
-        n_substrate = index_refraction_substrate + wavenumber_substrate * 1j
+        n_substrate = substrate.n(wavelength)
 
     if normal is None:
         normal = na.Cartesian3dVectorArray(0, 0, -1)
@@ -420,23 +410,16 @@ def index_refraction(
         self,
         rays: optika.rays.AbstractRayVectorArray,
     ) -> na.ScalarLike:
-        index_refraction = self._chemical.index_refraction
-        return na.interp(
-            x=rays.wavelength,
-            xp=index_refraction.inputs,
-            fp=index_refraction.outputs,
-        )
+        result = self._chemical.index_refraction(rays.wavelength)
+        return result
 
     def attenuation(
         self,
         rays: optika.rays.AbstractRayVectorArray,
     ) -> na.ScalarLike:
-        wavenumber = self._chemical.wavenumber
-        return na.interp(
-            x=rays.wavelength,
-            xp=wavenumber.inputs,
-            fp=4 * np.pi * wavenumber.outputs / wavenumber.inputs,
-        )
+        result = self._chemical.wavenumber(rays.wavelength)
+        result = 4 * np.pi * result / rays.wavelength
+        return result
 
     @property
     def is_mirror(self) -> bool: