Metasurface: microwave frequency selective surface (FSS)

import os
import numpy as np
import matplotlib.pyplot as plt

import gdsfactory as gf
from gdsfactory.technology import LogicalLayer, LayerLevel, LayerStack
import luminescent as lumi

# simulation folder
path = os.path.join("runs", "FSS")

# dims [mm]
# length units are arbitrary so long as they're consistent
P = 15  # periodicity
L = 9.4  # length of cross
W = 2.0  # width of cross stub
h_sub = 2.2  # substrate thickness
h_metal = 0.1  # metal thickness

# frequencies [GHz]
# spectrum computed from single FDTD run. more frequencies slightly increases runtime
frequency = 12
bandwidth = 4
nfs = 81
frequencies = np.linspace(frequency - bandwidth / 2, frequency + bandwidth / 2, nfs)

c0 = 300  # speed of light [mm GHz]

# incidence angle [degrees] from z in xz plane. E along y (TE)
theta = 0

def setup(path, frequency, frequencies, P, L, W, h_sub, h_metal, theta):
    # simulation parameters
    nres = 12
    dx = W / 16  # desired resolution near metal
    shutoff = 1e-3
    tmax = 1000  # max time [periods]
    PML_decay = 0.0001  # x y z
    courant = 0.9
    subpixel_smoothing = 1
    shot_every = 2.5  # movie [periods]
    views = [lumi.view("Ex",x=0, z=h_metal / 2)]

    # normal incidence: quarter unit cell [0, P/2]² with symmetry walls (E along y: PEC on x, PMC on y)
    # oblique incidence: one full unit cell, periodic in x y. Bloch wavevector sets incidence angle (exact at center frequency)
    sym = theta == 0
    if sym:
        boundaries = ["PEC", "PMC", "PML"]
        k_bloch = None
    else:
        boundaries = ["periodic", "periodic", "PML"]
        k_bloch = lumi.bloch_wavevector(theta, 0, dim=3, f=frequency, c0=c0)

    # margins
    # air gap keeps evanescent resonant field from coupling to PML
    gap = 30.0
    source_offset = 2.0

    # vertical dims
    z_o1 = h_metal + gap
    z_o2 = -h_sub - gap
    zmin = z_o2
    zmax = z_o1 + source_offset

    # materials
    # refine mesh in substrate and metal: n sets mesh density relative to vacuum
    materials = {
        "air": lumi.material(1),
        "sub": lumi.material(epsilon=2.5),
        "PEC": lumi.PEC_material(dx=dx),
    }
    background = "air"

    # layers
    METAL = 1, 0
    SUB = 2, 0
    layer_stack = LayerStack(
        layers={
            "metal": LayerLevel(
                layer=LogicalLayer(layer=METAL),
                thickness=h_metal,
                zmin=0,
                material="PEC",
                mesh_order=1,
            ),
            "sub": LayerLevel(
                layer=LogicalLayer(layer=SUB),
                thickness=h_sub,
                zmin=-h_sub,
                material="sub",
                mesh_order=2,
            ),
        }
    )

    # make component in gdsfactory. alternatively can import .gds into gdsfactory
    # cross unit cell centered at origin, clipped to quarter cell x, y >= 0 under symmetry
    lo = 0 if sym else -1

    def rect(a, b):
        return [(lo * a, lo * b), (a, lo * b), (a, b), (lo * a, b)]

    c = gf.Component()
    c.add_polygon(rect(W / 2, L / 2), layer=METAL)
    c.add_polygon(rect(L / 2, W / 2), layer=METAL)
    c.add_polygon(rect(P / 2, P / 2), layer=SUB)
    c.plot()

    bbox = c.bbox_np().tolist()
    bbox[0].append(zmin)
    bbox[1].append(zmax)

    # sources, monitors
    # plane wave on face frame spanning the cell. E along y, H along x
    # wave admittance in air: Hx/Ey = cos(theta) in magnitude (sign fixed automatically)
    frames = [lumi.face("plane_wave", modes=[[{"Ex": 1, "Hy": np.cos(np.deg2rad(theta))}]])]

    # dim: outward normal axis. source injects into cell, monitors record in (-) and out (+) waves
    sources = [
        lumi.source(
            "o1", z_o1 + source_offset, dim=3,
            frame="plane_wave",
            f=frequency,
            bandwidth=max(max(frequencies) - min(frequencies), frequency / 8),
        )
    ]

    monitors = [
        lumi.monitor("o1", z_o1, dim=3, frame="plane_wave"),
        lumi.monitor("o2", z_o2, dim=-3, frame="plane_wave"),
    ]

    # writes .obj mesh files from gdsfactory component
    lumi.mesh(path, c, layer_stack, materials, zmin, zmax)

    # writes setup file
    lumi.setup(
        path=path,
        f=frequency,
        fs=frequencies,
        c0=c0,
        materials=materials,
        background=background,
        boundaries=boundaries,
        # k_bloch=k_bloch,
        bbox=bbox,
        #
        frames=frames,
        sources=sources,
        monitors=monitors,
        #
        nres=nres,
        courant=courant,
        PML_decay=PML_decay,
        subpixel_smoothing=subpixel_smoothing,
        # 
        views=views,
        shot_every=shot_every,
        tmax=tmax,
        shutoff=shutoff,
    )


setup(path, frequency, frequencies, P, L, W, h_sub, h_metal, theta)

lumi.solve(path)

lumi.movie(path)

sol = lumi.load(path)
o1 = sol["monitors"]["o1"]
o2 = sol["monitors"]["o2"]

# + denotes outgoing mode (wrt port monitor outward normal), - incoming
S21 = o2["modes"][0]["+"] / o1["modes"][0]["-"]
S11 = o1["modes"][0]["+"] / o1["modes"][0]["-"]

plt.figure()
plt.plot(frequencies, 10 * np.log10(np.abs(S21)**2), label="S21")
plt.plot(frequencies, 10 * np.log10(np.abs(S11)**2), label="S11")
plt.xlabel("f [GHz]")
plt.ylabel("power [dB]")
plt.legend()
plt.show()