def get_bend_eme(
    bends=[(1, 1)],
    wvgIn=2,
    wvgOut=2,
    width=0.6,
    thickness=0.4,
    plane_size=(6, 6),
    num_modes=160,
    target_neff=1.5,
    core_index=2,
    precision="double",
    store_port_modes=False,
    constraint=None,
):
    """
    Creates and returns an EME simulation for a bent waveguide structure.

    Parameters:
        bends (list of tuple): A list of tuples where each tuple represents a bend in the waveguide.
                                Each tuple consists of two values: length and bend radius. Default is [(1, 1)].
        wvgIn (float): The length of the straight waveguide section at the input. Default is 2.
        wvgOut (float): The length of the straight waveguide section at the output. Default is 2.
        width (float): The width of the waveguide. Default is 0.6.
        thickness (float): The thickness of the waveguide core. Default is 0.4.
        plane_size (tuple): The size of the simulation plane (height, width). Default is (6, 6).
        num_modes (int): The number of modes to compute for the simulation. Default is 160.
        target_neff (float): The target effective index for the modes. Default is 1.5.
        core_index (float): The refractive index of the core material. Default is 2.
        precision (str): The numerical precision for the simulation, either "single" or "double". Default is "double".
        store_port_modes (bool): Whether to store port modes. Default is False.
    """

    # default mode spec
    mode_spec = td.EMEModeSpec(
        num_modes=num_modes,
        bend_radius=None,
        bend_axis=1,
        num_pml=(12, 12),
        target_neff=target_neff,
        precision=precision,
    )

    # size of the EME simulation. Start with 2 to compensate the port offsets
    size = 2

    # add straight sections
    if wvgIn > 0:
        bends = [(wvgIn, None)] + bends
    if wvgOut > 0:
        bends = bends + [(wvgOut, None)]

    structures = []
    mode_specs = []

    # iterate through each section and create each mode spec
    for i, bend in enumerate(bends):
        length, radius = bend

        mode_specs.append(mode_spec.updated_copy(bend_radius=radius))

        if i == 0:
            boundaries = [length]
        else:
            boundaries.append(boundaries[-1] + length)

        size += length

    eme_grid_spec = td.EMEExplicitGrid(boundaries=boundaries[:-1], mode_specs=mode_specs)

    # monitor for visualizing fields
    eme_field_mon = td.EMEFieldMonitor(name="field", size=(td.inf, 0, td.inf), num_modes=1)

    # creating structures

    sio2 = td.material_library["SiO2"]["Palik_LowLoss"]
    core = td.Medium(permittivity=core_index**2)

    # creating waveguide structure
    structures = [
        td.Structure(
            attrs={},
            geometry=td.PolySlab(
                axis=2,
                sidewall_angle=0.0,
                slab_bounds=(-thickness / 2, thickness / 2),
                vertices=np.array(
                    [
                        [-1e15, -width / 2],
                        [-1e15, width / 2],
                        [1e15, width / 2],
                        [1e15, -width / 2],
                    ]
                ),
            ),
            medium=core,
        )
    ]

    # creating the simulation object
    sim = td.EMESimulation(
        size=(size, plane_size[0], plane_size[1]),
        medium=sio2,
        center=(size / 2 - 1, 0, 0),
        structures=structures,
        axis=0,
        freqs=[td.C_0 / 1.55],
        eme_grid_spec=eme_grid_spec,
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=20),
        monitors=[eme_field_mon],
        port_offsets=(1, 1),
        symmetry=(0, 0, 1),
        store_port_modes=store_port_modes,
        constraint=constraint,
    )

    return sim

# function for creating points for an arc section
def arc(length, radius, center=(0, 0), theta0=0):
    theta = theta0 + np.linspace(0, length, 1001) / radius
    x = radius * np.cos(theta - np.pi / 2)
    y = radius * np.sin(theta - np.pi / 2)

    x = x + (center[0] - x[0])
    y = y + (center[1] - y[0])
    return x, y, theta


# function for creating points for a straight section at a given angle
def line(length, theta, center=(0, 0)):
    x = [center[0], center[0] + length * np.cos(theta)]
    y = [center[1], center[1] + length * np.sin(theta)]

    return x, y


# function for creating a list of points representing the bend waveguide from an EME simulation object
def eme2Curve(sim):
    mode_specs = sim.eme_grid_spec.mode_specs
    boundaries = sim.eme_grid_spec.boundaries
    boundaries = np.append(0, boundaries)
    boundaries = np.append(boundaries, sim.size[0] - sum(sim.port_offsets))
    lengths = np.diff(boundaries)
    radius = [i.bend_radius for i in mode_specs]

    x, y = [], []

    cx, cy = 0, 0
    theta = 0
    for i, r in enumerate(radius):
        if r:
            a, b, t = arc(lengths[i], r, (cx, cy), theta0=theta)
            theta = t[-1]
        else:
            a, b = line(lengths[i], theta, (cx, cy))
        x = np.append(x, a)
        y = np.append(y, b)

        cx = x[-1]
        cy = y[-1]

    return x, y

def get_fdtd(x, y, sim, width=0.6, thickness=0.4):
    # output angle
    theta = np.arctan2(y[-1] - y[-2], x[-1] - x[-2])

    # simulation size
    sx = max(x) - min(x)
    sy = max(y) - min(y)

    # simulation center
    centerX = min(x) + sx / 2
    centerY = min(y) + sy / 2

    # adding extra size so the structures will not start or end at the PML boundaries
    x = np.append([-5, 0], x)
    y = np.append([0, 0], y)

    cx = x[-1]
    cy = y[-1]

    a, b = line(5, theta, (cx, cy))
    x = np.append(x, a)
    y = np.append(y, b)

    # creating the structure with gdstk
    import gdstk

    cell = gdstk.Cell("bends")
    cell.add(gdstk.FlexPath(x + 1j * y, width, layer=1, datatype=0))

    geo = td.PolySlab.from_gds(
        cell,
        gds_layer=1,
        axis=2,
        slab_bounds=(-thickness / 2, thickness / 2),
    )[0]

    # simulation size
    size = (sx + 4, sy + 5, 3)

    # structure and source
    fdtd_structure = td.Structure(geometry=geo, medium=sim.structures[0].medium)

    source = td.ModeSource(
        center=(0, 0, 0),
        size=(0, 3, 3),
        direction="+",
        mode_spec=td.ModeSpec(num_modes=1),
        source_time=td.GaussianPulse(freq0=td.C_0 / 1.55, fwidth=0.1 * td.C_0 / 1.55),
    )

    while theta >= 2 * np.pi:
        theta -= 2 * np.pi

    thetaDeg = 180 * theta / np.pi

    # adjusting monitor size for the output angle
    if (thetaDeg > -45) and (thetaDeg <= 45):
        mon_size = (0, 3, 3)
    elif (thetaDeg > 45) and (thetaDeg <= 135):
        mon_size = (3, 0, 3)
    elif (thetaDeg > 135) and (thetaDeg <= 225):
        mon_size = (0, 3, 3)
    else:
        mon_size = (3, 0, 3)

    # monitors
    fieldMon = td.FieldMonitor(
        name="field", center=(0, 0, 0), size=(td.inf, td.inf, 0), freqs=[td.C_0 / 1.55]
    )

    # adjusting the theta angle of the output monitor
    if abs(theta - np.pi / 2) < 0.05 or abs(theta - np.pi) < 0.05:
        modeMon = td.ModeMonitor(
            name="modeMon",
            size=mon_size,
            center=(cx, cy, 0),
            freqs=[td.C_0 / 1.55],
            mode_spec=source.mode_spec,
        )

    else:
        # adjusting the correct angle for each quadrant
        if thetaDeg <= 90:
            t = theta
        if (thetaDeg > 90) and (thetaDeg <= 180):
            t = -(theta - np.pi / 2)
        if (thetaDeg > 180) and (thetaDeg <= 270):
            t = (theta - np.pi) + theta - np.pi / 2
        if (thetaDeg > 270) and (thetaDeg <= 360):
            t = -(theta - 3 * np.pi / 2)

        modeMon = td.ModeMonitor(
            center=(cx, cy, 0),
            size=mon_size,
            freqs=[td.C_0 / 1.55],
            mode_spec=td.ModeSpec(
                num_modes=5,
                angle_theta=t,
                angle_phi=0,
            ),
            name="modeMon",
        )

    sim_fdtd = td.Simulation(
        size=size,
        center=(centerX, centerY, sim.center[2]),
        medium=sim.medium,
        grid_spec=sim.grid_spec,
        sources=[source],
        monitors=[fieldMon, modeMon],
        run_time=10e-12,
        structures=[fdtd_structure],
        boundary_spec=td.BoundarySpec(
            x=td.Boundary.absorber(), y=td.Boundary.absorber(), z=td.Boundary.pml()
        ),
        symmetry=(0, 0, 1),
    )

    return sim_fdtd

def get_theta(length, offset):
    from scipy.optimize import fsolve

    def get_params(params, *args):
        theta = params
        length, offset = args
        r = (length / 2) / np.cos(theta - np.pi / 2)
        return abs(offset / 2 - (r * np.sin(theta - np.pi / 2) + r))

    theta = fsolve(get_params, (1,), args=(length, offset))[0]
    return (length / 2) / np.cos(theta - np.pi / 2), theta

lengths = [20, 40, 60, 80, 100]

eme_sims = []
fdtd_sims = []
for length in lengths:
    r, theta = get_theta(length, 20)
    eme_sim = get_bend_eme(
        wvgIn=2,
        wvgOut=2,
        bends=[(theta * r, r), (theta * r, -r)],
        plane_size=(17, 17),
        num_modes=160,
    )
    x, y = eme2Curve(eme_sim)
    fdtd_sim = get_fdtd(x, y, eme_sim)
    eme_sims.append(eme_sim)
    fdtd_sims.append(fdtd_sim)

eme_data, fdtd_data = web.run([eme_sims, fdtd_sims], verbose=True)

transmittance_eme = [
    float(d.smatrix.S21.isel(mode_index_in=0, mode_index_out=0, f=0).abs.item() ** 2) for d in eme_data
]
transmittance_fdtd = [
    float(d["modeMon"].amps.sel(direction="+", mode_index=0).abs.item() ** 2) for d in fdtd_data
]

fig, ax = plt.subplots()
ax.plot([20, 40, 60, 80, 100], transmittance_fdtd, "o", label="FDTD")
ax.plot([20, 40, 60, 80, 100], transmittance_eme, "*", label="EME")
ax.set_ylabel("transmittance")
ax.set_xlabel("bend length")

ax.legend()

plt.show()

def euler_curve(A=2.4, L=10, num_points=1000):
    from scipy import integrate

    Ls = np.linspace(0, L, num_points)  # L at (x1,y1)
    x1 = np.zeros(len(Ls))  # x coordinate of the clothoid curve
    y1 = np.zeros(len(Ls))  # y coordinate of the clothoid curve

    # compute x1 and y1 using the above integral equations

    y1 = np.array([integrate.quad(lambda theta: A * np.sin(theta**2 / 2), 0, L / A)[0] for L in Ls])
    x1 = np.array([integrate.quad(lambda theta: A * np.cos(theta**2 / 2), 0, L / A)[0] for L in Ls])

    return x1, y1


def euler_bend(offset, length, num_points=1000):
    from scipy.optimize import fsolve

    def get_params(params):
        A, L = params

        x1, y1 = euler_curve(A=A, num_points=num_points, L=L)

        return (abs(x1[-1] - length / 2) ** 2, abs(y1[-1] - offset / 2) ** 2)

    A, L = fsolve(get_params, (1, 1))
    x1, y1 = euler_curve(A=A, num_points=num_points, L=L)

    x = np.append(x1, np.flip(-x1 + 2 * x1[-1]))
    y = np.append(y1, np.flip(-y1 + 2 * y1[-1]))

    return x, y

def get_local_radius(x, y):
    yd = np.gradient(y, x, edge_order=2)
    ydd = np.gradient(yd, x, edge_order=2)
    R = ((1 + yd**2) ** (3 / 2)) / ydd

    return R

x_euler, y_euler = euler_bend(offset=5, length=30)
local_radius = get_local_radius(x_euler, y_euler)

fig, ax = plt.subplots()
ax.set_title("Euler S-bend")
ax.plot(x_euler, y_euler)

plt.show()

fig, ax = plt.subplots()
ax.set_title("Local curvature radius")
ax.plot(x_euler, local_radius)

plt.show()

dR = np.gradient(local_radius, x_euler, edge_order=2)
fig, ax = plt.subplots()
ax.set_title("Local radius derivative")
ax.plot(x_euler, abs(dR))

plt.show()

# function to create evenly distributed grid cells into one particular section
def get_positions(points, n_grid):
    if n_grid >= len(points):
        return np.arange(0, len(points), 1)
    indices = np.linspace(0, len(points) - 1, n_grid, dtype=int)
    return indices


def get_bends(x_euler, y_euler, grid_per_section=5, plot=False):
    # local radius and its derivatives
    local_radius = get_local_radius(x_euler, y_euler)
    dR = np.gradient(local_radius, x_euler, edge_order=2)

    # slice the data in chunks based on mean and standard deviation
    mean = np.nanmean(abs(dR))
    std = np.nanstd(abs(dR))

    section1 = abs(dR) >= mean + std
    section2 = (abs(dR) < mean + std) & (abs(dR) >= mean)
    section3 = abs(dR) < mean

    section1_points = get_positions(x_euler[section1], grid_per_section)
    section2_points = get_positions(x_euler[section2], grid_per_section)
    section3_points = get_positions(x_euler[section3], grid_per_section)

    # discontinuity points in the derivative that are continuous in the local radius curve
    discontinuity = x_euler[np.isnan(dR)][~np.isnan(local_radius[np.isnan(dR)])]
    discontinuity_y = y_euler[np.isnan(dR)][~np.isnan(local_radius[np.isnan(dR)])]
    discontinuity_radius = local_radius[np.isnan(dR)][~np.isnan(local_radius[np.isnan(dR)])]

    # retrieving the x position and local radius of the EME grid
    x_positions = np.concatenate(
        [
            x_euler[section1][section1_points],
            x_euler[section2][section2_points],
            x_euler[section3][section3_points],
            discontinuity,
            [x_euler[-1]],
        ]
    )
    sort = np.argsort(x_positions)
    x_positions = x_positions[sort]
    radius = np.concatenate(
        [
            local_radius[section1][section1_points],
            local_radius[section2][section2_points],
            local_radius[section3][section3_points],
            discontinuity_radius,
            [local_radius[-1]],
        ]
    )
    radius = radius[sort][1:]

    # transform positions in length
    length = np.diff(x_positions)

    # creating the list of tuples
    bends = list(zip(length, radius))

    # plot the grid points
    if plot:
        fig, ax = plt.subplots()
        ax.plot(x_euler, y_euler)
        ax.plot(
            x_euler[section1][section1_points],
            y_euler[section1][section1_points],
            "o",
            label="section 1",
        )
        ax.plot(
            x_euler[section2][section2_points],
            y_euler[section2][section2_points],
            "o",
            label="section 2",
        )
        ax.plot(
            x_euler[section3][section3_points],
            y_euler[section3][section3_points],
            "o",
            label="section 3",
        )
        ax.plot(discontinuity, discontinuity_y, "*", label="discontinuity")

        ax.legend()

        plt.show()

    return bends

x_euler, y_euler = euler_bend(offset=5, length=10)
bends = get_bends(x_euler, y_euler, grid_per_section=10, plot=True)

grids = 10
euler_lengths = [30, 50, 80, 100]

eme_sims = []
fdtd_sims = []
for l in euler_lengths:
    x_euler, y_euler = euler_bend(offset=5, length=l)
    bends = get_bends(x_euler, y_euler, grid_per_section=grids, plot=False)

    eme_sim = get_bend_eme(wvgIn=2, wvgOut=2, bends=bends, plane_size=(15, 15), num_modes=160)
    fdtd_sim = get_fdtd(x_euler, y_euler, eme_sim)
    eme_sims.append(eme_sim)
    fdtd_sims.append(fdtd_sim)

eme_data, fdtd_data = web.run([eme_sims, fdtd_sims], verbose=True)

transmittance_eme = [
    float(d.smatrix.S21.isel(mode_index_in=0, mode_index_out=0, f=0).abs.item() ** 2) for d in eme_data
]
transmittance_fdtd = [
    float(d["modeMon"].amps.sel(direction="+", mode_index=0).abs.item() ** 2) for d in fdtd_data
]

# plot results
fig, ax = plt.subplots()

ax.plot([30, 50, 80, 100], transmittance_fdtd, "o", label="FDTD")
ax.plot([30, 50, 80, 100], transmittance_eme, "*", label="EME")
ax.set_ylabel("transmittance")
ax.set_xlabel("bend length")
ax.legend()

plt.show()

x_euler, y_euler = euler_bend(offset=80, length=150)
bends = get_bends(x_euler, y_euler, grid_per_section=grids, plot=False)

eme_sim = get_bend_eme(wvgIn=2, wvgOut=2, bends=bends, plane_size=(15, 15), num_modes=160)
fdtd_sim = get_fdtd(x_euler, y_euler, eme_sim)

eme_data, fdtd_data = web.run([eme_sim, fdtd_sim], verbose=True)

fdtd = float(fdtd_data["modeMon"].amps.sel(direction="+", mode_index=0).abs.item() ** 2)
eme = float(eme_data.smatrix.S21.isel(mode_index_in=0, mode_index_out=0, f=0).abs.item() ** 2)

print(f"fdtd transmittance: {fdtd}\neme transmittance: {eme}")