// Special mathematical functions for fiber optics and mode field distributions.

///|
/// Error function erf(x) approximation via Abramowitz and Stegun formula
pub fn erf(x : Double) -> Double {
  if x == 0.0 {
    return 0.0
  }
  let sign = if x < 0.0 { -1.0 } else { 1.0 }
  let abs_x = if x < 0.0 { -x } else { x }

  let p = 0.3275911
  let a1 = 0.254829592
  let a2 = -0.284496736
  let a3 = 1.421413741
  let a4 = -1.453152027
  let a5 = 1.061405429

  let t = 1.0 / (1.0 + p * abs_x)
  let poly = ((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t
  let res = 1.0 - poly * @math.exp(-abs_x * abs_x)
  sign * res
}

///|
/// Complementary error function erfc(x) = 1 - erf(x)
pub fn erfc(x : Double) -> Double {
  1.0 - erf(x)
}

///|
/// Bessel function of the first kind J0(x) polynomial approximation
pub fn bessel_j0(x : Double) -> Double {
  let abs_x = if x < 0.0 { -x } else { x }
  if abs_x < 3.75 {
    let y = x / 3.75 * (x / 3.75)
    1.0 +
    y *
    (
      -2.2499997 +
      y *
      (
        1.2656208 +
        y * (-0.3163866 + y * (0.0444479 + y * (-0.0039444 + y * 0.0002100)))
      )
    )
  } else {
    let y = 3.75 / abs_x
    let f0 = 0.79788456 +
      y *
      (
        -0.00000077 +
        y *
        (-0.00552740 + y * (0.00009512 + y * (0.00137237 + y * -0.00072805)))
      )
    let theta0 = abs_x -
      0.78539816 +
      y *
      (
        -0.04166397 +
        y *
        (-0.00003954 + y * (0.00262573 + y * (-0.00054125 + y * -0.00029333)))
      )
    f0 / @math.pow(abs_x, 0.5) * @math.cos(theta0)
  }
}

///|
/// Bessel function of the first kind J1(x) polynomial approximation
pub fn bessel_j1(x : Double) -> Double {
  let abs_x = if x < 0.0 { -x } else { x }
  let sign = if x < 0.0 { -1.0 } else { 1.0 }
  if abs_x < 3.75 {
    let y = x / 3.75 * (x / 3.75)
    let ans = x *
      (
        0.5 +
        y *
        (
          -0.56249985 +
          y *
          (0.21093573 + y * (-0.03954289 + y * (0.00443319 + y * -0.00031761)))
        )
      )
    ans
  } else {
    let y = 3.75 / abs_x
    let f1 = 0.79788456 +
      y *
      (
        0.00000156 +
        y * (0.00165966 + y * (0.00017105 + y * (-0.00249511 + y * 0.00113653)))
      )
    let theta1 = abs_x -
      2.35619449 +
      y *
      (
        0.04166397 +
        y *
        (-0.00003954 + y * (-0.00262573 + y * (-0.00054125 + y * 0.00029333)))
      )
    sign * (f1 / @math.pow(abs_x, 0.5)) * @math.cos(theta1)
  }
}

///|
/// Modified Bessel function of the first kind I0(x)
pub fn bessel_i0(x : Double) -> Double {
  let abs_x = if x < 0.0 { -x } else { x }
  if abs_x < 3.75 {
    let y = x / 3.75 * (x / 3.75)
    1.0 +
    y *
    (
      3.5156229 +
      y *
      (
        3.0899424 +
        y * (1.2067492 + y * (0.2659732 + y * (0.0360768 + y * 0.0045813)))
      )
    )
  } else {
    let y = 3.75 / abs_x
    let ans = 0.39894228 +
      y *
      (
        0.01328592 +
        y *
        (0.00225319 + y * (-0.00157565 + y * (0.00916281 + y * -0.02057706)))
      )
    ans / @math.pow(abs_x, 0.5) * @math.exp(abs_x)
  }
}

///|
/// Gamma function Lanczos approximation
pub fn gamma_function(z : Double) -> Double {
  if z < 0.5 {
    PI / (@math.sin(PI * z) * gamma_function(1.0 - z))
  } else {
    let p = [
      676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059,
      12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7,
    ]
    let zm1 = z - 1.0
    let mut x = 0.99999999999980993
    for i = 0; i < p.length(); i = i + 1 {
      x = x + p[i] / (zm1 + i.to_double() + 1.0)
    }
    let t = zm1 + p.length().to_double() - 0.5
    @math.pow(2.0 * PI, 0.5) * @math.pow(t, zm1 + 0.5) * @math.exp(-t) * x
  }
}