oxiphoton 0.1.1

//! 3D FDTD waveguide propagation example.
//!
//! Simulates a TE-mode Gaussian pulse propagating through a silicon-on-insulator
//! strip waveguide using the OxiPhoton 3D FDTD engine with CPML absorbing boundaries.
//!
//! Waveguide geometry:
//!   - Si core: 400 nm × 220 nm, n_core = 3.476 (ε_r = 12.08)
//!   - SiO₂ substrate: n_sub = 1.444 (ε_r = 2.09)
//!   - Air cladding: n_air = 1.0
//!   - Operating wavelength: 1550 nm
//!
//! Simulation parameters:
//!   - Domain: 3.2 μm × 2.4 μm × 6.0 μm (x × y × z), dx = 40 nm
//!   - CPML: 8 cells thick on all boundaries
//!   - Time step: set by Courant condition
//!   - Source: Gaussian pulse injected into waveguide cross-section
//!   - Duration: 400 time steps
//!
//! Run with: cargo run --example fdtd_3d_waveguide --all-features

use std::f64::consts::PI;

use oxiphoton::fdtd::{BoundaryConfig, Fdtd3d};

fn main() {
    println!("=== 3D FDTD Si Waveguide Propagation Example ===");
    println!();

    // ── Grid parameters ──────────────────────────────────────────────────────
    let dx = 40e-9; // 40 nm cell size
    let lambda = 1550e-9; // target wavelength
    let n_max = 3.476_f64; // maximum refractive index (Si)
    let c = 2.997_924_58e8_f64;

    // Grid dimensions (small for fast demo, but physically meaningful)
    let nx = 80usize; // 80 × 40nm = 3.2 μm
    let ny = 60usize; // 60 × 40nm = 2.4 μm
    let nz = 150usize; // 150 × 40nm = 6.0 μm

    let n_steps = 400;
    let pml_cells = 8usize;

    // BoundaryConfig determines PML
    let boundary = BoundaryConfig::pml(pml_cells);

    // ── Build FDTD solver ────────────────────────────────────────────────────
    let mut fdtd = Fdtd3d::new(nx, ny, nz, dx, dx, dx, &boundary);
    let dt = fdtd.dt;

    // Courant number (computed after building solver)
    let courant_number = c * dt / dx * 3.0_f64.sqrt();

    println!("Grid:          {}×{}×{} cells", nx, ny, nz);
    println!("Cell size:     {:.0} nm", dx * 1e9);
    println!(
        "Domain:        {:.1} × {:.1} × {:.1} μm",
        nx as f64 * dx * 1e6,
        ny as f64 * dx * 1e6,
        nz as f64 * dx * 1e6
    );
    println!("Time step:     {:.4e} s", dt);
    println!("Steps:         {}", n_steps);
    println!("CPML cells:    {}", pml_cells);
    println!();
    println!("Courant number (×n_max): {:.4}", courant_number * n_max);
    println!("Courant stable:          {}", courant_number * n_max < 1.0);
    println!();

    // ── Define materials ─────────────────────────────────────────────────────
    // SiO₂ substrate: lower half (y < ny/3)
    let eps_sio2 = 2.09_f64;
    let y_sub_top = ny / 3;
    fdtd.fill_box(0, nx, 0, y_sub_top, 0, nz, eps_sio2, 1.0);

    // Si waveguide core: 400nm × 220nm strip at center of domain
    let eps_si = 12.08_f64; // n_Si = 3.476 → ε = n²
                            // Convert physical dimensions to grid cells
    let wg_width_cells = ((400e-9) / dx).round() as usize; // 10 cells
    let wg_height_cells = ((220e-9) / dx).round() as usize; // 6 cells

    let wg_x_start = (nx / 2).saturating_sub(wg_width_cells / 2);
    let wg_x_end = wg_x_start + wg_width_cells;
    let wg_y_start = y_sub_top;
    let wg_y_end = y_sub_top + wg_height_cells;

    fdtd.fill_box(
        wg_x_start,
        wg_x_end,
        wg_y_start,
        wg_y_end,
        pml_cells,
        nz - pml_cells,
        eps_si,
        1.0,
    );

    println!("Material regions:");
    println!(
        "  SiO₂ substrate (ε={:.2}): y = 0..{} cells",
        eps_sio2, y_sub_top
    );
    println!(
        "  Si waveguide (ε={:.2}):  x={}..{}, y={}..{} cells",
        eps_si, wg_x_start, wg_x_end, wg_y_start, wg_y_end
    );
    println!("  Air cladding (ε=1.0):   y > {} cells", wg_y_end);
    println!();

    // ── Add field probe ──────────────────────────────────────────────────────
    // Record Ez at the waveguide center, near the output face
    let probe_i = nx / 2;
    let probe_j = wg_y_start + wg_height_cells / 2;
    let probe_k = nz - pml_cells - 5; // 5 cells from exit CPML

    // ── Source: Gaussian pulse injected into waveguide ───────────────────────
    // Source position: near input face (k = pml_cells + 5)
    let src_k = pml_cells + 5;
    let t0 = 50.0 * dt; // pulse center time
    let sigma = 20.0 * dt; // pulse width
    let f0_hz = c / lambda; // center frequency (193 THz)
    let amp = 1.0_f64; // field amplitude (V/m normalized)

    // ── Run simulation ───────────────────────────────────────────────────────
    println!("Running {} time steps...", n_steps);
    let t_start = std::time::Instant::now();

    let mut probe_ez = Vec::with_capacity(n_steps);
    let mut energy_vs_step = Vec::with_capacity(n_steps / 10);

    for step in 0..n_steps {
        let t = step as f64 * dt;

        // Gaussian pulse source: envelope × carrier
        let envelope = (-(t - t0).powi(2) / (2.0 * sigma * sigma)).exp();
        let carrier = (2.0 * PI * f0_hz * t).sin();
        let src_amp = amp * envelope * carrier;

        // Inject across the waveguide cross-section for better mode overlap
        for j in wg_y_start..wg_y_end {
            for i in wg_x_start..wg_x_end {
                // Spatial mode profile: 2D Gaussian approximation
                let xi = (i as f64 - (wg_x_start + wg_x_end) as f64 / 2.0)
                    / (wg_width_cells as f64 / 2.0);
                let yj = (j as f64 - (wg_y_start + wg_y_end) as f64 / 2.0)
                    / (wg_height_cells as f64 / 2.0);
                let spatial = (-xi * xi - yj * yj).exp();
                fdtd.inject_ez(i, j, src_k, src_amp * spatial);
            }
        }

        // Record probe value (Ez is the 3rd field component in field_at output [ex,ey,ez,hx,hy,hz])
        let fields = fdtd.field_at(probe_i, probe_j, probe_k);
        probe_ez.push(fields[2]); // Ez

        // Record energy every 10 steps
        if step % 10 == 0 {
            energy_vs_step.push((step, fdtd.total_energy()));
        }

        fdtd.step();
    }

    let elapsed = t_start.elapsed();
    println!("Simulation complete in {:.2} s", elapsed.as_secs_f64());
    println!();

    // ── Analysis and reporting ───────────────────────────────────────────────

    // Peak field at probe
    let peak_ez = probe_ez
        .iter()
        .cloned()
        .fold(0.0_f64, |a, v| a.max(v.abs()));
    let rms_ez = (probe_ez.iter().map(|&v| v * v).sum::<f64>() / probe_ez.len() as f64).sqrt();

    println!("=== Field Analysis at Output Probe ===");
    println!(
        "Probe location: ({}, {}, {}) cells = ({:.1}, {:.1}, {:.1}) μm",
        probe_i,
        probe_j,
        probe_k,
        probe_i as f64 * dx * 1e6,
        probe_j as f64 * dx * 1e6,
        probe_k as f64 * dx * 1e6,
    );
    println!("Peak |Ez|:   {:.6e} (normalized)", peak_ez);
    println!("RMS Ez:      {:.6e} (normalized)", rms_ez);
    println!();

    // Energy evolution
    let e_initial = energy_vs_step.first().map(|&(_, e)| e).unwrap_or(0.0);
    let e_peak = energy_vs_step
        .iter()
        .map(|&(_, e)| e)
        .fold(0.0_f64, f64::max);
    let e_final = energy_vs_step.last().map(|&(_, e)| e).unwrap_or(0.0);

    println!("=== Energy Analysis ===");
    println!("Initial energy: {:.4e} J/m (normalized)", e_initial);
    println!("Peak energy:    {:.4e} J/m (normalized)", e_peak);
    println!("Final energy:   {:.4e} J/m (normalized)", e_final);
    println!();

    // Physical parameters summary
    let wavelength_cells = lambda / dx;
    println!("=== Physical Parameters ===");
    println!(
        "Wavelength:            {:.0} nm ({:.1} cells)",
        lambda * 1e9,
        wavelength_cells
    );
    println!("n_Si (core):           {:.4}", n_max);
    println!("n_SiO2 (substrate):    {:.4}", eps_sio2.sqrt());
    println!("λ_eff in Si (approx):  {:.0} nm", lambda / n_max * 1e9);
    println!("Courant number:        {:.4}", courant_number);
    println!(
        "Simulated time:        {:.4} ps",
        n_steps as f64 * dt * 1e12
    );
    println!();

    // Propagation length
    let prop_length_um = (probe_k - src_k) as f64 * dx * 1e6;
    println!("=== Propagation Analysis ===");
    println!("Source-to-probe distance: {:.2} μm", prop_length_um);
    println!("Approximate propagation loss:  < 0.1 dB/μm (ideal FDTD)");

    if peak_ez > 1e-15 {
        println!("✓ Waveguide mode successfully excited and propagated");
    } else {
        println!("⚠ Low probe field — check source/probe alignment");
    }

    println!();
    println!("=== Summary ===");
    println!("OxiPhoton 3D FDTD Si waveguide simulation complete.");
    println!(
        "Grid: {}×{}×{}, {:.0}nm resolution, {} steps",
        nx,
        ny,
        nz,
        dx * 1e9,
        n_steps
    );
    println!(
        "Estimated memory: ~{:.1} MB",
        3 * nx * ny * nz * 8 / (1024 * 1024)
    );
}