predict_rs/
sun.rs

1use core::f64::consts::PI;
2
3use crate::consts::{ASTRONOMICAL_UNIT_KM, JULIAN_TIME_DIFF, SECONDS_PER_DAY};
4use crate::geodetic::Geodetic;
5use crate::julian_date::predict_to_julian_double;
6use crate::math::Vec3;
7use crate::orbit::calculate_obs;
8use crate::predict::{PredictObservation, PredictObserver};
9
10/// The function `delta_et` has been added to allow calculations on the position of the sun.  It provides the difference between UT (approximately the same as UTC) and ET (now referred to as TDT). This function is based on a least squares fit of data from 1950 to 1991 and will need to be updated periodically. Values determined using data from 1950-1991 in the 1990 Astronomical Almanac.  See `DELTA_ET.WQ1` for details.
11#[must_use]
12pub fn delta_et(year: f64) -> f64 {
13    26.465 + 0.747_622 * (year - 1950.0) + 1.886_913 * (2.0 * PI * (year - 1975.0) / 33.0).sin()
14}
15
16/// Returns angle in radians from argument in degrees.
17#[must_use]
18pub fn radians(degrees: f64) -> f64 {
19    degrees * PI / 180.0
20}
21
22/// Returns angle in degrees from argument in radians.
23#[must_use]
24pub fn degrees(radians: f64) -> f64 {
25    radians * 180.0 / PI
26}
27
28#[must_use]
29#[allow(clippy::many_single_char_names)] // Try not to mess with original codes
30pub fn sun_predict(time: f64) -> Vec3 {
31    let jul_utc = time + JULIAN_TIME_DIFF;
32    let mjd = jul_utc - 2_415_020.0;
33    let year = 1900.0 + mjd / 365.25;
34    let t = (mjd + delta_et(year) / SECONDS_PER_DAY) / 36525.0;
35    let m = radians(
36        358.475_83 + (35_999.049_75 * t % 360.0) - (0.000_150 + 0.000_003_3 * t) * t.sqrt() % 360.0,
37    );
38    let l = radians(279.696_68 + (36_000.768_92 * t % 360.0) + 0.000_302_5 * t.sqrt() % 360.0);
39    let e = 0.016_751_04 - (0.000_041_8 + 0.000_000_126 * t) * t;
40    let c = radians(
41        (1.919_460 - (0.004_789 + 0.000_014 * t) * t) * m.sin()
42            + (0.020_094 - 0.000_100 * t) * (2.0 * m).sin()
43            + 0.000_293 * (3.0 * m).sin(),
44    );
45    let o = radians(259.18 - 1934.142 * t % 360.0);
46    let lsa = (l + c - radians(0.00569 - 0.00479 * o.sin())) % (2.0 * PI);
47    let nu = (m + c) % (2.0 * PI);
48    let mut r = 1.000_000_2 * (1.0 - e.sqrt()) / (1.0 + e * nu.cos());
49    let eps = radians(
50        23.452_294 - (0.013_012_5 + (0.000_001_64 - 0.000_000_503 * t) * t) * t + 0.00256 * o.cos(),
51    );
52    r *= ASTRONOMICAL_UNIT_KM;
53
54    Vec3(
55        r * lsa.cos(),
56        r * lsa.sin() * eps.cos(),
57        r * lsa.sin() * eps.sin(),
58    )
59}
60
61#[must_use]
62pub fn predict_observe_sun(observer: &PredictObserver, time: f64) -> PredictObservation {
63    // Find sun position
64    let jultime = predict_to_julian_double(time) + JULIAN_TIME_DIFF;
65
66    let solar_vector = sun_predict(jultime);
67
68    let geodetic = Geodetic {
69        lat: observer.latitude,
70        lon: observer.longitude,
71        alt: observer.altitude / 1000.0,
72        theta: 0.0,
73    };
74
75    /* Zero vector for initializations */
76    let zero_vector = Vec3(0.0, 0.0, 0.0);
77
78    /* Solar observed azimuth and elevation vector  */
79    let solar_set = calculate_obs(jultime, solar_vector, zero_vector, geodetic);
80
81    let sun_azi = solar_set.0;
82    let sun_ele = solar_set.1;
83
84    let sun_range = 1.0 + ((solar_set.2 - ASTRONOMICAL_UNIT_KM) / ASTRONOMICAL_UNIT_KM);
85    let sun_range_rate = 1000.0 * solar_set.3;
86
87    PredictObservation {
88        time,
89        azimuth: sun_azi,
90        elevation: sun_ele,
91        range: sun_range,
92        range_rate: sun_range_rate,
93        azimuth_rate: 0.0,
94        elevation_rate: 0.0,
95        range_x: 0.0,
96        range_y: 0.0,
97        range_z: 0.0,
98        revolutions: 0.0,
99        visible: false,
100    }
101}
predict_rs/sun.rs

predict_rs/
sun.rs
