#![feature(inclusive_range_syntax)]

extern crate rulinalg;
extern crate electro;

use rulinalg::matrix::BaseMatrix;
use electro::fd::FiniteDifference;
use electro::ElectroGrid;

fn main() {
  // create potential matrix for input to finite diff function

  // grid parameters
  let grid_size = 1;
  let rect_dims: (usize, usize) = (10, 10);
  let row_count = rect_dims.0 / grid_size - 1;
  let col_count = rect_dims.1 / grid_size - 1;

  // rectangular boundary conditions in ESWN order
  let edge_pot: (f64, f64, f64, f64) = (0.0, 0.0, 0.0, 0.0);

  // 2. create field grid matrix (i.e., "physical" mapping)
  let mut grid =
    ElectroGrid::<f64,f64>::uniform_freespace_zero(row_count as f64,
                                                   col_count as f64,
                                                   1.0);
  // Equivalent to
  //
  // let mut grid = ElectroGrid::<f64> {
  //   space: uniform_grid(1.0, row_count as f64, col_count as f64),
  //   spacing: [1,1],
  //   potential: Matrix::<f64>::zeros(row_count, col_count),
  //   permittivity: Matrix::<f64>::ones(row_count, col_count)
  // };
  
  // set boundary values
  for row_ind in 0..row_count {
    grid.potential[[row_ind, col_count - 1]] = edge_pot.0;
    grid.potential[[row_ind, 0]] = edge_pot.2;
  }
  for col_ind in 0..col_count {
    grid.potential[[row_count - 1, col_ind]] = edge_pot.1;
    grid.potential[[0, col_ind]] = edge_pot.3;
  }
  // 2a. Set corners to averaged corner boundaries if both are nonzero
  if edge_pot.0 != 0.0 && edge_pot.1 != 0.0 && edge_pot.0 != edge_pot.1 { 
    grid.potential[[row_count - 1, col_count - 1]] = (edge_pot.0 + edge_pot.1) / 2.0;
  }
  if edge_pot.2 != 0.0 && edge_pot.2 != 0.0 && edge_pot.2 != edge_pot.3 { 
    grid.potential[[0, 0]] = (edge_pot.2 + edge_pot.3) / 2.0;
  }
  if edge_pot.0 != 0.0 && edge_pot.3 != 0.0 && edge_pot.0 != edge_pot.3 { 
    grid.potential[[0, col_count - 1]] = (edge_pot.0 + edge_pot.3) / 2.0;
  }
  if edge_pot.1 != 0.0 && edge_pot.2 != 0.0 && edge_pot.1 != edge_pot.2 { 
    grid.potential[[row_count - 1, 0]] = (edge_pot.1 + edge_pot.2) / 2.0;
  }

  // create center conductor (single row)
  
  // conductor params
  let conductor_voltage = 100.0;
  let conductor_width = 4;
  let conductor_col_width = conductor_width / grid_size;
  let conductor_row = row_count / 2;
  let middle_col = col_count / 2;
  let conductor_start = middle_col - conductor_col_width / 2;
  let conductor_end = middle_col + conductor_col_width / 2;
  let conductor_col_inds = conductor_start...conductor_end;

  for col_ind in conductor_col_inds {
    grid.potential[[conductor_row, col_ind]] = conductor_voltage;
  }

  let computed_grid = grid.finite_diff_rect().unwrap();

  for row in computed_grid.potential.iter_rows() {
    println!("{:?}", row);
    println!("");
  }

  
  // Compute the capacitance by evaluating a contour integral around the
  // conductive element
  let conductor_region: [[usize; 2]; 4] = [
    [conductor_row + 2, conductor_start - 2],
    [conductor_row + 2, conductor_end + 2],
    [conductor_row - 2, conductor_end + 2],
    [conductor_row - 2, conductor_start - 2]
  ];

  let capacitance = 0.5 / conductor_voltage *
    computed_grid.contour_int_rect(&conductor_region, false);

  println!("region: {:?}", &conductor_region);
  println!("capacitance: {}", capacitance);
}