rust-roche 0.4.4

use crate::errors::RocheError;
use crate::{Star, Vec3};
use crate::{dbrent, ref_sphere, rpot_grad, rpot_val, sphere_eclipse_vector};
use pyo3::prelude::*;

///
/// fblink works out whether or not a given point is eclipsed by a Roche-distorted star by searching along the line
/// of sight to the point to see if the Roche potential ever drops below the value at the stellar surface.
///
/// Arguments:
///
/// * `q`:      mass ratio = M2/M1
/// * `star`:   star concerned
/// * `spin`:   ratio of spin to orbital frequency
/// * `star`:   which star is doing the eclipsing, primary or secondary
/// * `ffac`:   the filling factor of the star
/// * `acc`:    accuracy of location of minimum potential, units of separation. The accuracy in height relative to the
/// Roche potential is acc*acc/(2*R) where R is the radius of curvature of the Roche potential surface, so don't be too
/// picky. 1.e-4 would be more than good enough in most cases.
/// * `earth`:  vector pointing towards earth
/// * `p`:      point of interest
///
/// Returns:
///
/// * true if minimum potential is below the potential at stellar surface
///
#[pyfunction]
pub fn fblink(
    q: f64,
    star: Star,
    spin: f64,
    ffac: f64,
    acc: f64,
    earth: &Vec3,
    p: &Vec3,
) -> Result<bool, RocheError> {
    let (rref, pref) = ref_sphere(q, star, spin, ffac)?;

    let cofm: Vec3 = match star {
        Star::Primary => Vec3::cofm1(),
        Star::Secondary => Vec3::cofm2(),
    };

    // First compute the multipliers cutting the reference sphere (if any)
    let mut lam1 = 0.0;
    let mut lam2 = 0.0;
    if !sphere_eclipse_vector(earth, p, &cofm, rref, &mut lam1, &mut lam2) {
        return Ok(false);
    }
    if lam1 == 0.0 {
        return Ok(true);
    }

    // Create function objects for 1D minimisation in lambda direction
    let func = |lam: f64| rpot_val(q, star, spin, earth, p, lam);

    // Now try to bracket a minimum. We just crudely compute function at regularly spaced intervals filling in the
    // gaps until the step size between the points drops below the threshold. Take every opportunity to jump out early
    // either if the potential is below the threshold or if we have bracketed a minimum.
    let mut nstep: i32 = 1;
    let mut step: f64 = lam2 - lam1;

    let mut f1: f64 = 0.0;
    let mut f2: f64 = 0.0;
    let mut flam: f64 = 1.0;
    let mut lam: f64 = lam1;

    while step > acc {
        lam = lam1 + step / 2.0;

        for _ in 0..nstep {
            flam = func(lam)?;
            if flam <= pref {
                return Ok(true);
            }

            // Calculate these as late as possible because they may often not be needed
            if nstep == 1 {
                f1 = func(lam1)?;
                f2 = func(lam2)?;
            }

            if flam < f1 && flam < f2 {
                break;
            }

            lam += step;
        }
        if flam < f1 && flam < f2 {
            break;
        }
        step /= 2.0;
        nstep *= 2;
    }

    if flam < f1 && flam < f2 {
        // OK, minimum bracketted, so finally pin it down accurately
        // Possible that multiple minima could cause problems but I have
        // never seen this in practice.
        let dfunc = |lam: f64| {
            let (_dp, dl) = rpot_grad(q, star, spin, earth, p, lam)?;
            Ok(dl)
        };

        // this may fail
        let (_xmin, flam) = dbrent(lam1, lam, lam2, func, dfunc, acc, true, pref)?;

        Ok(flam < pref)
    } else {
        // Not bracketted even after a detailed search, and we have not jumped
        // out either, so assume no eclipse
        Ok(false)
    }
}