add Inkstone to tests. reformat tests/examples

.github/workflows/testing.yml

@@ -16,10 +16,10 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: ['3.9', '3.10']
-        exclude:
-          - os: windows-latest
-            python-version: '3.9'
+        python-version: ['3.9', '3.10', '3.11', '3.12']
+#        exclude:
+#          - os: windows-latest
+#            python-version: '3.9'
 
     steps:
     - uses: actions/checkout@v1

examples/grating_pyramids_OPTOS.py

@@ -42,7 +42,7 @@
 options.nx = 25
 options.ny = 25
 options.parallel = True
-options.RCWA_method = "S4"
+options.RCWA_method = "Inkstone"
 
 # materials with constant n, zero k
 x = 1000

rayflare/rigorous_coupled_wave_analysis/rcwa.py

@@ -839,10 +839,23 @@ def initialise_S_inkstone(
                     )
                 # ? argument are different than S4
                 elif shape["type"] == "polygon":
+
+                    vertices = shape["vertices"]
+                    if type(vertices) is tuple:
+                        vertices = list(vertices)
+
+                    if 'center' in shape:
+                        center = shape['center']
+                        for i2, v in enumerate(vertices):
+                            vertices[i2] = [v[0] + center[0], v[1] + center[1]]
+
+                    if 'angle' in shape:
+                        logger.warn('Angle not implemented for polygon shapes in Inkstone')
+
                     S.AddPatternPolygon(
                         layer_name,
                         mat_name,
-                        shape["vertices"],
+                        vertices,
                     )
 
     return S
@@ -1335,36 +1348,94 @@ def get_fourier_epsilon(
 
         wl_ind = np.argmin(np.abs(self.current_wavelengths * 1e9 - wavelength))
 
-        if extent is None:
-            xdim = np.max(abs(np.array(self.size)[:, 0]))
-            ydim = np.max(abs(np.array(self.size)[:, 1]))
-            xs = np.linspace(-1.5 * xdim, 1.5 * xdim, n_points)
-            ys = np.linspace(-1.5 * ydim, 1.5 * ydim, n_points)
+        S = self.make_S(options, wl_ind)
 
-        else:
-            xs = np.linspace(extent[0][0], extent[0][1], n_points)
-            ys = np.linspace(extent[1][0], extent[1][1], n_points)
+        if options["RCWA_method"].lower() == "s4":
 
-        xys = np.meshgrid(xs, ys, indexing="ij")
+            if extent is None:
+                xdim = np.max(abs(np.array(self.size)[:, 0]))
+                ydim = np.max(abs(np.array(self.size)[:, 1]))
+                xs = np.linspace(-1.5 * xdim, 1.5 * xdim, n_points)
+                ys = np.linspace(-1.5 * ydim, 1.5 * ydim, n_points)
+
+            else:
+                xs = np.linspace(extent[0][0], extent[0][1], n_points)
+                ys = np.linspace(extent[1][0], extent[1][1], n_points)
 
-        a_r = np.zeros(len(xs) * len(ys))
-        a_i = np.zeros(len(xs) * len(ys))
+            xys = np.meshgrid(xs, ys, indexing="ij")
 
-        S = self.make_S(options, wl_ind)
+            a_r = np.zeros(len(xs) * len(ys))
+            a_i = np.zeros(len(xs) * len(ys))
 
-        if layer_index > 0:
-            depth = np.cumsum([0] + self.widths[1:-1] + [0])[layer_index - 1] + 1e-10
+            if layer_index > 0:
+                depth = np.cumsum([0] + self.widths[1:-1] + [0])[layer_index - 1] + 1e-10
+
+            else:
+                depth = -1
+
+            for i, (xi, yi) in enumerate(zip(xys[0].flatten(), xys[1].flatten())):
+                calc = S.GetEpsilon(xi, yi, depth)
+                a_r[i] = np.real(calc)
+                a_i[i] = np.imag(calc)
+
+            a_r = a_r.reshape((len(xs), len(ys)))
+            a_i = a_i.reshape((len(xs), len(ys)))
 
         else:
-            depth = -1
+            if extent is not None:
+                logger.warn('Plotting extent not implemented for Inkstone RCWA. Returns results for one unit cell.')
+
+            xs, ys, epsilon, mu = S.ReconstructLayer(f'layer_{layer_index + 1}', n_points, n_points)
+
+            # assume material is isotropic! top left entry of dielectric tensor
+            a_r = np.real(epsilon)[:, :, 0, 0]
+            a_i = np.imag(epsilon)[:, :, 0, 0]
+
+        if plot:
+            fig, axs = plt.subplots(1, 2, figsize=(7, 2.6))
+            im1 = axs[0].pcolormesh(xs, ys, a_r.T, cmap="magma")
+            fig.colorbar(im1, ax=axs[0])
+            axs[0].set_xlabel("x (nm)")
+            axs[0].set_ylabel("y (nm)")
+            axs[0].set_aspect(aspect=1)
+
+            im2 = axs[1].pcolormesh(xs, ys, a_i.T, cmap="magma")
+            fig.colorbar(im2, ax=axs[1])
+            axs[1].set_xlabel("x (nm)")
+            axs[1].set_ylabel("y (nm)")
+            axs[1].set_aspect(aspect=1)
+            plt.show()
+
+        return xs, ys, a_r, a_i
+
+    def get_fourier_epsilon_inkstone(
+        self, layer_index, wavelength, options, extent=None, n_points=200, plot=True
+    ):
+        """
+        Get the Fourier-decomposed epsilon scanning across x-y points for some layer in the structure for the number
+        of order specified in the options for the structure. Can also plot this automatically.
 
-        for i, (xi, yi) in enumerate(zip(xys[0].flatten(), xys[1].flatten())):
-            calc = S.GetEpsilon(xi, yi, depth)
-            a_r[i] = np.real(calc)
-            a_i[i] = np.imag(calc)
+        :param layer_index: index of the layer in which to get epsilon. layer 0 is the incidence medium, layer 1 is the first layer in the stack, etc.
+        :param wavelength: wavelength (in nm) at which to get epsilon
+        :param options: dictionary or State object containing user options
+        :param extent: range of x/y values in format [[x_min, x_max], [y_min, y_max]]. Default is 'None', will choose a reasonable area based \
+        on the unit cell size by default
+        :param n_points: number of points to scan across in the x and y directions
+        :param plot: plot the results (True or False, default True)
+
+        :return: xs, ys, a_r, a_i. The x points, y points, and the real and imaginary parts of the dielectric function.
+        """
+
+        wl_ind = np.argmin(np.abs(self.current_wavelengths * 1e9 - wavelength))
+
+        S = self.make_S(options, wl_ind)
+
+        xs, ys, epsilon, mu = S.ReconstructLayer(f'layer_{layer_index + 1}', n_points, n_points)
+
+        # assume material is isotropic! top left entry of dielectric tensor
+        a_r = np.real(epsilon)[:,:, 0, 0]
+        a_i = np.imag(epsilon)[:,:, 0, 0]
 
-        a_r = a_r.reshape((len(xs), len(ys)))
-        a_i = a_i.reshape((len(xs), len(ys)))
 
         if plot:
             fig, axs = plt.subplots(1, 2, figsize=(7, 2.6))
@@ -1383,6 +1454,8 @@ def get_fourier_epsilon(
 
         return xs, ys, a_r, a_i
 
+
+
     def get_fields(
         self,
         layer_index,

rayflare/utilities.py

@@ -287,7 +287,7 @@ def get_matrices_or_paths(
                 prof_dataset = xr.load_dataset(prof_mat_path)
                 return [existing, [full_mat, A_mat, prof_dataset]]
             else:
-                print("Recalculating with absorption profile information")
+                logger.info("Recalculating with absorption profile information")
                 existing = False
                 return [existing, [savepath_RT, savepath_A, prof_mat_path]]