//! this example creates a "set/reset" latch from a non-linear ferromagnetic device.
//! this is quite a bit like a "core memory" device.
//! the SR latch in this example is wired to a downstream latch, mostly to show that it's
//! possible to transfer the state (with some limitation) from one latch to another.

use coremem::{Driver, mat, meas};
use coremem::geom::{Coord as _, Meters, Torus};
use coremem::sim::spirv::{SpirvSim, WgpuBackend};
use coremem::sim::units::Seconds;
use coremem::stim::{CurlVectorField, ModulatedVectorField, Sinusoid, TimeVaryingExt as _};


fn main() {
    coremem::init_logging();
    // feature size: the side-length of each discrete grid cell to model (in Meters)
    let feat_size = 20e-6f32;

    // parameters used below to describe the components we construct below. units are (M, A or S).
    let depth = 1600e-6;
    // closest distance between the non-vacuum component and the dissipating boundary
    // longer distances cause boundary reflections to be more dissipated (generally good).
    let buffer_xy = 240e-6;
    // length of our energy-dissipating boundary. longer distances decrease boundary reflections (good)
    let boundary_xy = 500e-6;
    let boundary_z = 300e-6;
    // geometry parameters for the ferrite cores (modeled as torii)
    let ferro_major = 320e-6;
    let ferro_minor = 60e-6;
    let ferro_buffer = 60e-6;
    // geometry parameters for the coupling, drive, and sense wires (modeled as torii)
    let wire_minor = 40e-6;
    let wire_major = 160e-6;
    let wire_coupling_major = 280e-6; // (ferro_minor*4 + ferro_buffer)/2 + wire_minor = 190
    let peak_current = 7.5e6;
    let current_duration = 6.0e-9; // half-wavelength of the sine wave
    let drive_conductivity = 5e6f32;
    let sense_conductivity = 5e3f32;

    let half_depth = depth * 0.5;

    // intermediate computed geometric parameters
    let ferro_top_mid = boundary_xy + buffer_xy + wire_minor + wire_major;
    let ferro_center_y = ferro_top_mid + ferro_major;
    let ferro_bot_mid = ferro_center_y + ferro_major;
    let height = ferro_bot_mid + wire_major + wire_minor + buffer_xy + boundary_xy;
    let ferro1_left_edge = boundary_xy + buffer_xy;
    let ferro1_center = ferro1_left_edge + ferro_minor + ferro_major;
    let ferro1_right_edge = ferro1_center + ferro_major + ferro_minor;
    let ferro2_left_edge = ferro1_right_edge + ferro_buffer;
    let ferro2_center = ferro2_left_edge + ferro_minor + ferro_major;
    let ferro2_right_edge = ferro2_center + ferro_major + ferro_minor;
    let width = ferro2_right_edge + wire_major + wire_minor + buffer_xy + boundary_xy;

    // create actual Regions from the computed parameters
    let ferro1_region = Torus::new_xy(Meters::new(ferro1_center, ferro_center_y, half_depth), ferro_major, ferro_minor);
    let ferro2_region = Torus::new_xy(Meters::new(ferro2_center, ferro_center_y, half_depth), ferro_major, ferro_minor);
    let set_region = Torus::new_yz(Meters::new(ferro1_center, ferro_center_y - ferro_major, half_depth), wire_major, wire_minor);
    let reset_region = Torus::new_yz(Meters::new(ferro1_center, ferro_center_y + ferro_major, half_depth), wire_major, wire_minor);
    let coupling_region = Torus::new_xz(Meters::new(0.5*(ferro1_center + ferro2_center), ferro_center_y, half_depth), wire_coupling_major, wire_minor);
    let sense_region = Torus::new_xz(Meters::new(ferro2_center + ferro_major, ferro_center_y, half_depth), wire_major, wire_minor);

    let mut driver = Driver::new(SpirvSim::<f32, mat::GenericMaterial<f32>, WgpuBackend>::new(
        Meters::new(width, height, depth).to_index(feat_size), feat_size
    ));

    // mu_r=881.33, starting at H=25 to H=75.
    driver.fill_region(&ferro1_region, mat::MHPgram::new(25.0, 881.33, 44000.0));
    driver.fill_region(&ferro2_region, mat::MHPgram::new(25.0, 881.33, 44000.0));
    driver.fill_region(&set_region, mat::IsomorphicConductor::new(drive_conductivity));
    driver.fill_region(&reset_region, mat::IsomorphicConductor::new(drive_conductivity));
    driver.fill_region(&coupling_region, mat::IsomorphicConductor::new(drive_conductivity));
    driver.fill_region(&sense_region, mat::IsomorphicConductor::new(sense_conductivity));

    // fill the edge of the simulation with a graded conductor:
    // on the inside of the simulation it matches vacuum (conductivity of 0), and
    // ramps up conductivity (to dissipate energy) as it approaches the edge of the simulation.
    driver.add_classical_boundary(Meters::new(boundary_xy, boundary_xy, boundary_z));

    // helper to schedule a stimulus at the provided start time/duration.
    let mut add_drive_pulse = |region: &Torus, start, duration, amp| {
        let wave = Sinusoid::from_wavelength(duration * 2.0)
            .half_cycle()
            .scaled(amp)
            .shifted(start);
        driver.add_stimulus(ModulatedVectorField::new(
            CurlVectorField::new(region.clone()),
            wave,
        ));
    };

    // stimuli apply some delta_E to the simulation, so we need to map our current to E:
    // J=\sigma E
    // dJ/dt = \sigma dE/dT
    // dE/dt = dJ/dt / \sigma
    // dE/dt = dI/dt / (A*\sigma)
    // if I = k*sin(w t) then dE/dt = k*w cos(w t) / (A*\sigma)
    // i.e. dE/dt is proportional to I/(A*\sigma), multiplied by w (or, divided by wavelength)
    let peak_stim1 = peak_current/current_duration / (set_region.cross_section() * drive_conductivity);

    // pulse the SET wire near the start of the simulation
    add_drive_pulse(&set_region, 0.01*current_duration, current_duration, peak_stim1);

    // RESET
    add_drive_pulse(&reset_region, 4.0*current_duration, current_duration, peak_stim1);

    // TOGGLE
    add_drive_pulse(&set_region, 8.0*current_duration, current_duration, peak_stim1);
    add_drive_pulse(&reset_region, 10.0*current_duration, current_duration, peak_stim1);

    // TOGETHER
    add_drive_pulse(&reset_region, 14.0*current_duration, current_duration, peak_stim1);
    add_drive_pulse(&set_region, 14.0*current_duration, current_duration, peak_stim1);

    // SET TWICE
    add_drive_pulse(&set_region, 18.0*current_duration, current_duration, peak_stim1);
    add_drive_pulse(&set_region, 20.0*current_duration, current_duration, peak_stim1);

    let duration = 25.0*current_duration;

    // measure a bunch of items of interest throughout the whole simulation duration:
    driver.add_measurement(meas::CurrentLoop::new("coupling", coupling_region.clone()));
    driver.add_measurement(meas::Current::new("coupling", coupling_region.clone()));
    driver.add_measurement(meas::CurrentLoop::new("sense", sense_region.clone()));
    driver.add_measurement(meas::Current::new("sense", sense_region.clone()));
    driver.add_measurement(meas::MagneticLoop::new("mem1", ferro1_region.clone()));
    driver.add_measurement(meas::Magnetization::new("mem1", ferro1_region.clone()));
    driver.add_measurement(meas::MagneticFlux::new("mem1", ferro1_region.clone()));
    driver.add_measurement(meas::MagneticLoop::new("mem2", ferro2_region.clone()));
    driver.add_measurement(meas::CurrentLoop::new("set", set_region.clone()));
    driver.add_measurement(meas::Current::new("set", set_region.clone()));
    driver.add_measurement(meas::Power::new("set", set_region.clone()));
    driver.add_measurement(meas::CurrentLoop::new("reset", reset_region.clone()));

    // XXX: if you change any parameters (above), then change this prefix. otherwise simulations
    // will try to load state generated by an earlier run and use it to compute the current
    // (differently-parameterized) run.
    let prefix = "out/applications/sr-latch/";
    let _ = std::fs::create_dir_all(&prefix);

    // add a state file for easy resumption
    driver.add_state_file(&*format!("{}state.bc", prefix), 9600);
    // serialize frames for later viewing with `cargo run --release --bin viewer`
    driver.add_serializer_renderer(&*format!("{}frame-", prefix), 36000, None);
    // render a couple CSV files: one very detailed and the other more sparsely detailed
    driver.add_csv_renderer(&*format!("{}meas.csv", prefix), 200, None);
    driver.add_csv_renderer(&*format!("{}meas-sparse.csv", prefix), 1600, None);

    driver.step_until(Seconds(duration));
}