lib = gdstk.Library()

cell = lib.new_cell("BULLSEYE")

central_inner_radius = 0.355 + 0.027190402264497166
ring_width = 0.07 + np.array(
    [
        -0.04835278,
        -0.01293301,
        0.03173962,
        0.02891523,
        0.02110662,
        0.02439804,
        0.02123075,
    ]
)  # rings' width
ring_height = 0.2
theta1 = 3 * np.pi / 180
theta2 = 87 * np.pi / 180
ringsNumbers = 7
vertexIncrement = 0.1075
ringNumber = 7
xOrigin = 0
yOrigin = 0

inner_radius = central_inner_radius

for i in range(ringsNumbers):
    outer_radius = inner_radius + ring_width[i]
    for j in range(4):
        gcSlice = gdstk.ellipse(
            (xOrigin, yOrigin),
            inner_radius,
            outer_radius,
            theta1 + j * np.pi / 2,
            theta2 + j * np.pi / 2,
            layer=1,
            datatype=1,
            tolerance=0.0001,
        )
        cell.add(gcSlice)

    inner_radius = outer_radius + vertexIncrement

gc_etch = td.Geometry.from_gds(
    cell, gds_layer=1, axis=2, slab_bounds=(-ring_height / 2, ring_height / 2)
)
fig = plt.figure()
ax = gc_etch.plot(z=0, ax=fig.add_subplot(111))

ax.set_xlabel("x (µm)")
ax.set_ylabel("y (µm)")
plt.show()

sim = td.Simulation(
    center=(0, 0, 0),
    size=sim_size,
    grid_spec=td.GridSpec.auto(
        wavelength=(lambda_min + lambda_max) / 2, min_steps_per_wvl=automesh_per_wvl
    ),
    run_time=t_stop,
    sources=[dipole_source],
    monitors=time_monitors + [far_field_monitor],
    structures=[slab, etch],
    medium=vacuum,
    shutoff=1e-7,
    symmetry=[-1, 1, 0],
)

sim_data = web.run(simulation=sim, task_name="Aperiodic ring cavity")
print(sim_data.log)

resonance_finder = ResonanceFinder(freq_window=(monitor_freq[-1], monitor_freq[0]))
resonance_data = resonance_finder.run(signals=sim_data.data)
resonance_data.to_dataframe()

fig, ax = plt.subplots(1, 2, tight_layout=True, figsize=(12, 4))

time_monitor_names = [monitor.name for monitor in time_monitors]

time_response_Ex = np.zeros_like(sim_data["time_monitor_1"].Ex.squeeze())
time_response_Ey = np.zeros_like(sim_data["time_monitor_1"].Ey.squeeze())

for monitor_name in time_monitor_names:
    time_response_Ex += sim_data[monitor_name].Ex.squeeze()
    time_response_Ey += sim_data[monitor_name].Ey.squeeze()

time_response_Ex /= len(time_monitor_names)
time_response_Ey /= len(time_monitor_names)

# Perform FFT and get the frequency response
freq_response_Ex = np.abs(np.fft.fft(time_response_Ex))
freq_response_Ey = np.abs(np.fft.fft(time_response_Ey))

# Calculate the frequency values
freqs_Ex = np.linspace(0, 1 / sim_data.simulation.dt, len(time_response_Ex))
freqs_Ey = np.linspace(0, 1 / sim_data.simulation.dt, len(time_response_Ey))

# Calculate the corresponding wavelengths
with np.errstate(divide="ignore"):
    wavelengths_Ex = td.constants.C_0 / freqs_Ex
    wavelengths_Ey = td.constants.C_0 / freqs_Ey

# Select the indices within the monitor frequency range
plot_inds_Ex = np.where((monitor_freq[-1] < freqs_Ex) & (freqs_Ex < monitor_freq[0]))
plot_inds_Ey = np.where((monitor_freq[-1] < freqs_Ey) & (freqs_Ey < monitor_freq[0]))

# Plot the frequency response with respect to wavelength
ax[0].plot(wavelengths_Ex[plot_inds_Ex], freq_response_Ex[plot_inds_Ex])
ax[0].set_xlabel("Wavelength (um)")
ax[0].set_ylabel("Amplitude")
ax[0].set_title("Ex Frequency Response vs. Wavelength")
ax[0].grid(True)

ax[1].plot(wavelengths_Ey[plot_inds_Ey], freq_response_Ey[plot_inds_Ey])
ax[1].set_xlabel("Wavelength (um)")
ax[1].set_ylabel("Amplitude")
ax[1].set_title("Ey Frequency Response vs. Wavelength")
ax[1].grid(True)

plt.show()

make_field_plot(phi_vals, theta_vals, intensity)

# Parameters
NA = 0.68
n = 1.0
theta_max = np.arcsin(NA / n)  # max angle for the Gaussian beam
theta_0 = theta_max  # beam divergence angle

# Grids: theta_vals, phi_vals, intensity are assumed already defined
phi_grid, theta_grid = np.meshgrid(phi_vals, theta_vals)

# Differential elements (assume uniform)
dtheta = np.diff(theta_vals)[0]
dphi = np.diff(phi_vals)[0]
dOmega = np.sin(theta_grid) * dtheta * dphi  # solid angle element

# Define normalized Gaussian beam intensity profile
E_gauss = np.exp(-((theta_grid / theta_0) ** 2))  # (amplitude)
I_gauss = E_gauss**2  # intensity (|E|^2)

# --- Normalize E_sim (optional if already normalized) ---
I_sim = intensity  # if intensity = |E|^2
E_sim = np.sqrt(I_sim)

# --- Numerator: |<E_sim | E_gauss>|^2 ---
overlap = np.sum(E_sim * E_gauss * dOmega)
numerator = np.abs(overlap) ** 2

# --- Denominator: energy normalization ---
norm_sim = np.sum(I_sim * dOmega)
norm_gauss = np.sum(I_gauss * dOmega)

# --- Efficiency ---
eta = numerator / (norm_sim * norm_gauss)
percentage_overlap = eta * 100

print(
    f"Percentage of power coupled into Gaussian beam (NA={NA}): {percentage_overlap:.2f}%"
)