Can Meep be used to simulate FSS (Frequency Selective Surface)?

[Brianye99]

I am currently researching FSS. Can I simulate plane waves in meep? In the final scene of my FSS, there are many periodic arrangements. How to set boundary conditions to simulate them? I need the S parameter of FSS, how can I implement it?
I have tried many times on my own, but the results have not been able to converge,here is an example code:
Guys, can you help me check this out? Please : )

import meep as mp
import numpy as np
import matplotlib.pyplot as plt
import gc  # Garbage collection

# ==============================================================================
# 1. Basic Parameters and Pattern
# ==============================================================================
unit_length = 1e-3  # 1 unit = 1 mm
resolution = 5      # Resolution
period_x = 30.0
period_y = 30.0
sub_h = 0.762       # Substrate thickness

# Generate FSS pattern (cross-shaped example)
fss_pattern = np.zeros((30, 30))
center_idx = 15
half_width = 3
fss_pattern[center_idx-half_width:center_idx+half_width, :] = 1
fss_pattern[:, center_idx-half_width:center_idx+half_width] = 1

# Frequency band settings (1-6 GHz)
GHz_to_Meep = 1e9 * unit_length / 299792458
f_min = 1.0 * GHz_to_Meep
f_max = 6.0 * GHz_to_Meep
f_cen = 0.5 * (f_min + f_max)
df = f_max - f_min
nfreq = 100

# Material definitions
fr4_material = mp.Medium(epsilon=4.4) # FR4 dielectric material
fss_material_grid = mp.MaterialGrid(mp.Vector3(period_x, period_y, 0),
                                    weights=fss_pattern,
                                    medium1=mp.air,
                                    medium2=mp.metal,
                                    do_averaging=False)

# ==============================================================================
# 2. Core Simulation Function (Runs for air and structure cases)
# ==============================================================================
def run_simulation(is_empty_run=False, flux_data_to_load=None):
    """
    is_empty_run: True = run air case (normalization), False = run structure case
    flux_data_to_load: If running structure case, pass the incident wave data from air run
    """
    
    # --- A. Dimensions and Boundaries ---
    dpml = 20
    pad_air = 40
    sz = 2 * dpml + 2 * pad_air + sub_h
    cell_size = mp.Vector3(period_x, period_y, sz)
    boundary_layers = [mp.PML(dpml, direction=mp.Z)]
    
    # --- B. Geometry Construction (Key Logic) ---
    geometry = []
    
    if is_empty_run:
        # Case 1: Air run
        # Add nothing (pure air) to calculate the pure incident wave intensity
        pass 
        
    else:
        # Case 2: Structure run
        # 1. Add substrate (slightly raised at the center to align the bottom with Z=0 for easier reasoning)
        #    Substrate range: Z in [0, 0.762]
        geometry.append(
            mp.Block(mp.Vector3(mp.inf, mp.inf, sub_h),
                     center=mp.Vector3(0, 0, sub_h/2), 
                     material=fr4_material)
        )
        
        # 2. Add FSS metal (use effective thickness!)
        eff_thickness = 1.0 / resolution  # 0.2mm, ensuring the grid can capture it
        
        # Metal is placed on top of the substrate surface (Z=0.762)
        # Metal center = 0.762 + 0.2/2 = 0.862
        metal_z = sub_h + eff_thickness/2
        
        geometry.append(
            mp.Block(mp.Vector3(period_x, period_y, eff_thickness),
                     center=mp.Vector3(0, 0, metal_z),
                     material=fss_material_grid)
        )

    # --- C. Dynamically Calculate Positions (Prevent Overlapping) ---
    # The value 0 here refers to the absolute coordinate origin. Since we place the substrate bottom at 0,
    # it's safer to readjust the relative positions based on sz (total height)
    # Simpler approach:
    
    # The structure region is approximately between [-2, 2] (including substrate and metal)
    # We place the light source at the bottom and the monitors in the middle
    
    source_pos = -sz/2 + dpml + 5.0      # 5mm above the bottom PML
    refl_pos   = source_pos + 10.0       # 10mm above the light source
    trans_pos  = sz/2 - dpml - 10.0      # 10mm below the top PML
    
    # --- D. Light Source and Simulation Object ---
    sources = [
        mp.Source(
            mp.GaussianSource(f_cen, fwidth=df),
            component=mp.Ex,
            center=mp.Vector3(0, 0, source_pos),
            size=mp.Vector3(period_x, period_y, 0)
        )
    ]

    sim = mp.Simulation(
        cell_size=cell_size,
        boundary_layers=boundary_layers,
        geometry=geometry,
        sources=sources,
        resolution=resolution,
        k_point=mp.Vector3(0,0,0),
        courant=0.3 # Stability factor
    )

    # --- E. Monitors ---
    # Define regions
    flux_region_r = mp.FluxRegion(center=mp.Vector3(0, 0, refl_pos), size=mp.Vector3(period_x, period_y, 0))
    flux_region_t = mp.FluxRegion(center=mp.Vector3(0, 0, trans_pos), size=mp.Vector3(period_x, period_y, 0))

    # Add monitors
    refl = sim.add_flux(f_cen, df, nfreq, flux_region_r)
    trans = sim.add_flux(f_cen, df, nfreq, flux_region_t)

    # If it's the second run, load the incident wave data from the first run (to subtract for pure reflection)
    if not is_empty_run and flux_data_to_load is not None:
        sim.load_minus_flux_data(refl, flux_data_to_load)

    # --- F. Run Simulation ---
    # Can use 1e-3 to accelerate here, as the effective thickness makes the results obvious
    sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(0, 0, trans_pos), 1e-3))

    # --- G. Data Collection and Cleanup ---
    if is_empty_run:
        # For air run, save the "incident wave data" for the next run
        inc_data = sim.get_flux_data(refl)
        inc_flux = mp.get_fluxes(trans)
        
        sim.reset_meep() # [Key] Reset the simulation environment
        return inc_data, inc_flux
    else:
        # For structure run, get the results directly
        refl_flux = mp.get_fluxes(refl)
        trans_flux = mp.get_fluxes(trans)
        
        sim.reset_meep() # [Key] Reset the simulation environment
        return refl_flux, trans_flux


# ==============================================================================
# 3. Main Execution Logic
# ==============================================================================
if __name__ == "__main__":
    
    # --- Phase 1: Run air case (normalization) ---
    if mp.am_master():
        print(">>> Phase 1: Normalization (Air run)...")
        
    inc_data, inc_flux = run_simulation(is_empty_run=True)

    # --- Phase 2: Run structure case ---
    if mp.am_master():
        print(">>> Phase 2: Device Simulation...")
        
    refl_flux, trans_flux = run_simulation(is_empty_run=False, flux_data_to_load=inc_data)

    # --- Phase 3: Calculation and Plotting (Master node only) ---
    if mp.am_master():
        print(">>> Calculating S-parameters...")
        
        freqs = np.linspace(f_min, f_max, nfreq) / GHz_to_Meep # Convert back to GHz
        
        # S11 = reflected / incident
        # Note: load_minus_flux is applied in run_simulation, so refl_flux is already (total field - incident field) = reflected field
        # We take the negative sign here due to Meep's flux direction definition
        R = -np.array(refl_flux) / np.array(inc_flux)
        T = np.array(trans_flux) / np.array(inc_flux)
        
        S11_dB = 10 * np.log10(np.abs(R)**2 + 1e-10) # Power spectrum uses square, or directly 20*log10(|E|)
        S21_dB = 10 * np.log10(np.abs(T)**2 + 1e-10) # Flux itself is power (Poynting), so use 10log directly

        # Plotting
        plt.figure(figsize=(10, 6))
        plt.plot(freqs, S11_dB, 'r-', linewidth=2, label='S11 (Reflection)')
        plt.plot(freqs, S21_dB, 'b-', linewidth=2, label='S21 (Transmission)')
        
        plt.xlabel('Frequency (GHz)')
        plt.ylabel('Magnitude (dB)')
        plt.title('FSS S-Parameters Simulation')
        plt.legend()
        plt.grid(True, which='both', alpha=0.3)
        plt.ylim(-50, 5)
        
        save_name = 'fss_final_result.png'
        plt.savefig(save_name)
        print(f">>> Done! Result saved to {save_name}")