Silicon Microring Resonator

import gdsfactory as gf
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Step 1: Create parametric ring resonator
def create_ring_filter(radius=7, gap=0.15, width=0.45):
    """Create a single ring resonator drop filter"""
    c = gf.Component()
    
    # Create ring and bus waveguides
    ring = c.add_ref(gf.components.ring(radius=radius, width=width))
    bus = c.add_ref(gf.components.straight(length=2*radius+10))
    
    # Position bus waveguide with specified gap
    bus.movey(-radius - width/2 - gap - width/2)
    
    # Add ports
    c.add_port("in", port=bus.ports["o1"])
    c.add_port("through", port=bus.ports["o2"])
    c.add_port("drop", port=ring.ports["o2"])
    c.add_port("add", port=ring.ports["o1"])
    
    return c

# Step 2: FDTD simulation setup (using Tidy3D)
def simulate_ring(component, wavelengths):
    """Run FDTD simulation and extract S-parameters"""
    import tidy3d as td
    
    # Convert GDS to simulation geometry
    sim = gf.plugins.tidy3d.get_simulation(
        component=component,
        wavelength=1.55,
        wavelength_span=0.1,
        port_source="in",
        port_monitors=["through", "drop"],
        mesh_accuracy=3
    )
    
    # Run simulation
    data = td.web.run(sim, task_name="ring_resonator")
    
    # Extract S-parameters
    S21 = data["through"].amps
    S31 = data["drop"].amps
    
    return S21, S31

# Step 3: Analyze resonances
def lorentzian(x, x0, gamma, A, offset):
    """Lorentzian function for resonance fitting"""
    return A * gamma**2 / ((x - x0)**2 + gamma**2) + offset

def extract_metrics(wavelengths, transmission):
    """Extract Q-factor, FSR, and extinction ratio"""
    # Find resonances
    from scipy.signal import find_peaks
    peaks, _ = find_peaks(-transmission, height=0.1)
    
    # Fit Lorentzian to each resonance
    metrics = []
    for peak in peaks:
        # Select data around peak
        idx_range = slice(max(0, peak-20), min(len(wavelengths), peak+20))
        x_fit = wavelengths[idx_range]
        y_fit = transmission[idx_range]
        
        # Initial guess
        x0_guess = wavelengths[peak]
        gamma_guess = 0.0004  # ~0.8nm FWHM
        A_guess = 1 - np.min(y_fit)
        offset_guess = np.max(y_fit)
        
        # Fit
        popt, _ = curve_fit(lorentzian, x_fit, y_fit, 
                           p0=[x0_guess, gamma_guess, A_guess, offset_guess])
        
        # Calculate metrics
        resonance_wavelength = popt[0]
        fwhm = 2 * popt[1]
        Q_factor = resonance_wavelength / fwhm
        extinction_ratio_dB = -10 * np.log10(np.min(y_fit) / np.max(y_fit))
        
        metrics.append({
            'wavelength': resonance_wavelength,
            'Q': Q_factor,
            'ER_dB': extinction_ratio_dB,
            'FWHM_nm': fwhm * 1000  # Convert to nm
        })
    
    # Calculate FSR
    if len(metrics) > 1:
        fsr_nm = np.diff([m['wavelength'] for m in metrics]).mean() * 1000
        fsr_ghz = 3e8 / (1.55e-6)**2 * fsr_nm * 1e-9 * 1e9  # Convert to GHz
    else:
        fsr_nm = fsr_ghz = None
    
    return metrics, fsr_ghz

# Step 4: Parameter sweep
def sweep_design_space():
    """Sweep radius and gap to optimize performance"""
    radii = np.linspace(6, 8, 5)
    gaps = np.linspace(0.12, 0.25, 5)
    
    results = []
    for radius in radii:
        for gap in gaps:
            # Create design
            ring = create_ring_filter(radius=radius, gap=gap)
            
            # Simulate (simplified - would use actual FDTD)
            # Here we use analytical approximation
            Q_loaded = estimate_Q(radius, gap)
            ER = estimate_ER(gap)
            
            results.append({
                'radius': radius,
                'gap': gap,
                'Q': Q_loaded,
                'ER_dB': ER
            })
    
    return results

# Helper functions for analytical estimates
def estimate_Q(radius, gap):
    """Estimate loaded Q-factor"""
    # Coupling coefficient (empirical model)
    kappa = np.exp(-2 * gap / 0.1)  # Simplified exponential model
    
    # Intrinsic Q (loss-limited)
    Q_intrinsic = 50000  # Typical for low-loss SOI
    
    # Loaded Q
    Q_loaded = Q_intrinsic / (1 + kappa * Q_intrinsic)
    
    return Q_loaded

def estimate_ER(gap):
    """Estimate extinction ratio"""
    kappa = np.exp(-2 * gap / 0.1)
    ER_linear = (1 - kappa)**2 / (4 * kappa)
    return 10 * np.log10(ER_linear)

# Step 5: Generate compact model
def create_compact_model(ring_params, s_params):
    """Generate Caphe/INTERCONNECT compatible model"""
    model = {
        'type': 'ring_resonator',
        'parameters': ring_params,
        'ports': ['in', 'through', 'drop', 'add'],
        's_parameters': s_params,
        'center_wavelength': 1.55e-6,
        'ng': 4.2,  # Group index
        'loss_dB_per_cm': 2.0
    }
    
    return model

# Main execution
if __name__ == "__main__":
    # Create optimized design
    ring = create_ring_filter(radius=7.2, gap=0.18)
    
    # Export GDS
    ring.write_gds("ring_filter.gds")
    
    # Run parameter sweep
    results = sweep_design_space()
    
    # Find optimal design
    optimal = max(results, key=lambda x: x['ER_dB'] if x['Q'] > 1500 else 0)
    
    print(f"Optimal Design:")
    print(f"  Radius: {optimal['radius']:.1f} µm")
    print(f"  Gap: {optimal['gap']:.3f} µm")
    print(f"  Q-factor: {optimal['Q']:.0f}")
    print(f"  Extinction Ratio: {optimal['ER_dB']:.1f} dB")