1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

use crate::binarydata as pa_bd;
use crate::macros as pa_m;
use crate::util as pa_u;

/// Calculate orbital data for binary star.
///
/// ## Arguments
/// * `greenwich_date_day` -- Greenwich date (day)
/// * `greenwich_date_month` -- Greenwich date (month)
/// * `greenwich_date_year` -- Greenwich date (year)
/// * `binary_name` -- Abbreviated name of binary
///
/// ## Returns
/// * `position_angle_deg` -- Position angle (degrees)
/// * `separation_arcsec` -- Separation of binary members (arcseconds)
pub fn binary_star_orbit(
    greenwich_date_day: f64,
    greenwich_date_month: u32,
    greenwich_date_year: u32,
    binary_name: String,
) -> (f64, f64) {
    let (binary_info, _binary_info_status) = pa_bd::get_binary_info_vector(binary_name);

    let y_years = (greenwich_date_year as f64
        + (pa_m::cd_jd(
            greenwich_date_day,
            greenwich_date_month,
            greenwich_date_year,
        ) - pa_m::cd_jd(0.0, 1, greenwich_date_year))
            / 365.242191)
        - binary_info.epoch_peri;
    let m_deg = 360.0 * y_years / binary_info.period;
    let m_rad = (m_deg - 360.0 * (m_deg / 360.0).floor()).to_radians();
    let eccentricity = binary_info.ecc;
    let true_anomaly_rad = pa_m::true_anomaly(m_rad, eccentricity);
    let r_arcsec = (1.0 - eccentricity * (pa_m::eccentric_anomaly(m_rad, eccentricity)).cos())
        * binary_info.axis;
    let ta_peri_rad = true_anomaly_rad + binary_info.long_peri.to_radians();

    let y = (ta_peri_rad).sin() * ((binary_info.incl).to_radians()).cos();
    let x = (ta_peri_rad).cos();
    let a_deg = pa_m::degrees(y.atan2(x));
    let theta_deg1 = a_deg + binary_info.pa_node;
    let theta_deg2 = theta_deg1 - 360.0 * (theta_deg1 / 360.0).floor();
    let rho_arcsec =
        r_arcsec * (ta_peri_rad).cos() / ((theta_deg2 - binary_info.pa_node).to_radians()).cos();

    let position_angle_deg = pa_u::round_f64(theta_deg2, 1);
    let separation_arcsec = pa_u::round_f64(rho_arcsec, 2);

    return (position_angle_deg, separation_arcsec);
}