// LU and Cholesky decomposition for 3x3 and 4x4 optical Jones matrices.

///|
/// 3x3 Double matrix container
pub struct Matrix3x3 {
  m00 : Double
  m01 : Double
  m02 : Double
  m10 : Double
  m11 : Double
  m12 : Double
  m20 : Double
  m21 : Double
  m22 : Double
} derive(Debug, Eq)

///|
/// Create 3x3 Double matrix
pub fn Matrix3x3::new(
  m00 : Double,
  m01 : Double,
  m02 : Double,
  m10 : Double,
  m11 : Double,
  m12 : Double,
  m20 : Double,
  m21 : Double,
  m22 : Double,
) -> Matrix3x3 {
  { m00, m01, m02, m10, m11, m12, m20, m21, m22 }
}

///|
/// Identity 3x3 Double matrix
pub fn Matrix3x3::identity() -> Matrix3x3 {
  {
    m00: 1.0,
    m01: 0.0,
    m02: 0.0,
    m10: 0.0,
    m11: 1.0,
    m12: 0.0,
    m20: 0.0,
    m21: 0.0,
    m22: 1.0,
  }
}

///|
/// Compute determinant of 3x3 matrix
pub fn Matrix3x3::det(self : Matrix3x3) -> Double {
  self.m00 * (self.m11 * self.m22 - self.m12 * self.m21) -
  self.m01 * (self.m10 * self.m22 - self.m12 * self.m20) +
  self.m02 * (self.m10 * self.m21 - self.m11 * self.m20)
}

///|
/// Inverse of 3x3 matrix
pub fn Matrix3x3::inverse(self : Matrix3x3) -> Matrix3x3 {
  let d = self.det()
  if d == 0.0 {
    Matrix3x3::identity()
  } else {
    let inv_d = 1.0 / d
    {
      m00: (self.m11 * self.m22 - self.m12 * self.m21) * inv_d,
      m01: (self.m02 * self.m21 - self.m01 * self.m22) * inv_d,
      m02: (self.m01 * self.m12 - self.m02 * self.m11) * inv_d,
      m10: (self.m12 * self.m20 - self.m10 * self.m22) * inv_d,
      m11: (self.m00 * self.m22 - self.m02 * self.m20) * inv_d,
      m12: (self.m02 * self.m10 - self.m00 * self.m12) * inv_d,
      m20: (self.m10 * self.m21 - self.m11 * self.m20) * inv_d,
      m21: (self.m01 * self.m20 - self.m00 * self.m21) * inv_d,
      m22: (self.m00 * self.m11 - self.m01 * self.m10) * inv_d,
    }
  }
}

///|
/// 4x4 Complex matrix for multimode spatial-polarization optics
pub struct Matrix4x4Complex {
  m00 : Complex
  m01 : Complex
  m02 : Complex
  m03 : Complex
  m10 : Complex
  m11 : Complex
  m12 : Complex
  m13 : Complex
  m20 : Complex
  m21 : Complex
  m22 : Complex
  m23 : Complex
  m30 : Complex
  m31 : Complex
  m32 : Complex
  m33 : Complex
} derive(Debug, Eq)

///|
/// Identity 4x4 Complex matrix
pub fn Matrix4x4Complex::identity() -> Matrix4x4Complex {
  {
    m00: Complex::one(),
    m01: Complex::zero(),
    m02: Complex::zero(),
    m03: Complex::zero(),
    m10: Complex::zero(),
    m11: Complex::one(),
    m12: Complex::zero(),
    m13: Complex::zero(),
    m20: Complex::zero(),
    m21: Complex::zero(),
    m22: Complex::one(),
    m23: Complex::zero(),
    m30: Complex::zero(),
    m31: Complex::zero(),
    m32: Complex::zero(),
    m33: Complex::one(),
  }
}