fix tests, examples

examples/aSi_pillar_grating.py

@@ -9,7 +9,6 @@
 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)
@@ -26,15 +25,28 @@
 # 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.detailed_rcwa = True # this makes the RCWA implementation (S4) calculate the power per diffraction order
+options.A_per_order = 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, 11) # list of pillar radii to scan through
+n_radii = 101
+n_lat = 81
 
-lattice_constant = np.linspace(1000, 3000, 21)
+rad_list = np.linspace(100, 1000, n_radii) # list of pillar radii to scan through
+
+lattice_constant = np.linspace(1000, 3000, n_lat)
+
+rad = 1000/2
+x = 1000
+grating_circles = [{"type": "circle", "mat": aSi, "center": (0, 0), "radius": rad}]
+layers = [Layer(material=Air, width=si("2000nm"), geometry=grating_circles)]
+size = ((x, 0), (x / 2, np.sin(np.pi / 3) * x))
+S4_setup = rcwa_structure(layers, size=size, options=options, incidence=MgF2, transmission=Air)
+RAT = S4_setup.calculate(options)
+
+S4_setup.get_fourier_epsilon(1, 4000, options)
 # list of lattice constants to scan through (distance between centre of pillars)
 
 # 21 radii, 31 lattice constants: miniconda python takes 101 s
@@ -45,6 +57,8 @@
 # specific lattice constant. Will use this below to execute in parallel.
 def fixed_lattice_constant(x, rad_list):
 
+    print(x)
+
     # Unit cell vectors (real space) for the hexagonal layout. The unit cell is a triangle.
     size = ((x, 0), (x / 2, np.sin(np.pi / 3) * x))
 
@@ -60,17 +74,24 @@ def fixed_lattice_constant(x, rad_list):
         grating_circles = [{"type": "circle", "mat": aSi, "center": (0, 0), "radius": rad}]
         layers = [Layer(material=Air, width=si("2000nm"), geometry=grating_circles)]
 
+        if rad > x/2:
+            T_result[j1] = 0
+            T_amp[j1] = 0
+            inc_amp[j1] = 0
+
+
         # Make the object to do RCWA calculations. Assume light is incident from the MgF2 which is below
         # the pillars in the paper.
-        S4_setup = rcwa_structure(layers, size=size, options=options, incidence=MgF2, transmission=Air)
-        RAT = S4_setup.calculate(options)
+        else:
+            S4_setup = rcwa_structure(layers, size=size, options=options, incidence=MgF2, transmission=Air)
+            RAT = S4_setup.calculate(options)
 
-        # Store the results; we care about total transmission, and want the amplitudes of the incident
-        # light and the transmitted light (first entry in the last is always (0, 0) to calculate the
-        # phase change.
-        T_result[j1] = RAT["T"][0]
-        T_amp[j1] = RAT["T_amplitudes"][0][0][0]
-        inc_amp[j1] = RAT['R_amplitudes'][0][0][0]
+            # Store the results; we care about total transmission, and want the amplitudes of the incident
+            # light and the transmitted light (first entry in the last is always (0, 0) to calculate the
+            # phase change.
+            T_result[j1] = RAT["T"][0]
+            T_amp[j1] = RAT["T_amplitudes"][0][0][0]
+            inc_amp[j1] = RAT['R_amplitudes'][0][0][0]
 
     return T_result, T_amp, inc_amp
 
@@ -85,19 +106,34 @@ def fixed_lattice_constant(x, rad_list):
 T_amp = np.stack([item[1] for item in allres])
 inc_amp = np.stack([item[2] for item in allres])
 
+# save results for a specific number of radii/lattice constants:
+name_base = f'{n_radii}_rads_{n_lat}_lc'
+
+np.save(f'{name_base}_T_result.npy', T_result)
+np.save(f'{name_base}_T_amp.npy', T_amp)
+np.save(f'{name_base}_inc_amp.npy', inc_amp)
+
 # calculate phase and put it in the range 0, 2pi:
 phase_res = np.angle(T_amp) # this is in range -pi to pi
 phase_res[phase_res < 0] = 2*np.pi - np.abs(phase_res[phase_res < 0])
 
+fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 10))
+c = ax1.imshow(T_result, aspect=0.3,
+           extent=[rad_list[0], rad_list[-1], lattice_constant[0], lattice_constant[-1]],
+           origin="lower", cmap="Spectral_r")
+plt.colorbar(c, shrink=0.8)
+ax1.set_xlabel("Radius (nm)")
+ax1.set_ylabel("Lattice constant (nm)")
+ax1.set_title("Transmission")
 
-plt.figure()
-plt.imshow(phase_res/np.pi, aspect="auto",
+c = ax2.imshow(phase_res/np.pi, aspect=0.3,
            extent=[rad_list[0], rad_list[-1], lattice_constant[0], lattice_constant[-1]],
            origin="lower", cmap="Spectral_r")
-plt.colorbar()
-plt.xlabel("Radius (nm)")
-plt.ylabel("Lattice constant (nm)")
-plt.title(r"Phase (units of $\pi$)")
+plt.colorbar(c, shrink=0.8)
+ax2.set_xlabel("Radius (nm)")
+ax2.set_ylabel("Lattice constant (nm)")
+ax2.set_title(r"Phase (units of $\pi$)")
+plt.tight_layout()
 plt.show()
 
 max_radius = 610