// Complex numbers for numerical optics calculations.

///|
pub struct Complex {
  re : Double
  im : Double
} derive(Debug, Eq)

///|
pub fn Complex::new(re : Double, im : Double) -> Complex {
  { re, im }
}

///|
pub fn Complex::zero() -> Complex {
  { re: 0.0, im: 0.0 }
}

///|
pub fn Complex::one() -> Complex {
  { re: 1.0, im: 0.0 }
}

///|
pub fn Complex::i() -> Complex {
  { re: 0.0, im: 1.0 }
}

///|
pub fn Complex::from_real(r : Double) -> Complex {
  { re: r, im: 0.0 }
}

///|
pub fn Complex::from_imag(i : Double) -> Complex {
  { re: 0.0, im: i }
}

///|
pub fn Complex::add(self : Complex, other : Complex) -> Complex {
  { re: self.re + other.re, im: self.im + other.im }
}

///|
pub fn Complex::op_add(self : Complex, other : Complex) -> Complex {
  self.add(other)
}

///|
pub fn Complex::sub(self : Complex, other : Complex) -> Complex {
  { re: self.re - other.re, im: self.im - other.im }
}

///|
pub fn Complex::op_sub(self : Complex, other : Complex) -> Complex {
  self.sub(other)
}

///|
pub fn Complex::mul(self : Complex, other : Complex) -> Complex {
  {
    re: self.re * other.re - self.im * other.im,
    im: self.re * other.im + self.im * other.re,
  }
}

///|
pub fn Complex::op_mul(self : Complex, other : Complex) -> Complex {
  self.mul(other)
}

///|
pub fn Complex::scale(self : Complex, s : Double) -> Complex {
  { re: self.re * s, im: self.im * s }
}

///|
pub fn Complex::div(self : Complex, other : Complex) -> Complex {
  let denom = other.re * other.re + other.im * other.im
  if denom == 0.0 {
    { re: 0.0, im: 0.0 }
  } else {
    {
      re: (self.re * other.re + self.im * other.im) / denom,
      im: (self.im * other.re - self.re * other.im) / denom,
    }
  }
}

///|
pub fn Complex::op_div(self : Complex, other : Complex) -> Complex {
  self.div(other)
}

///|
pub fn Complex::div_real(self : Complex, s : Double) -> Complex {
  if s == 0.0 {
    { re: 0.0, im: 0.0 }
  } else {
    { re: self.re / s, im: self.im / s }
  }
}

///|
pub fn Complex::neg(self : Complex) -> Complex {
  { re: -self.re, im: -self.im }
}

///|
pub fn Complex::conj(self : Complex) -> Complex {
  { re: self.re, im: -self.im }
}

///|
pub fn Complex::abs2(self : Complex) -> Complex {
  { re: self.re * self.re + self.im * self.im, im: 0.0 }
}

///|
pub fn Complex::magnitude2(self : Complex) -> Double {
  self.re * self.re + self.im * self.im
}

///|
pub fn Complex::abs(self : Complex) -> Double {
  let a2 = self.magnitude2()
  if a2 <= 0.0 {
    0.0
  } else {
    @math.pow(a2, 0.5)
  }
}

///|
pub fn Complex::arg(self : Complex) -> Double {
  let y = self.im
  let x = self.re
  if x == 0.0 && y == 0.0 {
    0.0
  } else {
    @math.atan2(y, x)
  }
}

///|
pub fn Complex::polar(r : Double, theta : Double) -> Complex {
  { re: r * @math.cos(theta), im: r * @math.sin(theta) }
}

///|
pub fn Complex::exp(self : Complex) -> Complex {
  let ex = @math.exp(self.re)
  { re: ex * @math.cos(self.im), im: ex * @math.sin(self.im) }
}

///|
pub fn Complex::log(self : Complex) -> Complex {
  let r = self.abs()
  let theta = self.arg()
  if r <= 0.0 {
    { re: 0.0, im: 0.0 }
  } else {
    { re: @math.ln(r), im: theta }
  }
}

///|
pub fn Complex::pow_real(self : Complex, p : Double) -> Complex {
  let r = self.abs()
  let theta = self.arg()
  if r <= 0.0 {
    Complex::zero()
  } else {
    let new_r = @math.pow(r, p)
    Complex::polar(new_r, theta * p)
  }
}

///|
pub fn Complex::sqrt(self : Complex) -> Complex {
  self.pow_real(0.5)
}

///|
pub fn Complex::sin(self : Complex) -> Complex {
  let x = self.re
  let y = self.im
  let cosh_y = (@math.exp(y) + @math.exp(-y)) / 2.0
  let sinh_y = (@math.exp(y) - @math.exp(-y)) / 2.0
  { re: @math.sin(x) * cosh_y, im: @math.cos(x) * sinh_y }
}

///|
pub fn Complex::cos(self : Complex) -> Complex {
  let x = self.re
  let y = self.im
  let cosh_y = (@math.exp(y) + @math.exp(-y)) / 2.0
  let sinh_y = (@math.exp(y) - @math.exp(-y)) / 2.0
  { re: @math.cos(x) * cosh_y, im: -@math.sin(x) * sinh_y }
}

///|
pub fn Complex::tan(self : Complex) -> Complex {
  self.sin().div(self.cos())
}

///|
pub fn Complex::sinh(self : Complex) -> Complex {
  let ez = self.exp()
  let emz = self.neg().exp()
  ez.sub(emz).scale(0.5)
}

///|
pub fn Complex::cosh(self : Complex) -> Complex {
  let ez = self.exp()
  let emz = self.neg().exp()
  ez.add(emz).scale(0.5)
}

///|
pub fn Complex::tanh(self : Complex) -> Complex {
  self.sinh().div(self.cosh())
}

///|
pub fn Complex::sech(self : Complex) -> Complex {
  Complex::one().div(self.cosh())
}