# All dimensions in µm

# --- GeSi waveguide (from paper, Fig. 1a) ---
gesi_w = 0.60  # GeSi waveguide width (600 nm)
gesi_h = 0.35  # GeSi waveguide height (350 nm)
gesi_intrinsic_w = 0.45  # intrinsic region width (450 nm, Section 4)
device_length_um = 40.0  # device length along propagation (Section 4)

# --- Si rib waveguide (SOI, 220 nm platform) ---
si_rib_h = 0.22  # full Si rib height (220 nm)
si_rib_w = 0.80  # Si rib width (800 nm)
si_slab_h = 0.14  # Si slab height outside the rib (140 nm) # estimated
# GeSi is partially embedded in the rib; lateral contacts are in the 140 nm slab region
si_half_w = 2.5  # half-width of Si slab (estimated: ~5 µm total)  # estimated

# --- Doping concentrations ---
# GeSi doping (from paper, Section 4)
gesi_p_doping = 7e16  # acceptor concentration in GeSi p-side [cm^-3]
gesi_n_doping = 7e16  # donor concentration in GeSi n-side [cm^-3]
# Acceptors present in the as-grown GeSi without any intentional doping step.
# The paper does not report their concentration; this value is an assumption.
gesi_bg_doping = 1e16  # background acceptor concentration in the GeSi [cm^-3]

# Si slab doping (imec 220 nm SOI platform)
si_p_doping = 5e17  # Si p+ slab region [cm^-3]
si_pp_doping = 1e19  # Si p++ contact region [cm^-3]
si_n_doping = 5e17  # Si n+ slab region [cm^-3]
si_nn_doping = 1e19  # Si n++ contact region [cm^-3]
si_pp_width = 0.5  # p++ region width from left edge [µm]
si_nn_width = 0.5  # n++ region width from right edge [µm]
si_junction_x = 0.0  # x where the Si p implant meets the Si n implant [µm]

# --- Defective layer at the GeSi/Si growth interface ---
# The paper reports pre-existing defects at the Ge/Si and Ge/Ox interfaces but
# gives no thickness or density, so both values below are assumptions. Only the
# bottom (growth) Ge/Si interface is modeled here, not the GeSi sidewalls.
damaged_h = 0.02  # assumed thickness of the defective layer [µm]
tau_damaged = 1e-11  # assumed SRH lifetime in that layer [s]

# --- BOX and cladding ---
box_h = 2.0  # SiO2 BOX thickness (typical 220 nm SOI)  # estimated
clad_h = 0.92  # SiO2 cladding above GeSi

# --- Metal contacts ---
contact_w = 0.3  # metal contact width  # estimated

# --- Bias voltages to simulate ---
# the biases of the measured trace of Fig. 8(a), plus 0 V for the unbiased state
voltages = np.arange(-3.0, 1.0, 0.25)

# --- Derived quantities ---
sim_half_w = si_half_w + contact_w  # simulation half-width
z_size = 0.5  # structure z-extent for 2D simulation

# y-coordinates of key interfaces
y_bot = 0.0  # BOX / Si interface (bottom of Si)
y_slab_top = si_slab_h  # top of Si slab (= 0.14 µm)
y_rib_top = si_rib_h  # top of Si rib (= 0.22 µm)
y_gesi_top = si_rib_h / 2 + gesi_h  # top of GeSi (= 0.46 µm)

# --- GeSi doping (from paper) ---
# The boxes span the full GeSi height, from the embedded bottom (y_rib_top / 2) to the
# top (y_gesi_top), so the lower GeSi section inside the Si rib is doped as well.
p_doping_gesi = td.ConstantDoping.from_bounds(
    concentration=gesi_p_doping,
    rmin=(-gesi_w / 2, y_rib_top / 2 - 1e-4, -z_size),
    rmax=(-gesi_intrinsic_w / 2, y_gesi_top + 1e-4, z_size),
)
n_doping_gesi = td.ConstantDoping.from_bounds(
    concentration=gesi_n_doping,
    rmin=(+gesi_intrinsic_w / 2, y_rib_top / 2 - 1e-4, -z_size),
    rmax=(+gesi_w / 2, y_gesi_top + 1e-4, z_size),
)
# unintentional acceptor background over the whole GeSi film; the wing implants add on top
bg_doping_gesi = td.ConstantDoping.from_bounds(
    concentration=gesi_bg_doping,
    rmin=(-gesi_w / 2 - 1e-4, y_rib_top / 2 - 1e-4, -z_size),
    rmax=(+gesi_w / 2 + 1e-4, y_gesi_top + 1e-4, z_size),
)

# --- Si doping ---

p_doping_si = td.ConstantDoping.from_bounds(
    concentration=si_p_doping,
    rmin=(-si_half_w, y_bot - 1e-4, -z_size),
    rmax=(si_junction_x, y_rib_top + 1e-4, z_size),
)
pp_doping_si = td.ConstantDoping.from_bounds(
    concentration=si_pp_doping,
    rmin=(-si_half_w, y_bot - 1e-4, -z_size),
    rmax=(-si_half_w + si_pp_width, y_rib_top + 1e-4, z_size),  # full Si height
)
n_doping_si = td.ConstantDoping.from_bounds(
    concentration=si_n_doping,
    rmin=(si_junction_x, y_bot - 1e-4, -z_size),
    rmax=(+si_half_w, y_rib_top + 1e-4, z_size),
)
nn_doping_si = td.ConstantDoping.from_bounds(
    concentration=si_nn_doping,
    rmin=(+si_half_w - si_nn_width, y_bot - 1e-4, -z_size),
    rmax=(+si_half_w, y_rib_top + 1e-4, z_size),  # full Si height
)

# --- Build doped semiconductor mediums ---
doped_gesi = intrinsic_gesi.updated_copy(
    N_a=[p_doping_gesi, bg_doping_gesi],
    N_d=[n_doping_gesi],
)
doped_si = intrinsic_si.updated_copy(
    N_a=[p_doping_si, pp_doping_si],
    N_d=[n_doping_si, nn_doping_si],
)

# --- MultiPhysicsMedium ---
GeSi = td.MultiPhysicsMedium(
    charge=doped_gesi, optical=td.material_library["Ge"]["Nunley"], name="GeSi"
)
# same medium with only the SRH lifetimes shortened
tau_dam = td.FossumCarrierLifetime(
    tau_300=tau_damaged, alpha_T=0, A=1, B=0, C=0, N0=1e16, alpha=1
)
GeSiDamaged = td.MultiPhysicsMedium(
    charge=doped_gesi.updated_copy(
        R=[
            r.updated_copy(tau_n=tau_dam, tau_p=tau_dam)
            if isinstance(r, td.ShockleyReedHallRecombination)
            else r
            for r in doped_gesi.R
        ]
    ),
    optical=td.material_library["Ge"]["Nunley"],
    name="GeSi_damaged",
)
Si = td.MultiPhysicsMedium(
    charge=doped_si, optical=td.material_library["cSi"]["Palik_Lossless"], name="Si"
)
SiO2 = td.MultiPhysicsMedium(
    charge=td.ChargeInsulatorMedium(permittivity=3.9),
    optical=td.material_library["SiO2"]["Horiba"],
    name="SiO2",
)
Metal = td.MultiPhysicsMedium(
    charge=td.ChargeConductorMedium(conductivity=1),
    optical=td.Medium(permittivity=1),
    name="Metal",
)

# --- Structure objects ---
# SiO2 background (BOX + cladding)
sio2_struct = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-sim_half_w, -box_h, -z_size),
        rmax=(+sim_half_w, y_gesi_top + clad_h, z_size),
    ),
    medium=SiO2,
    name="sio2_struct",
)
# Si SLAB: thin (140 nm), full lateral width — carries P/N doping on sides
si_slab = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-si_half_w, y_bot, -z_size), rmax=(+si_half_w, y_slab_top, z_size)
    ),
    medium=Si,
    name="si_slab",
)
# Si RIB: full 220 nm, 800 nm wide — overlaps slab in the center, intrinsic (no doping boxes here)
si_rib = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-si_rib_w / 2, y_bot, -z_size), rmax=(+si_rib_w / 2, y_rib_top, z_size)
    ),
    medium=Si,
    name="si_rib",
)
# GeSi waveguide: 600 nm wide, 350 nm tall, partially embedded in the Si rib
gesi_struct = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-gesi_w / 2, y_rib_top / 2, -z_size),
        rmax=(+gesi_w / 2, y_gesi_top, z_size),
    ),
    medium=GeSi,
    name="gesi_struct",
)
# damaged GeSi layer over the first 20 nm above the growth interface
gesi_damaged_struct = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-gesi_w / 2, y_rib_top / 2, -z_size),
        rmax=(+gesi_w / 2, y_rib_top / 2 + damaged_h, z_size),
    ),
    medium=GeSiDamaged,
    name="gesi_damaged",
)
# Metal contacts ON TOP of p++/n++ Si slab regions (above y = y_slab_top)
contact_p = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(-si_half_w, y_slab_top, -z_size),
        rmax=(-si_half_w + si_pp_width, y_slab_top + contact_w, z_size),
    ),
    medium=Metal,
    name="contact_p",
)
contact_n = td.Structure(
    geometry=td.Box.from_bounds(
        rmin=(+si_half_w - si_nn_width, y_slab_top, -z_size),
        rmax=(+si_half_w, y_slab_top + contact_w, z_size),
    ),
    medium=Metal,
    name="contact_n",
)

# Order matters: later structures override earlier ones in overlapping regions
structures = [
    sio2_struct,
    si_slab,
    si_rib,
    gesi_struct,
    gesi_damaged_struct,
    contact_p,
    contact_n,
]

# Build Scene and visualize doping + structure layout
scene = td.Scene(
    structures=structures,
)

fig, axes = plt.subplots(1, 2, figsize=(14, 4))
scene.plot(z=0, ax=axes[0])
axes[0].set_title("Structure cross-section")

scene.plot_structures_property(z=0, property="doping", ax=axes[1], scale="symlog")
axes[1].set_title("Doping map")
plt.tight_layout()
plt.show()

bc_p = td.HeatChargeBoundarySpec(
    placement=td.StructureStructureInterface(structures=[contact_p.name, si_slab.name]),
    condition=td.VoltageBC(source=td.DCVoltageSource(voltage=voltages)),
)
bc_n = td.HeatChargeBoundarySpec(
    placement=td.StructureStructureInterface(structures=[contact_n.name, si_slab.name]),
    condition=td.VoltageBC(source=td.DCVoltageSource(voltage=0.0)),
)
bcs = [bc_p, bc_n]

capacitance_mnt = td.SteadyCapacitanceMonitor(
    center=(0, 0, 0),
    size=(td.inf, td.inf, 0),
    name="capacitance_mnt",
)
# Electric field monitor over the GeSi region — directly provides E = -grad(phi) [V/µm]
efield_mnt = td.SteadyElectricFieldMonitor(
    center=(0, (y_rib_top / 2 + y_gesi_top) / 2, 0),
    size=(gesi_w + 0.2, gesi_h + 0.1, 0),
    name="efield_mnt",
)
monitors = [capacitance_mnt, efield_mnt]

# Mesh: uniform grid inside the semiconductors, refined near material interfaces.
mesh_resolution = 0.01
mesh_spec = td.AutoUnstructuredGrid(
    dl_reference=mesh_resolution,
    uniform_grid_mediums=[GeSi.name, GeSiDamaged.name, Si.name],
)

# the small reverse-bias currents need a tight rel_tol, and at this tolerance the
# linear solver needs more preconditioner_iterations than the default
convergence_settings = td.ChargeToleranceSpec(rel_tol=1e-13, preconditioner_iterations=100)

analysis_type = td.IsothermalSteadyChargeDCAnalysis(
    temperature=300.0,
    convergence_dv=0.5,
    tolerance_settings=convergence_settings,
    fermi_dirac=True,
)

sim = td.HeatChargeSimulation.from_scene(
    scene=scene,
    size=(2 * sim_half_w, box_h + y_gesi_top + clad_h, 0),
    center=(0, (y_gesi_top + clad_h - box_h) / 2, 0),
    monitors=monitors,
    grid_spec=mesh_spec,
    boundary_spec=bcs,
    analysis_spec=analysis_type,
)

# Visualize simulation cross-section (zoomed into the waveguide region)
fig, ax = plt.subplots(figsize=(8, 5))
sim.plot_property(z=0, property="electric_conductivity", ax=ax, monitor_alpha=0)
ax.set_ylim(-0.1, 0.8)
ax.set_title("Simulation cross-section (conductivity map)")
plt.show()

charge_data = web.run(
    sim,
    task_name="GeSiEAM_charge",
    path="gesi_eam_charge_data.hdf5",
)

E_data = charge_data[efield_mnt.name].E  # unstructured; components along 'axis'

# |E| maps on a regular grid over the GeSi cross-section, for a subset of biases
x_map = np.linspace(-0.4, 0.4, 161)
y_map = np.linspace(y_rib_top / 2 - 0.05, y_gesi_top + 0.05, 101)
E_map = E_data.interp(x=x_map, y=y_map, z=0.0)
E_map_mag = np.sqrt(E_map.sel(axis=0) ** 2 + E_map.sel(axis=1) ** 2).squeeze("z") * 10  # kV/cm

# log color scale: the field spans two decades between the intrinsic region and
# the Ge/Si interface
plot_voltages = [-3.0, -1.0, 0.0, 0.75]
fig, axes = plt.subplots(2, 2, figsize=(9, 7))
for ax, v in zip(axes.ravel(), plot_voltages):
    E_map_mag.sel(voltage=v, method="nearest").plot(
        x="x",
        y="y",
        ax=ax,
        cmap="viridis",
        norm=LogNorm(vmin=1, vmax=1e3),
        cbar_kwargs={"label": "|E| [kV/cm]"},
    )
    ax.set_aspect("equal")
    ax.set_title(f"V = {v} V")
plt.suptitle("Electric field in the GeSi region", y=1.02)
plt.tight_layout()
plt.show()

# regular grid over the intrinsic GeSi, full height from the embedded bottom to the top
x_int = np.linspace(-gesi_intrinsic_w / 2, gesi_intrinsic_w / 2, 60)
y_int = np.linspace(y_rib_top / 2, y_gesi_top, 60)

# Tidy3D native interpolation of the unstructured field onto the grid
Eg = E_data.interp(x=x_int, y=y_int, z=0.0)
E_mag = np.sqrt(Eg.sel(axis=0) ** 2 + Eg.sel(axis=1) ** 2).squeeze("z")  # |E|(x,y,voltage)

print(f"{'V_bias (V)':>12} | {'E_avg (kV/cm)':>16}")
print("-" * 32)
for v in voltages:
    e_avg_kv_cm = float(E_mag.sel(voltage=v, method="nearest").mean()) * 1e6 / 1e5  # V/µm -> kV/cm
    print(f"{v:>12.2f} | {e_avg_kv_cm:>16.1f}")

# dark current digitized from Fig. 8(a) of the paper
V_meas = np.array(
    [-3.0, -2.75, -2.5, -2.25, -2.0, -1.75, -1.5, -1.25, -1.0, -0.75, -0.5, -0.25, 0.25, 0.5, 0.75]
)
I_meas_A = np.array(
    [2.23e-7, 1.90e-7, 1.57e-7, 1.22e-7, 9.06e-8, 6.55e-8, 4.46e-8, 2.89e-8,
     1.79e-8, 9.90e-9, 6.14e-9, 4.16e-9, 1.43e-6, 1.01e-4, 7.39e-4]
)

I_per_um = np.abs(charge_data.device_characteristics.steady_dc_current_voltage)  # [A/µm]
V_all = I_per_um.coords["v"].values
I_all_A = I_per_um.values * device_length_um  # [A]
# 0 V is excluded: the net terminal current vanishes at equilibrium
nonzero = V_all != 0.0
V_sim, I_sim_A = V_all[nonzero], I_all_A[nonzero]

fig, ax = plt.subplots(figsize=(6.5, 4.5))
ax.semilogy(V_meas, I_meas_A, "ro", ms=5, mfc="none", label="measured (Fig. 8a)")
ax.semilogy(V_sim, I_sim_A, "ks-", ms=5, label="simulated")
ax.set_xlabel("Applied bias (V)")
ax.set_ylabel("Dark current (A)")
ax.set_title("Dark current vs bias")
ax.set_ylim(1e-10, 1e-3)
ax.legend()
ax.grid(True, which="both", alpha=0.3)
plt.tight_layout()
plt.show()

i_sim_3 = float(I_sim_A[np.isclose(V_sim, -3.0)][0])
i_sim_1 = float(I_sim_A[np.isclose(V_sim, -1.0)][0])
i_meas_3, i_meas_1 = I_meas_A[0], I_meas_A[8]
print(f"I(-3 V): simulated {i_sim_3 * 1e9:.1f} nA, measured {i_meas_3 * 1e9:.1f} nA")
print(f"I(-3)/I(-1): simulated {i_sim_3 / i_sim_1:.2f}, measured {i_meas_3 / i_meas_1:.2f}")

def FranzKeldysh(omega, E_field, n_r, E_p, m_e_Gamma, m_lh, m_hh, E_gd, DeltaE_g_lh, DeltaE_g_hh):
    # E_field     : applied electric field                     [V/µm]

    # Material parameters
    # omega       : optical angular frequency                  [rad/s]
    # n_r         : refractive index                           [dimensionless]
    # E_p         : Kane energy (matrix constant)              [eV]
    # m_e_Gamma   : electron effective mass at the Γ valley    [eV·s²/µm²]  (= (m*/m0)·M_E_EV)
    # m_lh        : light-hole effective mass                  [eV·s²/µm²]
    # m_hh        : heavy-hole effective mass                  [eV·s²/µm²]
    # E_gd        : direct bandgap                            [eV]
    # DeltaE_g_lh : strain-induced gap shift (light hole)      [eV]
    # DeltaE_g_hh : strain-induced gap shift (heavy hole)      [eV]

    # constants
    e = 1  # elementary charge = 1 (eV unit system)
    c = td.C_0  # speed of light [µm/s]
    m_0 = td.constants.M_E_EV  # free-electron mass [eV·s²/µm²]  (= M_E_C_SQUARE / C_0²)
    hbar = td.HBAR  # reduced Planck constant [eV·s]
    eps0 = td.EPSILON_0 / td.Q_e  # F/µm (Coulomb) -> elementary-charge unit (factor 1/e)

    # Airy functions
    from scipy.special import airy

    def Airy(x):
        return airy(x)[0]  # Ai(x)

    def Airy2(x):
        return airy(x)[1]  # Ai'(x)

    # effective mass
    m_r_lh = m_e_Gamma * m_lh / (m_e_Gamma + m_lh)  # reduced e-h mass, light  [eV·s²/µm²]
    m_r_hh = m_e_Gamma * m_hh / (m_e_Gamma + m_hh)  # reduced e-h mass, heavy  [eV·s²/µm²]

    # hbarTheta
    hbarTheta_cte = (hbar**2) * (e**2) * (E_field**2)
    hbarTheta_lh = ((hbarTheta_cte) / (2 * m_r_lh)) ** (1 / 3)  # electro-optic energy ħθ (lh) [eV]
    hbarTheta_hh = ((hbarTheta_cte) / (2 * m_r_hh)) ** (1 / 3)  # electro-optic energy ħθ (hh) [eV]

    # eta
    eta_lh = (E_gd + DeltaE_g_lh - hbar * omega) / hbarTheta_lh  # [dimensionless]
    eta_hh = (E_gd + DeltaE_g_hh - hbar * omega) / hbarTheta_hh  # [dimensionless]

    A = (e**2 * E_p) / (12 * n_r * c * eps0 * m_0 * omega)  # Equation (1) prefactor

    cte_lh = (2 * m_r_lh / hbar**2) ** (3 / 2)
    cte_hh = (2 * m_r_hh / hbar**2) ** (3 / 2)

    B_lh = cte_lh * np.sqrt(hbarTheta_lh) * (-eta_lh * Airy(eta_lh) ** 2 + Airy2(eta_lh) ** 2)
    B_hh = cte_hh * np.sqrt(hbarTheta_hh) * (-eta_hh * Airy(eta_hh) ** 2 + Airy2(eta_hh) ** 2)

    return A * (B_lh + B_hh)  # absorption coefficient α  [1/µm]