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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
/*
* ANISE Toolkit
* Copyright (C) 2021-onward Christopher Rabotin <christopher.rabotin@gmail.com> et al. (cf. AUTHORS.md)
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at https://mozilla.org/MPL/2.0/.
*
* Documentation: https://nyxspace.com/
*/
use crate::{
astro::Aberration, constants::frames::SUN_J2000, ephemerides::EphemerisError, prelude::Frame,
NaifId,
};
use super::Almanac;
use hifitime::Epoch;
#[cfg(feature = "python")]
use pyo3::prelude::*;
#[cfg_attr(feature = "python", pymethods)]
impl Almanac {
/// Returns the angular separation (between 0 and 180 degrees) between the observer and the Sun, and the observer and the target body ID.
/// This is formally known as the "solar elongation".
/// This computes the Sun Probe Earth angle (SPE) if the probe is in a loaded SPK, its ID is the "observer_id", and the target is set to its central body.
///
/// # Geometry
/// If the SPE is greater than 90 degrees, then the celestial object below the probe is in sunlight.
///
/// This angle determines the illumination phase of the target as seen by the observer:
/// * **~0° (Conjunction):** The Target is in the same direction as the Sun. The observer sees the unlit side ("New Moon").
/// * **~180° (Opposition):** The Target is in the opposite direction of the Sun. The observer sees the fully lit side ("Full Moon").
/// * **> 90°:** The observer is generally on the "day" side of the target.
///
///
/// ## Sunrise at nadir
/// ```text
/// Sun
/// | \
/// | \
/// | \
/// Obs. -- Target
/// ```
/// ## Sun high at nadir
/// ```text
/// Sun
/// \
/// \ __ θ > 90
/// \ \
/// Obs. ---------- Target
/// ```
///
/// ## Sunset at nadir
/// ```text
/// Sun
/// /
/// / __ θ < 90
/// / /
/// Obs. -- Target
/// ```
///
/// # Algorithm
/// 1. Compute the position of the Sun as seen from the observer
/// 2. Compute the position of the target as seen from the observer
/// 3. Return the arccosine of the dot product of the norms of these vectors.
///
/// :type target_id: int
/// :type observer_id: int
/// :type epoch: Epoch
/// :type ab_corr: Aberration
/// :rtype: float
pub fn sun_angle_deg(
&self,
target_id: NaifId,
observer_id: NaifId,
epoch: Epoch,
ab_corr: Option<Aberration>,
) -> Result<f64, EphemerisError> {
let obs_to_sun = self.translate(
SUN_J2000,
Frame::from_ephem_j2000(observer_id),
epoch,
ab_corr,
)?;
let obs_to_target = self.translate(
Frame::from_ephem_j2000(target_id),
Frame::from_ephem_j2000(observer_id),
epoch,
ab_corr,
)?;
Ok(obs_to_sun
.r_hat()
.dot(&obs_to_target.r_hat())
.acos()
.to_degrees())
}
/// Convenience function that calls `sun_angle_deg` with the provided frames instead of the ephemeris ID.
///
/// :type target: Frame
/// :type observer: Frame
/// :type epoch: Epoch
/// :type ab_corr: Aberration
/// :rtype: float
pub fn sun_angle_deg_from_frame(
&self,
target: Frame,
observer: Frame,
epoch: Epoch,
ab_corr: Option<Aberration>,
) -> Result<f64, EphemerisError> {
self.sun_angle_deg(target.ephemeris_id, observer.ephemeris_id, epoch, ab_corr)
}
}
#[cfg(test)]
mod ut_solar {
use crate::{
constants::{
celestial_objects::EARTH,
frames::{EARTH_J2000, IAU_EARTH_FRAME, SUN_J2000},
},
prelude::*,
};
/// Load a BSP of a spacecraft trajectory, compute the sun elevation at different points on the surface below that spacecraft
/// and ensure that it matches with the SPE calculation.
#[test]
fn verify_geometry() {
let ctx = Almanac::default()
.load("../data/de440s.bsp")
.and_then(|ctx| ctx.load("../data/gmat-hermite.bsp"))
.and_then(|ctx| ctx.load("../data/pck11.pca"))
.unwrap();
let epoch = Epoch::from_gregorian_hms(2000, 1, 1, 12, 0, 0, TimeScale::UTC);
let sc_id = -10000001;
let my_sc_j2k = Frame::from_ephem_j2000(sc_id);
// Grab the state in the J2000 frame
let state = ctx.transform(my_sc_j2k, EARTH_J2000, epoch, None).unwrap();
// We'll check at four different points in the orbit
for epoch in TimeSeries::inclusive(
epoch,
epoch + state.period().unwrap(),
0.05 * state.period().unwrap(),
) {
let spe_deg = ctx
.sun_angle_deg(EARTH, sc_id, epoch, Aberration::NONE)
.unwrap();
assert_eq!(
spe_deg,
ctx.sun_angle_deg_from_frame(EARTH_J2000, my_sc_j2k, epoch, Aberration::NONE)
.unwrap()
);
let iau_earth = ctx.frame_info(IAU_EARTH_FRAME).unwrap();
// Compute this state in the body fixed frame.
let state_bf = ctx
.transform(my_sc_j2k, IAU_EARTH_FRAME, epoch, None)
.unwrap();
// Build a local point on the surface below this spacecraft
let nadir_surface_point = Orbit::try_latlongalt(
state_bf.latitude_deg().unwrap(),
state_bf.longitude_deg(),
0.0,
epoch,
iau_earth,
)
.unwrap();
// Fetch the state of the sun at this time.
let sun_state = ctx
.transform(SUN_J2000, IAU_EARTH_FRAME, epoch, None)
.unwrap();
// Compute the Sun elevation from that point.
let sun_elevation_deg = ctx
.azimuth_elevation_range_sez(sun_state, nadir_surface_point, None, None)
.unwrap()
.elevation_deg;
println!(
"{epoch}\tsun el = {sun_elevation_deg:.3} deg\tlat = {:.3} deg\tlong = {:.3} deg\t-->\tSPE = {spe_deg:.3} deg",
nadir_surface_point.latitude_deg().unwrap(),
nadir_surface_point.longitude_deg()
);
// Test: the sun elevation from the ground should be about 90 degrees _less_ than the SPE
// The "about" is because the SPE effectively assumes a spherical celestial object, but the frame
// from the Almanac accounts for the flattening when computing the exact location below the vehicle.
// Hence, there is a small difference (we accept up to 0.05 degrees of difference).
assert!((sun_elevation_deg + 90.0 - spe_deg).abs() < 5e-2)
}
}
}