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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

//! This crate is a very incomplete implementation of the book
//! Astronomical Algorithms by Jean Meeus.

/// Utilities for working with angles
pub mod angles;

/// Functions related to dates
pub mod date;

/// General interpolation functions
pub mod interpolation;

/// Celestial coordinate systems
mod celestial_coordinates;

/// Astronomical constants
mod constants;

/// Regression (curve-fitting) algorithms
pub mod regression;

/// Iteration algorithms
pub mod iteration;

/// This computes the apparent magnitude of a comet
///
/// The distances to the sun and to earth are in astronomical units.
/// `absolute_magnitude` and `coefficient` need to be deduced from
/// observations.
pub fn comet_apparent_magnitude(
    absolute_magnitude: f64,
    coefficient: f64,
    dist_to_sun: f64,
    dist_to_earth: f64,
) -> f64 {
    absolute_magnitude + 5.0 * dist_to_earth.log10() + coefficient * dist_to_sun.log10()
}

/// Given a number of comet observations, this function performs
/// linear regression to find the absolute magnitude and coefficient of
/// a comet. Distances are in astronomical units.
pub fn comet_magnitude_regression(
    apparent_magnitudes: &[f64],
    dists_to_sun: &[f64],
    dists_to_earth: &[f64],
) -> (f64, f64) {
    assert!(!apparent_magnitudes.is_empty());
    assert_eq!(apparent_magnitudes.len(), dists_to_sun.len());
    assert_eq!(dists_to_sun.len(), dists_to_earth.len());
    let xs = dists_to_sun.iter().map(|d| d.log10());
    let ys = apparent_magnitudes
        .iter()
        .zip(dists_to_earth)
        .map(|(m, d)| m - 5.0 * d.log10());
    let (slope, intercept) = regression::linear(xs.zip(ys));
    let absolute_magnitude = intercept;
    let coefficient = slope;
    (absolute_magnitude, coefficient)
}

#[cfg(test)]
#[track_caller]
fn assert_approx_equal(a: f64, b: f64) {
    if (a - b).abs() > 1e-5 {
        panic!("{} != {}", a, b);
    }
}

#[cfg(test)]
mod tests {}