// 2D Matrix calculations and Jones matrices for vector optical fields.

///|
/// 2x2 Complex matrix struct
pub struct Matrix2x2 {
  m00 : Complex
  m01 : Complex
  m10 : Complex
  m11 : Complex
} derive(Debug, Eq)

///|
/// Create 2x2 complex matrix
pub fn Matrix2x2::new(
  m00 : Complex,
  m01 : Complex,
  m10 : Complex,
  m11 : Complex,
) -> Matrix2x2 {
  { m00, m01, m10, m11 }
}

///|
/// Identity 2x2 matrix
pub fn Matrix2x2::identity() -> Matrix2x2 {
  {
    m00: Complex::one(),
    m01: Complex::zero(),
    m10: Complex::zero(),
    m11: Complex::one(),
  }
}

///|
/// Create linear polarizer Jones matrix at angle theta (rad)
pub fn jones_polarizer(theta : Double) -> Matrix2x2 {
  let c = @math.cos(theta)
  let s = @math.sin(theta)
  {
    m00: Complex::from_real(c * c),
    m01: Complex::from_real(c * s),
    m10: Complex::from_real(c * s),
    m11: Complex::from_real(s * s),
  }
}

///|
/// Create waveplate (birefringence) Jones matrix with phase retardance phi and fast axis angle theta
pub fn jones_waveplate(phi : Double, theta : Double) -> Matrix2x2 {
  let c = @math.cos(theta)
  let s = @math.sin(theta)
  let exp_pos = Complex::polar(1.0, phi / 2.0)
  let exp_neg = Complex::polar(1.0, -phi / 2.0)

  let a00 = exp_pos.scale(c * c).add(exp_neg.scale(s * s))
  let a01 = exp_pos.sub(exp_neg).scale(c * s)
  let a10 = a01
  let a11 = exp_pos.scale(s * s).add(exp_neg.scale(c * c))

  { m00: a00, m01: a01, m10: a10, m11: a11 }
}

///|
/// Matrix multiplication of two 2x2 Jones matrices
pub fn mul_matrix2x2(a : Matrix2x2, b : Matrix2x2) -> Matrix2x2 {
  let r00 = a.m00.mul(b.m00).add(a.m01.mul(b.m10))
  let r01 = a.m00.mul(b.m01).add(a.m01.mul(b.m11))
  let r10 = a.m10.mul(b.m00).add(a.m11.mul(b.m10))
  let r11 = a.m10.mul(b.m01).add(a.m11.mul(b.m11))
  { m00: r00, m01: r01, m10: r10, m11: r11 }
}

///|
/// Transform 2D complex vector (Ex, Ey) by 2x2 Jones matrix
pub fn apply_jones_matrix(
  m : Matrix2x2,
  ex : Complex,
  ey : Complex,
) -> (Complex, Complex) {
  let new_ex = m.m00.mul(ex).add(m.m01.mul(ey))
  let new_ey = m.m10.mul(ex).add(m.m11.mul(ey))
  (new_ex, new_ey)
}