BESSELY

63 bessely "BESSELY" (2, 2) [V "x", V "n"] |s, a, sv| { bessel(s, a, k::bessel::bessely_kernel, false) }

fn y0(x: D) -> R {
    admit_phase(x, 0.785398164)?;
    if x.0 >= 8.0 {
        let value = super::bessel_rational::y0_large(x.0);
        return if value.is_finite() { Ok(D::new(value)) } else { Err(WorksheetErrorCode::Num) };
    }

    if x.0 > 0.0 && x.0 < 8.0 {
        let value = super::bessel_rational::y0(x.0);
        return if value.is_finite() { Ok(D::new(value)) } else { Err(WorksheetErrorCode::Num) };
    }
 Ok(D::new(libm::y0(x.0))) }

fn y1(x: D) -> R {
    admit_phase(x, 2.356194491)?;
    if x.0 >= 8.0 {
        let value = super::bessel_rational::y1_large(x.0);
        return if value.is_finite() { Ok(D::new(value)) } else { Err(WorksheetErrorCode::Num) };
    }

    if x.0 > 0.0 && x.0 < 8.0 {
        let value = super::bessel_rational::y1(x.0);
        return if value.is_finite() { Ok(D::new(value)) } else { Err(WorksheetErrorCode::Num) };
    }
 Ok(D::new(libm::y1(x.0))) }

/// `BESSELY(x, n)`.
pub fn bessely_kernel(x: f64, n: f64) -> Result<f64, WorksheetErrorCode> {
    if !x.is_finite() { return Err(WorksheetErrorCode::Num); }
    let n = order(n)?;
    let x = D::new(x);
    if !(x.0 > 0.0) {
        return Err(WorksheetErrorCode::Num);
    }
    finish(match n {
        0 => y0(x),
        1 => y1(x),
        _ => Ok(forward(x, n, y0(x)?, y1(x)?, false)),
    })
}

pub(super) fn y0(x: f64) -> f64 {
    let square = x * x;
    let rational = polynomial(square, Y0_NUM) / polynomial(square, Y0_DEN);
    rational + j0(x) * (W * excel_log(x))
}

pub(super) fn y1(x: f64) -> f64 {
    let square = x * x;
    let rational = (x * polynomial(square, Y1_NUM)) / polynomial(square, Y1_DEN);
    rational + W * (j1(x) * excel_log(x) - 1.0 / x)
}

pub(super) fn y1_large(x: f64) -> f64 {
    let square = x * x;
    let phase = x - 2.356194491;
    let reduction = 64.0 / square;
    let amplitude = polynomial(reduction, AMP1);
    let correction = polynomial(reduction, PHASE1);
    let leading = excel_sin(phase) * amplitude;
    let trailing = ((8.0 / x) * excel_cos(phase)) * correction;
    let combined = leading + trailing;
    let scale = (W / x).sqrt();
    scale * combined
}

pub(super) fn y0_large(x: f64) -> f64 {
    let square = x * x;
    let phase = x - 0.785398164;
    let reduction = 64.0 / square;
    let amplitude = polynomial(reduction, AMP0);
    let correction = polynomial(reduction, PHASE_Y0);
    let leading = excel_sin(phase) * amplitude;
    let trailing = ((8.0 / x) * excel_cos(phase)) * correction;
    let scale = (W / x).sqrt();
    scale * (leading + trailing)
}