// Array and matrix math utilities for optics fields.

///|
pub fn linspace(start : Double, stop : Double, num : Int) -> Array[Double] {
  if num <= 0 {
    return []
  }
  if num == 1 {
    return [start]
  }
  let res = Array::make(num, 0.0)
  let step = (stop - start) / (num.to_double() - 1.0)
  for i = 0; i < num; i = i + 1 {
    res[i] = start + i.to_double() * step
  }
  res
}

///|
pub fn arange(start : Double, stop : Double, step : Double) -> Array[Double] {
  if step == 0.0 || (stop - start) * step < 0.0 {
    return []
  }
  let count = ((stop - start) / step).to_int()
  let res = Array::make(count, 0.0)
  for i = 0; i < count; i = i + 1 {
    res[i] = start + i.to_double() * step
  }
  res
}

///|
pub fn zeros_complex(size : Int) -> Array[Complex] {
  Array::make(size, Complex::zero())
}

///|
pub fn zeros_double(size : Int) -> Array[Double] {
  Array::make(size, 0.0)
}

///|
pub fn copy_complex(arr : Array[Complex]) -> Array[Complex] {
  let len = arr.length()
  let res = Array::make(len, Complex::zero())
  for i = 0; i < len; i = i + 1 {
    res[i] = arr[i]
  }
  res
}

///|
pub fn copy_double(arr : Array[Double]) -> Array[Double] {
  let len = arr.length()
  let res = Array::make(len, 0.0)
  for i = 0; i < len; i = i + 1 {
    res[i] = arr[i]
  }
  res
}

///|
pub fn complex_magnitude_squared(arr : Array[Complex]) -> Array[Double] {
  let len = arr.length()
  let res = Array::make(len, 0.0)
  for i = 0; i < len; i = i + 1 {
    res[i] = arr[i].magnitude2()
  }
  res
}

///|
pub fn complex_magnitude(arr : Array[Complex]) -> Array[Double] {
  let len = arr.length()
  let res = Array::make(len, 0.0)
  for i = 0; i < len; i = i + 1 {
    res[i] = arr[i].abs()
  }
  res
}

///|
pub fn complex_phase(arr : Array[Complex]) -> Array[Double] {
  let len = arr.length()
  let res = Array::make(len, 0.0)
  for i = 0; i < len; i = i + 1 {
    res[i] = arr[i].arg()
  }
  res
}

///|
pub fn array_max_double(arr : Array[Double]) -> Double {
  if arr.length() == 0 {
    return 0.0
  }
  let mut m = arr[0]
  for i = 1; i < arr.length(); i = i + 1 {
    if arr[i] > m {
      m = arr[i]
    }
  }
  m
}

///|
pub fn array_min_double(arr : Array[Double]) -> Double {
  if arr.length() == 0 {
    return 0.0
  }
  let mut m = arr[0]
  for i = 1; i < arr.length(); i = i + 1 {
    if arr[i] < m {
      m = arr[i]
    }
  }
  m
}

///|
pub fn array_sum_double(arr : Array[Double]) -> Double {
  let mut s = 0.0
  for i = 0; i < arr.length(); i = i + 1 {
    s = s + arr[i]
  }
  s
}

///|
pub fn trapz_double(y : Array[Double], dx : Double) -> Double {
  let len = y.length()
  if len < 2 {
    return 0.0
  }
  let mut sum = (y[0] + y[len - 1]) / 2.0
  for i = 1; i < len - 1; i = i + 1 {
    sum = sum + y[i]
  }
  sum * dx
}

///|
pub fn trapz_complex_power(field : Array[Complex], dt : Double) -> Double {
  let powers = complex_magnitude_squared(field)
  trapz_double(powers, dt)
}

///|
pub fn scale_complex_array(arr : Array[Complex], factor : Double) -> Unit {
  for i = 0; i < arr.length(); i = i + 1 {
    arr[i] = arr[i].scale(factor)
  }
}

///|
pub fn add_complex_arrays(
  a : Array[Complex],
  b : Array[Complex],
) -> Array[Complex] {
  let len = a.length()
  let res = Array::make(len, Complex::zero())
  for i = 0; i < len; i = i + 1 {
    res[i] = a[i].add(b[i])
  }
  res
}

///|
pub fn sub_complex_arrays(
  a : Array[Complex],
  b : Array[Complex],
) -> Array[Complex] {
  let len = a.length()
  let res = Array::make(len, Complex::zero())
  for i = 0; i < len; i = i + 1 {
    res[i] = a[i].sub(b[i])
  }
  res
}

///|
pub fn mul_complex_arrays(
  a : Array[Complex],
  b : Array[Complex],
) -> Array[Complex] {
  let len = a.length()
  let res = Array::make(len, Complex::zero())
  for i = 0; i < len; i = i + 1 {
    res[i] = a[i].mul(b[i])
  }
  res
}