add detailed_rcwa option for R/T power and amplitudes per order

examples/aSi_rods.py

@@ -9,6 +9,7 @@
 from time import time
 
 # Paper: https://doi.org/10.1002/adom.201700585
+# We will replicate Figure 1 from this paper
 
 # Define materials:
 MgF2 = material("MgF2")()  # MgF2 (SOPRA database)
@@ -25,15 +26,15 @@
 # the actual number of orders used can be slightly higher or lower than this. A higher number should
 # give a more accurate result, since the Fourier transformed representation is more accurate, but
 # takes longer to calculate.
-options.A_per_order = True # this makes the RCWA implementation (S4) calculate the power per diffraction order
+options.detailed_rcwa = True # this makes the RCWA implementation (S4) calculate the power per diffraction order
 # and field amplitudes, which are required to calculate the phase. By default, if this option is not turned on,
 # only the total R and T will be calculated.
 options.pol = 's' # polarization of the incident plane wave. This can be 's', 'p', 'u' (equal mixture of s and p),
 # or a Jones vector such as options.pol = (1/np.sqrt(2), 1j/np.sqrt(2)) (circular polarization)
 
-rad_list = np.linspace(100, 1000, 41) # list of pillar radii to scan through
+rad_list = np.linspace(100, 1000, 11) # list of pillar radii to scan through
 
-lattice_constant = np.linspace(1000, 3000, 61)
+lattice_constant = np.linspace(1000, 3000, 21)
 # list of lattice constants to scan through (distance between centre of pillars)
 
 # 21 radii, 31 lattice constants: miniconda python takes 101 s

rayflare/options.py

@@ -32,6 +32,7 @@ def __init__(self):
 
         # RCWA options
         self.A_per_order = False
+        self.detailed_rcwa = False
         self.S4_options = dict(
             LatticeTruncation="Circular",
             DiscretizedEpsilon=False,

rayflare/rigorous_coupled_wave_analysis/rcwa.py

@@ -1128,6 +1128,7 @@ def calculate(self, options):
                 self.size,
                 options["orders"],
                 options["A_per_order"],
+                options["detailed_rcwa"],
                 self.user_options,
             )
             for i1 in range(len(wl))
@@ -1139,37 +1140,35 @@ def calculate(self, options):
 
         self.rat_output_A = np.sum(A_mat, 1)  # used for profile calculation
 
+        S_for_orders = self.make_S(options, 0)
+
+        basis_set = S_for_orders.GetBasisSet()
+        f_mat = S_for_orders.GetReciprocalLattice()
+
+        results = {
+            "R": R,
+            "T": T,
+            "A_per_layer": A_mat,
+            "basis_set": basis_set,
+            "reciprocal": f_mat,
+        }
+
         if options["A_per_order"]:
             A_order = np.real(np.stack([item[3] for item in allres]))
+            results["A_layer_order"] = A_order,
+
+        if options["detailed_rcwa"]:
             R_order = np.real(np.stack([item[4] for item in allres]))
             T_order = np.real(np.stack([item[5] for item in allres]))
             R_amplitudes = np.stack([item[6] for item in allres])
             T_amplitudes = np.stack([item[7] for item in allres])
 
-            S_for_orders = self.make_S(options, 0)
-
-            basis_set = S_for_orders.GetBasisSet()
-            f_mat = S_for_orders.GetReciprocalLattice()
-
-            results = {
-                "R": R,
-                "T": T,
-                "A_per_layer": A_mat,
-                "A_layer_order": A_order,
-                "R_per_order": R_order,
-                "T_per_order": T_order,
-                "R_amplitudes": R_amplitudes,
-                "T_amplitudes": T_amplitudes,
-                "basis_set": basis_set,
-                "reciprocal": f_mat,
-            }
-            self.results = results
-            return results
+            results["R_per_order"] = R_order
+            results["T_per_order"] = T_order
+            results["R_amplitudes"] = R_amplitudes
+            results["T_amplitudes"] = T_amplitudes
 
-        else:
-            results = {"R": R, "T": T, "A_per_layer": A_mat}
-            self.results = results
-            return results
+        return results
 
     def calculate_profile(self, options):
         """It calculates the absorbed energy density within the material.
@@ -1589,19 +1588,26 @@ def RCWA_structure_wl(
     size,
     orders,
     A_per_order,
+    detailed_rcwa,
     S4_options,
 ):
     def vs_pol(s, p):
         S.SetExcitationPlanewave((theta, phi), s, p, 0)
         S.SetFrequency(1 / wl)
         R, T, incpf = rcwa_rt(S, len(widths))
         A_layer = rcwa_absorption_per_layer(S, len(widths), incpf)
+
+        result_list = [R, T, A_layer]
+
         if A_per_order:
             A_per_layer_order = rcwa_absorption_per_layer_order(S, len(widths), incpf)
+            result_list.append(A_per_layer_order)
+
+        if detailed_rcwa:
             R_per_order, T_per_order, _, R_amplitudes, T_amplitudes = rcwa_rt_pfbo(S, len(widths), incpf)
-            return R, T, A_layer, A_per_layer_order, R_per_order, T_per_order, R_amplitudes, T_amplitudes
-        else:
-            return R, T, A_layer
+            result_list += [None, R_per_order, T_per_order, R_amplitudes, T_amplitudes]
+
+        return result_list
 
     S = initialise_S(size, orders, geom_list, layers_oc, shapes_oc, s_names, widths, S4_options)
 
@@ -1611,7 +1617,8 @@ def vs_pol(s, p):
 
 
 def RCWA_structure_wl_inkstone(
-    wl, geom_list, layers_oc, shapes_oc, s_names, pol, theta, phi, widths, size, orders, *args
+    wl, geom_list, layers_oc, shapes_oc, s_names, pol,
+        theta, phi, widths, size, orders, *args
 ):
     def vs_pol(s, p):
         S.SetExcitation(theta, phi, s, p)