Skip to content

Commit

Permalink
Replace OrbitalState with ANISE::Orbit
Browse files Browse the repository at this point in the history
Signed-off-by: Guillaume W. Bres <[email protected]>
  • Loading branch information
gwbres committed Aug 13, 2024
1 parent d618a48 commit c707d21
Show file tree
Hide file tree
Showing 7 changed files with 36 additions and 114 deletions.
6 changes: 5 additions & 1 deletion src/bancroft.rs
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,11 @@ impl Bancroft {

for i in 0..4 {
let state = cd[i].state.ok_or(Error::UnresolvedState)?;
let (x_i, y_i, z_i) = (state.position[0], state.position[1], state.position[2]);
let (x_i, y_i, z_i) = (
state.radius_km.x * 1.0E3,
state.radius_km.y * 1.0E3,
state.radius_km.z * 1.0E3,
);
let r_i = cd[i]
.prefered_pseudorange()
.ok_or(Error::MissingPseudoRange)?
Expand Down
12 changes: 7 additions & 5 deletions src/candidate.rs
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,8 @@ use std::cmp::Ordering;
use crate::{
bias::RuntimeParams as BiasRuntimeParams,
prelude::{
Carrier, Config, Duration, Epoch, Error, IonoComponents, Method, OrbitalState,
TropoComponents, TropoModel, Vector3, SV,
Carrier, Config, Duration, Epoch, Error, IonoComponents, Method, Orbit, TropoComponents,
TropoModel, Vector3, SV,
},
};

Expand Down Expand Up @@ -94,7 +94,7 @@ pub struct Candidate {
/// Sampling [Epoch]
pub t: Epoch,
/// State that needs to be resolved
pub state: Option<OrbitalState>,
pub state: Option<Orbit>,
/// t_tx Epoch
pub(crate) t_tx: Epoch,
// SV group delay expressed as a [Duration]
Expand Down Expand Up @@ -199,8 +199,10 @@ impl Candidate {
method: Method,
tropo_modeling: bool,
iono_modeling: bool,
apriori_geo_ddeg: (f64, f64, f64),
apriori: Orbit,
almanac: &Almanac,
) {
let el_az_range = almanac.azimuth_elevation_range_sez(rx, tx)?;
if let Some(state) = self.state {
if let Some(obs) = self.prefered_pseudorange() {
let rtm = BiasRuntimeParams {
Expand Down Expand Up @@ -560,7 +562,7 @@ impl Candidate {
Ok((e_tx, dt))
}
#[cfg(test)]
pub fn set_state(&mut self, state: OrbitalState) {
pub fn set_state(&mut self, state: Orbit) {
self.state = Some(state);
}
}
Expand Down
4 changes: 2 additions & 2 deletions src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -35,15 +35,15 @@ pub mod prelude {
pub use crate::carrier::Carrier;
pub use crate::cfg::{Config, Method};
pub use crate::navigation::{Filter, InvalidationCause, PVTSolution, PVTSolutionType};
pub use crate::orbit::{OrbitalState, OrbitalStateProvider};
pub use crate::orbit::OrbitalStateProvider;
pub use crate::position::Position;
pub use crate::rtk::BaseStation;
pub use crate::solver::{Error, Solver};
// re-export
pub use anise::{
constants::frames::{EARTH_J2000, SUN_J2000},
naif::SPK,
prelude::{Aberration, Almanac, Frame},
prelude::{Aberration, Almanac, Frame, Orbit},
};
pub use gnss::prelude::{Constellation, SV};
pub use hifitime::{Duration, Epoch, TimeScale};
Expand Down
13 changes: 9 additions & 4 deletions src/navigation/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -133,14 +133,19 @@ impl Input {
sv_input.azimuth = azimuth;
sv_input.elevation = elevation;

let (sv_x, sv_y, sv_z) = (state.position[0], state.position[1], state.position[2]);
let (sv_x, sv_y, sv_z) = (
state.radius_km.x * 1.0E3,
state.radius_km.y * 1.0E3,
state.radius_km.z * 1.0E3,
);

let mut rho = ((sv_x - x0).powi(2) + (sv_y - y0).powi(2) + (sv_z - z0).powi(2)).sqrt();

if cfg.modeling.relativistic_path_range {
let mu = Constants::EARTH_GRAVITATION;
let r_sat = (state.position[0].powi(2)
+ state.position[1].powi(2)
+ state.position[2].powi(2))
let r_sat = ((state.radius_km.x * 1.0E3).powi(2)
+ (state.radius_km.y * 1.0E3).powi(2)
+ (state.radius_km.z * 1.0E3).powi(2))
.sqrt();
let r_0 = (x0.powi(2) + y0.powi(2) + z0.powi(2)).sqrt();
let r_sat_0 = r_0 - r_sat;
Expand Down
8 changes: 6 additions & 2 deletions src/navigation/solutions/validator.rs
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,12 @@ impl Validator {
x[3] / SPEED_OF_LIGHT_M_S,
);

let sv_pos = cd.state.unwrap().position;
let (sv_x, sv_y, sv_z) = (sv_pos[0], sv_pos[1], sv_pos[2]);
let sv_orbit = cd.state.unwrap();
let (sv_x, sv_y, sv_z) = (
sv_orbit.radius_km.x * 1.0E3,
sv_orbit.radius_km.y * 1.0E3,
sv_orbit.radius_km.z * 1.0E3,
);

let rho = ((sv_x - x).powi(2) + (sv_y - y).powi(2) + (sv_z - z).powi(2)).sqrt();

Expand Down
99 changes: 3 additions & 96 deletions src/orbit.rs
Original file line number Diff line number Diff line change
@@ -1,105 +1,12 @@
use anise::prelude::{Frame, Orbit};

use crate::prelude::{Epoch, Vector3, SV};
use map_3d::ecef2geodetic;
use std::f64::consts::PI;

/// [OrbitalState] must be provide of each candidate.
#[derive(Copy, Clone, Debug, Default, PartialEq)]
pub struct OrbitalState {
/// Elevation compared to reference position and horizon in [°]
pub elevation: f64,
/// Azimuth compared to reference position and magnetic North in [°]
pub azimuth: f64,
/// APC Position vector in [m] ECEF
pub position: Vector3<f64>,
// Velocity vector in [m/s] ECEF that we calculated ourselves
pub(crate) velocity: Option<Vector3<f64>>,
}
use crate::prelude::{Epoch, Orbit, SV};

/// OrbitalStateProvider must be implemented
/// and provide SV state at specified `t` for the solving process can proceed.
pub trait OrbitalStateProvider {
/// Provide [OrbitalState] at requested `t` so we can proceed.
/// Provide Antenna Phase Center state as [Orbit] at requested [Epoch] for requested [SV].
/// In case interpolation is used, we propose an interpolation order,
/// that would fit current setup, which you can choose to ignore.
/// If you can't provide (missing data?): simply return None.
/// If None is returned for too long, this [Epoch] will eventually get dropped out
/// and we will proceed to the next.
fn next_at(&mut self, t: Epoch, sv: SV, order: usize) -> Option<OrbitalState>;
}

impl OrbitalState {
/// Creates [Self] from Antenna Phase Center (APC) position coordinates,
/// expressed in ECEF [m].
pub fn from_position(position: (f64, f64, f64)) -> Self {
Self {
velocity: None,
azimuth: 0.0_f64,
elevation: 0.0_f64,
position: Vector3::<f64>::new(position.0, position.1, position.2),
}
}
/// Complete [Self] definition with both Elevation and Azimuth angles
pub(crate) fn with_elevation_azimuth(&self, position: (f64, f64, f64)) -> Self {
let (x, y, z) = position;
let (lat_rad, lon_rad, _) = ecef2geodetic(x, y, z, map_3d::Ellipsoid::WGS84);

let (sv_x, sv_y, sv_z) = (self.position[0], self.position[1], self.position[2]);
let a_i = (sv_x - x, sv_y - y, sv_z - z);
let norm = (a_i.0.powf(2.0) + a_i.1.powf(2.0) + a_i.2.powf(2.0)).sqrt();
let a_i = (a_i.0 / norm, a_i.1 / norm, a_i.2 / norm);

let ecef_to_ven = (
(
lat_rad.cos() * lon_rad.cos(),
lat_rad.cos() * lon_rad.sin(),
lat_rad.sin(),
),
(-lon_rad.sin(), lon_rad.cos(), 0.0_f64),
(
-lat_rad.sin() * lon_rad.cos(),
-lat_rad.sin() * lon_rad.sin(),
lat_rad.cos(),
),
);
// ECEF to VEN transform
let ven = (
(ecef_to_ven.0 .0 * a_i.0 + ecef_to_ven.0 .1 * a_i.1 + ecef_to_ven.0 .2 * a_i.2),
(ecef_to_ven.1 .0 * a_i.0 + ecef_to_ven.1 .1 * a_i.1 + ecef_to_ven.1 .2 * a_i.2),
(ecef_to_ven.2 .0 * a_i.0 + ecef_to_ven.2 .1 * a_i.1 + ecef_to_ven.2 .2 * a_i.2),
);
let elevation = (PI / 2.0 - ven.0.acos()).to_degrees();
let mut azimuth = (ven.1.atan2(ven.2)).to_degrees();
if azimuth < 0.0 {
azimuth += 360.0;
}
Self {
azimuth,
elevation,
position: self.position,
velocity: self.velocity,
}
}
pub(crate) fn velocity(&self) -> Option<Vector3<f64>> {
self.velocity
}
pub(crate) fn orbit(&self, dt: Epoch, frame: Frame) -> Orbit {
let p = self.position;
let v = self.velocity().unwrap_or_default();
Orbit::cartesian(
p[0] / 1.0E3,
p[1] / 1.0E3,
p[2] / 1.0E3,
v[0] / 1.0E3,
v[1] / 1.0E3,
v[2] / 1.0E3,
dt,
frame,
)
}
#[cfg(test)]
fn set_elevation(&mut self, elev: f64) {
self.elevation = elev;
}
fn next_at(&mut self, t: Epoch, sv: SV, order: usize) -> Option<Orbit>;
}
8 changes: 4 additions & 4 deletions src/solver.rs
Original file line number Diff line number Diff line change
Expand Up @@ -29,9 +29,9 @@ use crate::{
solutions::validator::{InvalidationCause, Validator as SolutionValidator},
Input as NavigationInput, Navigation, PVTSolution, PVTSolutionType,
},
orbit::{OrbitalState, OrbitalStateProvider},
orbit::OrbitalStateProvider,
position::Position,
prelude::{Duration, Epoch, Observation, SV},
prelude::{Duration, Epoch, Observation, Orbit, SV},
rtk::BaseStation,
};

Expand Down Expand Up @@ -606,7 +606,7 @@ impl<O: OrbitalStateProvider, B: BaseStation> Solver<O, B> {
}
}
/* rotate interpolated position */
fn rotate_position(rotate: bool, interpolated: OrbitalState, dt_tx: Duration) -> OrbitalState {
fn rotate_position(rotate: bool, interpolated: Orbit, dt_tx: Duration) -> Orbit {
let mut reworked = interpolated;
let rot = if rotate {
const EARTH_OMEGA_E_WGS84: f64 = 7.2921151467E-5;
Expand All @@ -627,7 +627,7 @@ impl<O: OrbitalStateProvider, B: BaseStation> Solver<O, B> {
/*
* Determine velocities
*/
fn velocities(&self, t_tx: Epoch, sv: SV, interpolated: OrbitalState) -> OrbitalState {
fn orbit_with_velocities(&self, t_tx: Epoch, sv: SV, interpolated: Orbit) -> Orbit {
let mut reworked = interpolated;
if let Some((p_ttx, p_pos)) = self.prev_sv_state.get(&sv) {
let dt = (t_tx - *p_ttx).to_seconds();
Expand Down

0 comments on commit c707d21

Please sign in to comment.