Type Alias Orbit

pub type Orbit = CartesianState;
Expand description

A helper type alias, but no assumptions are made on the underlying validity of the frame.

Aliased Type§

struct Orbit {
    pub radius_km: Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>,
    pub velocity_km_s: Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>,
    pub epoch: Epoch,
    pub frame: Frame,
}

Fields§

§radius_km: Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>

Position radius in kilometers

§velocity_km_s: Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>

Velocity in kilometers per second

§epoch: Epoch

Epoch with time scale at which this is valid.

§frame: Frame

Frame in which this Cartesian state lives.

Implementations

§

impl CartesianState

pub fn dcm_from_topocentric_to_body_fixed( &self, _from: i32, ) -> Result<DCM, PhysicsError>

Builds the rotation matrix that rotates from the topocentric frame (SEZ) into the body fixed frame of this state.

§Frame warnings
  • If the state is NOT in a body fixed frame (i.e. ITRF93), then this computation is INVALID.
  • (Usually) no time derivative can be computed: the orbit is expected to be a body fixed frame where the at_epoch function will fail. Exceptions for Moon body fixed frames.
§UNUSED Arguments
  • from: ID of this new frame. Only used to set the “from” frame of the DCM. – No longer used since 0.5.3
§Source

From the GMAT MathSpec, page 30 section 2.6.9 and from Calculate_RFT in TopocentricAxes.cpp, this returns the rotation matrix from the topocentric frame (SEZ) to body fixed frame. In the GMAT MathSpec notation, R_{IF} is the DCM from body fixed to inertial. Similarly, R{FT} is from topocentric to body fixed.

pub fn dcm3x3_from_topocentric_to_body_fixed(&self) -> Result<DCM, PhysicsError>

Builds the rotation matrix that rotates from the topocentric frame (SEZ) into the body fixed frame of this state.

§Frame warning

If the state is NOT in a body fixed frame (i.e. ITRF93), then this computation is INVALID.

§Source

From the GMAT MathSpec, page 30 section 2.6.9 and from Calculate_RFT in TopocentricAxes.cpp, this returns the rotation matrix from the topocentric frame (SEZ) to body fixed frame. In the GMAT MathSpec notation, R_{IF} is the DCM from body fixed to inertial. Similarly, R{FT} is from topocentric to body fixed.

pub fn dcm_from_ric_to_inertial(&self) -> Result<DCM, PhysicsError>

Builds the rotation matrix that rotates from this state’s inertial frame to this state’s RIC frame

§Frame warning

If the state is NOT in an inertial frame, then this computation is INVALID.

§Algorithm
  1. Compute the state data one millisecond before and one millisecond assuming two body dynamics
  2. Compute the DCM for this state, and the pre and post states
  3. Build the c vector as the normalized orbital momentum vector
  4. Build the i vector as the cross product of \hat{r} and c
  5. Build the RIC DCM as a 3x3 of the columns [\hat{r}, \hat{i}, \hat{c}], for the post, post, and current states
  6. Compute the difference between the DCMs of the pre and post states, to build the DCM time derivative
  7. Return the DCM structure with a 6x6 state DCM.
§Note on the time derivative

If the pre or post states cannot be computed, then the time derivative of the DCM will not be set. Further note that most astrodynamics tools do not account for the time derivative in the RIC frame.

pub fn dcm3x3_from_ric_to_inertial(&self) -> Result<DCM, PhysicsError>

Builds the rotation matrix that rotates from this state’s inertial frame to this state’s RIC frame

§Frame warning

If the state is NOT in an inertial frame, then this computation is INVALID.

§Algorithm
  1. Build the c vector as the normalized orbital momentum vector
  2. Build the i vector as the cross product of \hat{r} and c
  3. Build the RIC DCM as a 3x3 of the columns [\hat{r}, \hat{i}, \hat{c}]
  4. Return the DCM structure without accounting for the transport theorem.

pub fn dcm3x3_from_rcn_to_inertial(&self) -> Result<DCM, PhysicsError>

Builds the rotation matrix that rotates from this state’s inertial frame to this state’s RCN frame (radial, cross, normal)

§Frame warning

If the stattion is NOT in an inertial frame, then this computation is INVALID.

§Algorithm
  1. Compute \hat{r}, \hat{h}, the unit vectors of the radius and orbital momentum.
  2. Compute the cross product of these
  3. Build the DCM with these unit vectors
  4. Return the DCM structure

pub fn dcm_from_rcn_to_inertial(&self) -> Result<DCM, PhysicsError>

Builds the rotation matrix that rotates from this state’s inertial frame to this state’s RCN frame (radial, cross, normal)

§Frame warning

If the stattion is NOT in an inertial frame, then this computation is INVALID.

§Algorithm
  1. Compute \hat{r}, \hat{h}, the unit vectors of the radius and orbital momentum.
  2. Compute the cross product of these
  3. Build the DCM with these unit vectors
  4. Return the DCM structure with a 6x6 DCM with the time derivative of the VNC frame set.
§Note on the time derivative

If the pre or post states cannot be computed, then the time derivative of the DCM will not be set. Further note that most astrodynamics tools do not account for the time derivative in the RIC frame.

pub fn dcm3x3_from_vnc_to_inertial(&self) -> Result<DCM, PhysicsError>

Builds the rotation matrix that rotates from this state’s inertial frame to this state’s VNC frame (velocity, normal, cross)

§Frame warning

If the stattion is NOT in an inertial frame, then this computation is INVALID.

§Algorithm
  1. Compute \hat{v}, \hat{h}, the unit vectors of the radius and orbital momentum.
  2. Compute the cross product of these
  3. Build the DCM with these unit vectors
  4. Return the DCM structure.

pub fn dcm_from_vnc_to_inertial(&self) -> Result<DCM, PhysicsError>

Builds the rotation matrix that rotates from this state’s inertial frame to this state’s VNC frame (velocity, normal, cross)

§Frame warning

If the stattion is NOT in an inertial frame, then this computation is INVALID.

§Algorithm
  1. Compute \hat{v}, \hat{h}, the unit vectors of the radius and orbital momentum.
  2. Compute the cross product of these
  3. Build the DCM with these unit vectors
  4. Compute the difference between the DCMs of the pre and post states (+/- 1 ms), to build the DCM time derivative
  5. Return the DCM structure with a 6x6 DCM with the time derivative of the VNC frame set.
§Note on the time derivative

If the pre or post states cannot be computed, then the time derivative of the DCM will not be set. Further note that most astrodynamics tools do not account for the time derivative in the RIC frame.

pub fn hx(&self) -> Result<f64, PhysicsError>

Returns the orbital momentum value on the X axis

:rtype: float

pub fn hy(&self) -> Result<f64, PhysicsError>

Returns the orbital momentum value on the Y axis

:rtype: float

pub fn hz(&self) -> Result<f64, PhysicsError>

Returns the orbital momentum value on the Z axis

:rtype: float

pub fn hmag(&self) -> Result<f64, PhysicsError>

Returns the norm of the orbital momentum

:rtype: float

pub fn energy_km2_s2(&self) -> Result<f64, PhysicsError>

Returns the specific mechanical energy in km^2/s^2

:rtype: float

pub fn sma_km(&self) -> Result<f64, PhysicsError>

Returns the semi-major axis in km

:rtype: float

Examples found in repository?
examples/03_geo_analysis/drift.rs (line 118)
26fn main() -> Result<(), Box<dyn Error>> {
27    pel::init();
28    // Dynamics models require planetary constants and ephemerides to be defined.
29    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30    // This will automatically download the DE440s planetary ephemeris,
31    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
32    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
33    // planetary constants kernels.
34    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
35    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
36    // references to many functions.
37    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
38    // Define the orbit epoch
39    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
40
41    // Define the orbit.
42    // First we need to fetch the Earth J2000 from information from the Almanac.
43    // This allows the frame to include the gravitational parameters and the shape of the Earth,
44    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
45    // by loading a different set of planetary constants.
46    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
47
48    // Placing this GEO bird just above Colorado.
49    // In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
50    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
51    // Print in in Keplerian form.
52    println!("{orbit:x}");
53
54    let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
55    let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;
56
57    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
58    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
59    // models such as solar radiation pressure.
60
61    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
62    let sc = Spacecraft::builder()
63        .orbit(orbit)
64        .mass(Mass::from_dry_mass(9.60))
65        .srp(SRPData {
66            area_m2: 10e-4,
67            coeff_reflectivity: 1.1,
68        })
69        .build();
70    println!("{sc:x}");
71
72    // Set up the spacecraft dynamics.
73
74    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
75    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
76    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
77
78    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
79    // We're using the JGM3 model here, which is the default in GMAT.
80    let mut jgm3_meta = MetaFile {
81        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
82        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
83    };
84    // And let's download it if we don't have it yet.
85    jgm3_meta.process(true)?;
86
87    // Build the spherical harmonics.
88    // The harmonics must be computed in the body fixed frame.
89    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
90    let harmonics_21x21 = Harmonics::from_stor(
91        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
92        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
93    );
94
95    // Include the spherical harmonics into the orbital dynamics.
96    orbital_dyn.accel_models.push(harmonics_21x21);
97
98    // We define the solar radiation pressure, using the default solar flux and accounting only
99    // for the eclipsing caused by the Earth and Moon.
100    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
103    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
104    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
105
106    println!("{dynamics}");
107
108    // Finally, let's propagate this orbit to the same epoch as above.
109    // The first returned value is the spacecraft state at the final epoch.
110    // The second value is the full trajectory where the step size is variable step used by the propagator.
111    let (future_sc, trajectory) = Propagator::default(dynamics)
112        .with(sc, almanac.clone())
113        .until_epoch_with_traj(epoch + Unit::Century * 0.03)?;
114
115    println!("=== High fidelity propagation ===");
116    println!(
117        "SMA changed by {:.3} km",
118        orbit.sma_km()? - future_sc.orbit.sma_km()?
119    );
120    println!(
121        "ECC changed by {:.6}",
122        orbit.ecc()? - future_sc.orbit.ecc()?
123    );
124    println!(
125        "INC changed by {:.3e} deg",
126        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
127    );
128    println!(
129        "RAAN changed by {:.3} deg",
130        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
131    );
132    println!(
133        "AOP changed by {:.3} deg",
134        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
135    );
136    println!(
137        "TA changed by {:.3} deg",
138        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
139    );
140
141    // We also have access to the full trajectory throughout the propagation.
142    println!("{trajectory}");
143
144    println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");
145
146    // With the trajectory, let's build a few data products.
147
148    // 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.
149
150    let analysis_step = Unit::Minute * 5;
151
152    trajectory.to_parquet(
153        "./03_geo_hf_prop.parquet",
154        Some(vec![
155            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
156        ]),
157        ExportCfg::builder().step(analysis_step).build(),
158        almanac.clone(),
159    )?;
160
161    // 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.
162
163    // We iterate over the trajectory, grabbing a state every two minutes.
164    let mut offset_s = vec![];
165    let mut epoch_str = vec![];
166    let mut longitude_deg = vec![];
167    let mut latitude_deg = vec![];
168    let mut altitude_km = vec![];
169
170    for state in trajectory.every(analysis_step) {
171        // Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
172        // These define the GEO stationkeeping box.
173
174        let this_epoch = state.epoch();
175
176        offset_s.push((this_epoch - orbit.epoch).to_seconds());
177        epoch_str.push(this_epoch.to_isoformat());
178
179        let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
180        let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
181        longitude_deg.push(long_deg);
182        latitude_deg.push(lat_deg);
183        altitude_km.push(alt_km);
184    }
185
186    println!(
187        "Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
188        orig_long_deg - longitude_deg.last().unwrap()
189    );
190
191    println!(
192        "Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
193        orig_lat_deg - latitude_deg.last().unwrap()
194    );
195
196    println!(
197        "Altitude changed by {:.3} km -- Box is 30 km",
198        orig_alt_km - altitude_km.last().unwrap()
199    );
200
201    // Build the station keeping data frame.
202    let mut sk_df = df!(
203        "Offset (s)" => offset_s.clone(),
204        "Epoch (UTC)" => epoch_str.clone(),
205        "Longitude E-W (deg)" => longitude_deg,
206        "Latitude N-S (deg)" => latitude_deg,
207        "Altitude (km)" => altitude_km,
208
209    )?;
210
211    // Create a file to write the Parquet to
212    let file = File::create("./03_geo_lla.parquet").expect("Could not create file");
213
214    // Create a ParquetWriter and write the DataFrame to the file
215    ParquetWriter::new(file).finish(&mut sk_df)?;
216
217    Ok(())
218}
More examples
Hide additional examples
examples/01_orbit_prop/main.rs (line 67)
30fn main() -> Result<(), Box<dyn Error>> {
31    pel::init();
32    // Dynamics models require planetary constants and ephemerides to be defined.
33    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
34    // This will automatically download the DE440s planetary ephemeris,
35    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
36    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
37    // planetary constants kernels.
38    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
39    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
40    // references to many functions.
41    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
42    // Define the orbit epoch
43    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45    // Define the orbit.
46    // First we need to fetch the Earth J2000 from information from the Almanac.
47    // This allows the frame to include the gravitational parameters and the shape of the Earth,
48    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
49    // by loading a different set of planetary constants.
50    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
51
52    let orbit =
53        Orbit::try_keplerian_altitude(300.0, 0.015, 68.5, 65.2, 75.0, 0.0, epoch, earth_j2000)?;
54    // Print in in Keplerian form.
55    println!("{orbit:x}");
56
57    // There are two ways to propagate an orbit. We can make a quick approximation assuming only two-body
58    // motion. This is a useful first order approximation but it isn't used in real-world applications.
59
60    // This approach is a feature of ANISE.
61    let future_orbit_tb = orbit.at_epoch(epoch + Unit::Day * 3)?;
62    println!("{future_orbit_tb:x}");
63
64    // Two body propagation relies solely on Kepler's laws, so only the true anomaly will change.
65    println!(
66        "SMA changed by {:.3e} km",
67        orbit.sma_km()? - future_orbit_tb.sma_km()?
68    );
69    println!(
70        "ECC changed by {:.3e}",
71        orbit.ecc()? - future_orbit_tb.ecc()?
72    );
73    println!(
74        "INC changed by {:.3e} deg",
75        orbit.inc_deg()? - future_orbit_tb.inc_deg()?
76    );
77    println!(
78        "RAAN changed by {:.3e} deg",
79        orbit.raan_deg()? - future_orbit_tb.raan_deg()?
80    );
81    println!(
82        "AOP changed by {:.3e} deg",
83        orbit.aop_deg()? - future_orbit_tb.aop_deg()?
84    );
85    println!(
86        "TA changed by {:.3} deg",
87        orbit.ta_deg()? - future_orbit_tb.ta_deg()?
88    );
89
90    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
91    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
92    // models such as solar radiation pressure.
93
94    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
95    let sc = Spacecraft::builder()
96        .orbit(orbit)
97        .mass(Mass::from_dry_mass(9.60))
98        .srp(SRPData {
99            area_m2: 10e-4,
100            coeff_reflectivity: 1.1,
101        })
102        .build();
103    println!("{sc:x}");
104
105    // Set up the spacecraft dynamics.
106
107    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
108    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
109    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
110
111    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
112    // We're using the JGM3 model here, which is the default in GMAT.
113    let mut jgm3_meta = MetaFile {
114        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
115        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
116    };
117    // And let's download it if we don't have it yet.
118    jgm3_meta.process(true)?;
119
120    // Build the spherical harmonics.
121    // The harmonics must be computed in the body fixed frame.
122    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
123    let harmonics_21x21 = Harmonics::from_stor(
124        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
125        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
126    );
127
128    // Include the spherical harmonics into the orbital dynamics.
129    orbital_dyn.accel_models.push(harmonics_21x21);
130
131    // We define the solar radiation pressure, using the default solar flux and accounting only
132    // for the eclipsing caused by the Earth.
133    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
134
135    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
136    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
137    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
138
139    println!("{dynamics}");
140
141    // Finally, let's propagate this orbit to the same epoch as above.
142    // The first returned value is the spacecraft state at the final epoch.
143    // The second value is the full trajectory where the step size is variable step used by the propagator.
144    let (future_sc, trajectory) = Propagator::default(dynamics)
145        .with(sc, almanac.clone())
146        .until_epoch_with_traj(future_orbit_tb.epoch)?;
147
148    println!("=== High fidelity propagation ===");
149    println!(
150        "SMA changed by {:.3} km",
151        orbit.sma_km()? - future_sc.orbit.sma_km()?
152    );
153    println!(
154        "ECC changed by {:.6}",
155        orbit.ecc()? - future_sc.orbit.ecc()?
156    );
157    println!(
158        "INC changed by {:.3e} deg",
159        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
160    );
161    println!(
162        "RAAN changed by {:.3} deg",
163        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
164    );
165    println!(
166        "AOP changed by {:.3} deg",
167        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
168    );
169    println!(
170        "TA changed by {:.3} deg",
171        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
172    );
173
174    // We also have access to the full trajectory throughout the propagation.
175    println!("{trajectory}");
176
177    // With the trajectory, let's build a few data products.
178
179    // 1. Export the trajectory as a CCSDS OEM version 2.0 file and as a parquet file, which includes the Keplerian orbital elements.
180
181    trajectory.to_oem_file(
182        "./01_cubesat_hf_prop.oem",
183        ExportCfg::builder().step(Unit::Minute * 2).build(),
184    )?;
185
186    trajectory.to_parquet_with_cfg(
187        "./01_cubesat_hf_prop.parquet",
188        ExportCfg::builder().step(Unit::Minute * 2).build(),
189        almanac.clone(),
190    )?;
191
192    // 2. Compare the difference in the radial-intrack-crosstrack frame between the high fidelity
193    // and Keplerian propagation. The RIC frame is commonly used to compute the difference in position
194    // and velocity of different spacecraft.
195    // 3. Compute the azimuth, elevation, range, and range-rate data of that spacecraft as seen from Boulder, CO, USA.
196
197    let boulder_station = GroundStation::from_point(
198        "Boulder, CO, USA".to_string(),
199        40.014984,   // latitude in degrees
200        -105.270546, // longitude in degrees
201        1.6550,      // altitude in kilometers
202        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
203    );
204
205    // We iterate over the trajectory, grabbing a state every two minutes.
206    let mut offset_s = vec![];
207    let mut epoch_str = vec![];
208    let mut ric_x_km = vec![];
209    let mut ric_y_km = vec![];
210    let mut ric_z_km = vec![];
211    let mut ric_vx_km_s = vec![];
212    let mut ric_vy_km_s = vec![];
213    let mut ric_vz_km_s = vec![];
214
215    let mut azimuth_deg = vec![];
216    let mut elevation_deg = vec![];
217    let mut range_km = vec![];
218    let mut range_rate_km_s = vec![];
219    for state in trajectory.every(Unit::Minute * 2) {
220        // Try to compute the Keplerian/two body state just in time.
221        // This method occasionally fails to converge on an appropriate true anomaly
222        // from the mean anomaly. If that happens, we just skip this state.
223        // The high fidelity and Keplerian states diverge continuously, and we're curious
224        // about the divergence in this quick analysis.
225        let this_epoch = state.epoch();
226        match orbit.at_epoch(this_epoch) {
227            Ok(tb_then) => {
228                offset_s.push((this_epoch - orbit.epoch).to_seconds());
229                epoch_str.push(format!("{this_epoch}"));
230                // Compute the two body state just in time.
231                let ric = state.orbit.ric_difference(&tb_then)?;
232                ric_x_km.push(ric.radius_km.x);
233                ric_y_km.push(ric.radius_km.y);
234                ric_z_km.push(ric.radius_km.z);
235                ric_vx_km_s.push(ric.velocity_km_s.x);
236                ric_vy_km_s.push(ric.velocity_km_s.y);
237                ric_vz_km_s.push(ric.velocity_km_s.z);
238
239                // Compute the AER data for each state.
240                let aer = almanac.azimuth_elevation_range_sez(
241                    state.orbit,
242                    boulder_station.to_orbit(this_epoch, &almanac)?,
243                    None,
244                    None,
245                )?;
246                azimuth_deg.push(aer.azimuth_deg);
247                elevation_deg.push(aer.elevation_deg);
248                range_km.push(aer.range_km);
249                range_rate_km_s.push(aer.range_rate_km_s);
250            }
251            Err(e) => warn!("{} {e}", state.epoch()),
252        };
253    }
254
255    // Build the data frames.
256    let ric_df = df!(
257        "Offset (s)" => offset_s.clone(),
258        "Epoch" => epoch_str.clone(),
259        "RIC X (km)" => ric_x_km,
260        "RIC Y (km)" => ric_y_km,
261        "RIC Z (km)" => ric_z_km,
262        "RIC VX (km/s)" => ric_vx_km_s,
263        "RIC VY (km/s)" => ric_vy_km_s,
264        "RIC VZ (km/s)" => ric_vz_km_s,
265    )?;
266
267    println!("RIC difference at start\n{}", ric_df.head(Some(10)));
268    println!("RIC difference at end\n{}", ric_df.tail(Some(10)));
269
270    let aer_df = df!(
271        "Offset (s)" => offset_s.clone(),
272        "Epoch" => epoch_str.clone(),
273        "azimuth (deg)" => azimuth_deg,
274        "elevation (deg)" => elevation_deg,
275        "range (km)" => range_km,
276        "range rate (km/s)" => range_rate_km_s,
277    )?;
278
279    // Finally, let's see when the spacecraft is visible, assuming 15 degrees minimum elevation.
280    let mask = aer_df
281        .column("elevation (deg)")?
282        .gt(&Column::Scalar(ScalarColumn::new(
283            "elevation mask (deg)".into(),
284            Scalar::new(DataType::Float64, AnyValue::Float64(15.0)),
285            offset_s.len(),
286        )))?;
287    let cubesat_visible = aer_df.filter(&mask)?;
288
289    println!("{cubesat_visible}");
290
291    Ok(())
292}

pub fn set_sma_km(&mut self, new_sma_km: f64) -> Result<(), PhysicsError>

Mutates this orbit to change the SMA

:type new_sma_km: float :rtype: None

pub fn with_sma_km( &self, new_sma_km: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a new SMA

:type new_sma_km: float :rtype: Orbit

pub fn add_sma_km( &self, delta_sma_km: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a provided SMA added to the current one

:type delta_sma_km: float :rtype: Orbit

pub fn period(&self) -> Result<Duration, PhysicsError>

Returns the period in seconds

:rtype: Duration

pub fn ecc(&self) -> Result<f64, PhysicsError>

Returns the eccentricity (no unit)

:rtype: float

Examples found in repository?
examples/03_geo_analysis/drift.rs (line 122)
26fn main() -> Result<(), Box<dyn Error>> {
27    pel::init();
28    // Dynamics models require planetary constants and ephemerides to be defined.
29    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30    // This will automatically download the DE440s planetary ephemeris,
31    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
32    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
33    // planetary constants kernels.
34    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
35    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
36    // references to many functions.
37    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
38    // Define the orbit epoch
39    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
40
41    // Define the orbit.
42    // First we need to fetch the Earth J2000 from information from the Almanac.
43    // This allows the frame to include the gravitational parameters and the shape of the Earth,
44    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
45    // by loading a different set of planetary constants.
46    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
47
48    // Placing this GEO bird just above Colorado.
49    // In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
50    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
51    // Print in in Keplerian form.
52    println!("{orbit:x}");
53
54    let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
55    let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;
56
57    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
58    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
59    // models such as solar radiation pressure.
60
61    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
62    let sc = Spacecraft::builder()
63        .orbit(orbit)
64        .mass(Mass::from_dry_mass(9.60))
65        .srp(SRPData {
66            area_m2: 10e-4,
67            coeff_reflectivity: 1.1,
68        })
69        .build();
70    println!("{sc:x}");
71
72    // Set up the spacecraft dynamics.
73
74    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
75    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
76    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
77
78    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
79    // We're using the JGM3 model here, which is the default in GMAT.
80    let mut jgm3_meta = MetaFile {
81        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
82        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
83    };
84    // And let's download it if we don't have it yet.
85    jgm3_meta.process(true)?;
86
87    // Build the spherical harmonics.
88    // The harmonics must be computed in the body fixed frame.
89    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
90    let harmonics_21x21 = Harmonics::from_stor(
91        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
92        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
93    );
94
95    // Include the spherical harmonics into the orbital dynamics.
96    orbital_dyn.accel_models.push(harmonics_21x21);
97
98    // We define the solar radiation pressure, using the default solar flux and accounting only
99    // for the eclipsing caused by the Earth and Moon.
100    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
103    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
104    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
105
106    println!("{dynamics}");
107
108    // Finally, let's propagate this orbit to the same epoch as above.
109    // The first returned value is the spacecraft state at the final epoch.
110    // The second value is the full trajectory where the step size is variable step used by the propagator.
111    let (future_sc, trajectory) = Propagator::default(dynamics)
112        .with(sc, almanac.clone())
113        .until_epoch_with_traj(epoch + Unit::Century * 0.03)?;
114
115    println!("=== High fidelity propagation ===");
116    println!(
117        "SMA changed by {:.3} km",
118        orbit.sma_km()? - future_sc.orbit.sma_km()?
119    );
120    println!(
121        "ECC changed by {:.6}",
122        orbit.ecc()? - future_sc.orbit.ecc()?
123    );
124    println!(
125        "INC changed by {:.3e} deg",
126        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
127    );
128    println!(
129        "RAAN changed by {:.3} deg",
130        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
131    );
132    println!(
133        "AOP changed by {:.3} deg",
134        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
135    );
136    println!(
137        "TA changed by {:.3} deg",
138        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
139    );
140
141    // We also have access to the full trajectory throughout the propagation.
142    println!("{trajectory}");
143
144    println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");
145
146    // With the trajectory, let's build a few data products.
147
148    // 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.
149
150    let analysis_step = Unit::Minute * 5;
151
152    trajectory.to_parquet(
153        "./03_geo_hf_prop.parquet",
154        Some(vec![
155            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
156        ]),
157        ExportCfg::builder().step(analysis_step).build(),
158        almanac.clone(),
159    )?;
160
161    // 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.
162
163    // We iterate over the trajectory, grabbing a state every two minutes.
164    let mut offset_s = vec![];
165    let mut epoch_str = vec![];
166    let mut longitude_deg = vec![];
167    let mut latitude_deg = vec![];
168    let mut altitude_km = vec![];
169
170    for state in trajectory.every(analysis_step) {
171        // Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
172        // These define the GEO stationkeeping box.
173
174        let this_epoch = state.epoch();
175
176        offset_s.push((this_epoch - orbit.epoch).to_seconds());
177        epoch_str.push(this_epoch.to_isoformat());
178
179        let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
180        let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
181        longitude_deg.push(long_deg);
182        latitude_deg.push(lat_deg);
183        altitude_km.push(alt_km);
184    }
185
186    println!(
187        "Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
188        orig_long_deg - longitude_deg.last().unwrap()
189    );
190
191    println!(
192        "Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
193        orig_lat_deg - latitude_deg.last().unwrap()
194    );
195
196    println!(
197        "Altitude changed by {:.3} km -- Box is 30 km",
198        orig_alt_km - altitude_km.last().unwrap()
199    );
200
201    // Build the station keeping data frame.
202    let mut sk_df = df!(
203        "Offset (s)" => offset_s.clone(),
204        "Epoch (UTC)" => epoch_str.clone(),
205        "Longitude E-W (deg)" => longitude_deg,
206        "Latitude N-S (deg)" => latitude_deg,
207        "Altitude (km)" => altitude_km,
208
209    )?;
210
211    // Create a file to write the Parquet to
212    let file = File::create("./03_geo_lla.parquet").expect("Could not create file");
213
214    // Create a ParquetWriter and write the DataFrame to the file
215    ParquetWriter::new(file).finish(&mut sk_df)?;
216
217    Ok(())
218}
More examples
Hide additional examples
examples/01_orbit_prop/main.rs (line 71)
30fn main() -> Result<(), Box<dyn Error>> {
31    pel::init();
32    // Dynamics models require planetary constants and ephemerides to be defined.
33    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
34    // This will automatically download the DE440s planetary ephemeris,
35    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
36    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
37    // planetary constants kernels.
38    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
39    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
40    // references to many functions.
41    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
42    // Define the orbit epoch
43    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45    // Define the orbit.
46    // First we need to fetch the Earth J2000 from information from the Almanac.
47    // This allows the frame to include the gravitational parameters and the shape of the Earth,
48    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
49    // by loading a different set of planetary constants.
50    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
51
52    let orbit =
53        Orbit::try_keplerian_altitude(300.0, 0.015, 68.5, 65.2, 75.0, 0.0, epoch, earth_j2000)?;
54    // Print in in Keplerian form.
55    println!("{orbit:x}");
56
57    // There are two ways to propagate an orbit. We can make a quick approximation assuming only two-body
58    // motion. This is a useful first order approximation but it isn't used in real-world applications.
59
60    // This approach is a feature of ANISE.
61    let future_orbit_tb = orbit.at_epoch(epoch + Unit::Day * 3)?;
62    println!("{future_orbit_tb:x}");
63
64    // Two body propagation relies solely on Kepler's laws, so only the true anomaly will change.
65    println!(
66        "SMA changed by {:.3e} km",
67        orbit.sma_km()? - future_orbit_tb.sma_km()?
68    );
69    println!(
70        "ECC changed by {:.3e}",
71        orbit.ecc()? - future_orbit_tb.ecc()?
72    );
73    println!(
74        "INC changed by {:.3e} deg",
75        orbit.inc_deg()? - future_orbit_tb.inc_deg()?
76    );
77    println!(
78        "RAAN changed by {:.3e} deg",
79        orbit.raan_deg()? - future_orbit_tb.raan_deg()?
80    );
81    println!(
82        "AOP changed by {:.3e} deg",
83        orbit.aop_deg()? - future_orbit_tb.aop_deg()?
84    );
85    println!(
86        "TA changed by {:.3} deg",
87        orbit.ta_deg()? - future_orbit_tb.ta_deg()?
88    );
89
90    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
91    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
92    // models such as solar radiation pressure.
93
94    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
95    let sc = Spacecraft::builder()
96        .orbit(orbit)
97        .mass(Mass::from_dry_mass(9.60))
98        .srp(SRPData {
99            area_m2: 10e-4,
100            coeff_reflectivity: 1.1,
101        })
102        .build();
103    println!("{sc:x}");
104
105    // Set up the spacecraft dynamics.
106
107    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
108    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
109    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
110
111    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
112    // We're using the JGM3 model here, which is the default in GMAT.
113    let mut jgm3_meta = MetaFile {
114        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
115        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
116    };
117    // And let's download it if we don't have it yet.
118    jgm3_meta.process(true)?;
119
120    // Build the spherical harmonics.
121    // The harmonics must be computed in the body fixed frame.
122    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
123    let harmonics_21x21 = Harmonics::from_stor(
124        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
125        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
126    );
127
128    // Include the spherical harmonics into the orbital dynamics.
129    orbital_dyn.accel_models.push(harmonics_21x21);
130
131    // We define the solar radiation pressure, using the default solar flux and accounting only
132    // for the eclipsing caused by the Earth.
133    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
134
135    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
136    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
137    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
138
139    println!("{dynamics}");
140
141    // Finally, let's propagate this orbit to the same epoch as above.
142    // The first returned value is the spacecraft state at the final epoch.
143    // The second value is the full trajectory where the step size is variable step used by the propagator.
144    let (future_sc, trajectory) = Propagator::default(dynamics)
145        .with(sc, almanac.clone())
146        .until_epoch_with_traj(future_orbit_tb.epoch)?;
147
148    println!("=== High fidelity propagation ===");
149    println!(
150        "SMA changed by {:.3} km",
151        orbit.sma_km()? - future_sc.orbit.sma_km()?
152    );
153    println!(
154        "ECC changed by {:.6}",
155        orbit.ecc()? - future_sc.orbit.ecc()?
156    );
157    println!(
158        "INC changed by {:.3e} deg",
159        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
160    );
161    println!(
162        "RAAN changed by {:.3} deg",
163        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
164    );
165    println!(
166        "AOP changed by {:.3} deg",
167        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
168    );
169    println!(
170        "TA changed by {:.3} deg",
171        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
172    );
173
174    // We also have access to the full trajectory throughout the propagation.
175    println!("{trajectory}");
176
177    // With the trajectory, let's build a few data products.
178
179    // 1. Export the trajectory as a CCSDS OEM version 2.0 file and as a parquet file, which includes the Keplerian orbital elements.
180
181    trajectory.to_oem_file(
182        "./01_cubesat_hf_prop.oem",
183        ExportCfg::builder().step(Unit::Minute * 2).build(),
184    )?;
185
186    trajectory.to_parquet_with_cfg(
187        "./01_cubesat_hf_prop.parquet",
188        ExportCfg::builder().step(Unit::Minute * 2).build(),
189        almanac.clone(),
190    )?;
191
192    // 2. Compare the difference in the radial-intrack-crosstrack frame between the high fidelity
193    // and Keplerian propagation. The RIC frame is commonly used to compute the difference in position
194    // and velocity of different spacecraft.
195    // 3. Compute the azimuth, elevation, range, and range-rate data of that spacecraft as seen from Boulder, CO, USA.
196
197    let boulder_station = GroundStation::from_point(
198        "Boulder, CO, USA".to_string(),
199        40.014984,   // latitude in degrees
200        -105.270546, // longitude in degrees
201        1.6550,      // altitude in kilometers
202        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
203    );
204
205    // We iterate over the trajectory, grabbing a state every two minutes.
206    let mut offset_s = vec![];
207    let mut epoch_str = vec![];
208    let mut ric_x_km = vec![];
209    let mut ric_y_km = vec![];
210    let mut ric_z_km = vec![];
211    let mut ric_vx_km_s = vec![];
212    let mut ric_vy_km_s = vec![];
213    let mut ric_vz_km_s = vec![];
214
215    let mut azimuth_deg = vec![];
216    let mut elevation_deg = vec![];
217    let mut range_km = vec![];
218    let mut range_rate_km_s = vec![];
219    for state in trajectory.every(Unit::Minute * 2) {
220        // Try to compute the Keplerian/two body state just in time.
221        // This method occasionally fails to converge on an appropriate true anomaly
222        // from the mean anomaly. If that happens, we just skip this state.
223        // The high fidelity and Keplerian states diverge continuously, and we're curious
224        // about the divergence in this quick analysis.
225        let this_epoch = state.epoch();
226        match orbit.at_epoch(this_epoch) {
227            Ok(tb_then) => {
228                offset_s.push((this_epoch - orbit.epoch).to_seconds());
229                epoch_str.push(format!("{this_epoch}"));
230                // Compute the two body state just in time.
231                let ric = state.orbit.ric_difference(&tb_then)?;
232                ric_x_km.push(ric.radius_km.x);
233                ric_y_km.push(ric.radius_km.y);
234                ric_z_km.push(ric.radius_km.z);
235                ric_vx_km_s.push(ric.velocity_km_s.x);
236                ric_vy_km_s.push(ric.velocity_km_s.y);
237                ric_vz_km_s.push(ric.velocity_km_s.z);
238
239                // Compute the AER data for each state.
240                let aer = almanac.azimuth_elevation_range_sez(
241                    state.orbit,
242                    boulder_station.to_orbit(this_epoch, &almanac)?,
243                    None,
244                    None,
245                )?;
246                azimuth_deg.push(aer.azimuth_deg);
247                elevation_deg.push(aer.elevation_deg);
248                range_km.push(aer.range_km);
249                range_rate_km_s.push(aer.range_rate_km_s);
250            }
251            Err(e) => warn!("{} {e}", state.epoch()),
252        };
253    }
254
255    // Build the data frames.
256    let ric_df = df!(
257        "Offset (s)" => offset_s.clone(),
258        "Epoch" => epoch_str.clone(),
259        "RIC X (km)" => ric_x_km,
260        "RIC Y (km)" => ric_y_km,
261        "RIC Z (km)" => ric_z_km,
262        "RIC VX (km/s)" => ric_vx_km_s,
263        "RIC VY (km/s)" => ric_vy_km_s,
264        "RIC VZ (km/s)" => ric_vz_km_s,
265    )?;
266
267    println!("RIC difference at start\n{}", ric_df.head(Some(10)));
268    println!("RIC difference at end\n{}", ric_df.tail(Some(10)));
269
270    let aer_df = df!(
271        "Offset (s)" => offset_s.clone(),
272        "Epoch" => epoch_str.clone(),
273        "azimuth (deg)" => azimuth_deg,
274        "elevation (deg)" => elevation_deg,
275        "range (km)" => range_km,
276        "range rate (km/s)" => range_rate_km_s,
277    )?;
278
279    // Finally, let's see when the spacecraft is visible, assuming 15 degrees minimum elevation.
280    let mask = aer_df
281        .column("elevation (deg)")?
282        .gt(&Column::Scalar(ScalarColumn::new(
283            "elevation mask (deg)".into(),
284            Scalar::new(DataType::Float64, AnyValue::Float64(15.0)),
285            offset_s.len(),
286        )))?;
287    let cubesat_visible = aer_df.filter(&mask)?;
288
289    println!("{cubesat_visible}");
290
291    Ok(())
292}

pub fn set_ecc(&mut self, new_ecc: f64) -> Result<(), PhysicsError>

Mutates this orbit to change the ECC

:type new_ecc: float :rtype: None

pub fn with_ecc(&self, new_ecc: f64) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a new ECC

:type new_ecc: float :rtype: Orbit

pub fn add_ecc(&self, delta_ecc: f64) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a provided ECC added to the current one

:type delta_ecc: float :rtype: Orbit

pub fn inc_deg(&self) -> Result<f64, PhysicsError>

Returns the inclination in degrees

:rtype: float

Examples found in repository?
examples/03_geo_analysis/drift.rs (line 126)
26fn main() -> Result<(), Box<dyn Error>> {
27    pel::init();
28    // Dynamics models require planetary constants and ephemerides to be defined.
29    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30    // This will automatically download the DE440s planetary ephemeris,
31    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
32    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
33    // planetary constants kernels.
34    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
35    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
36    // references to many functions.
37    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
38    // Define the orbit epoch
39    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
40
41    // Define the orbit.
42    // First we need to fetch the Earth J2000 from information from the Almanac.
43    // This allows the frame to include the gravitational parameters and the shape of the Earth,
44    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
45    // by loading a different set of planetary constants.
46    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
47
48    // Placing this GEO bird just above Colorado.
49    // In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
50    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
51    // Print in in Keplerian form.
52    println!("{orbit:x}");
53
54    let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
55    let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;
56
57    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
58    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
59    // models such as solar radiation pressure.
60
61    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
62    let sc = Spacecraft::builder()
63        .orbit(orbit)
64        .mass(Mass::from_dry_mass(9.60))
65        .srp(SRPData {
66            area_m2: 10e-4,
67            coeff_reflectivity: 1.1,
68        })
69        .build();
70    println!("{sc:x}");
71
72    // Set up the spacecraft dynamics.
73
74    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
75    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
76    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
77
78    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
79    // We're using the JGM3 model here, which is the default in GMAT.
80    let mut jgm3_meta = MetaFile {
81        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
82        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
83    };
84    // And let's download it if we don't have it yet.
85    jgm3_meta.process(true)?;
86
87    // Build the spherical harmonics.
88    // The harmonics must be computed in the body fixed frame.
89    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
90    let harmonics_21x21 = Harmonics::from_stor(
91        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
92        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
93    );
94
95    // Include the spherical harmonics into the orbital dynamics.
96    orbital_dyn.accel_models.push(harmonics_21x21);
97
98    // We define the solar radiation pressure, using the default solar flux and accounting only
99    // for the eclipsing caused by the Earth and Moon.
100    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
103    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
104    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
105
106    println!("{dynamics}");
107
108    // Finally, let's propagate this orbit to the same epoch as above.
109    // The first returned value is the spacecraft state at the final epoch.
110    // The second value is the full trajectory where the step size is variable step used by the propagator.
111    let (future_sc, trajectory) = Propagator::default(dynamics)
112        .with(sc, almanac.clone())
113        .until_epoch_with_traj(epoch + Unit::Century * 0.03)?;
114
115    println!("=== High fidelity propagation ===");
116    println!(
117        "SMA changed by {:.3} km",
118        orbit.sma_km()? - future_sc.orbit.sma_km()?
119    );
120    println!(
121        "ECC changed by {:.6}",
122        orbit.ecc()? - future_sc.orbit.ecc()?
123    );
124    println!(
125        "INC changed by {:.3e} deg",
126        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
127    );
128    println!(
129        "RAAN changed by {:.3} deg",
130        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
131    );
132    println!(
133        "AOP changed by {:.3} deg",
134        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
135    );
136    println!(
137        "TA changed by {:.3} deg",
138        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
139    );
140
141    // We also have access to the full trajectory throughout the propagation.
142    println!("{trajectory}");
143
144    println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");
145
146    // With the trajectory, let's build a few data products.
147
148    // 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.
149
150    let analysis_step = Unit::Minute * 5;
151
152    trajectory.to_parquet(
153        "./03_geo_hf_prop.parquet",
154        Some(vec![
155            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
156        ]),
157        ExportCfg::builder().step(analysis_step).build(),
158        almanac.clone(),
159    )?;
160
161    // 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.
162
163    // We iterate over the trajectory, grabbing a state every two minutes.
164    let mut offset_s = vec![];
165    let mut epoch_str = vec![];
166    let mut longitude_deg = vec![];
167    let mut latitude_deg = vec![];
168    let mut altitude_km = vec![];
169
170    for state in trajectory.every(analysis_step) {
171        // Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
172        // These define the GEO stationkeeping box.
173
174        let this_epoch = state.epoch();
175
176        offset_s.push((this_epoch - orbit.epoch).to_seconds());
177        epoch_str.push(this_epoch.to_isoformat());
178
179        let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
180        let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
181        longitude_deg.push(long_deg);
182        latitude_deg.push(lat_deg);
183        altitude_km.push(alt_km);
184    }
185
186    println!(
187        "Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
188        orig_long_deg - longitude_deg.last().unwrap()
189    );
190
191    println!(
192        "Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
193        orig_lat_deg - latitude_deg.last().unwrap()
194    );
195
196    println!(
197        "Altitude changed by {:.3} km -- Box is 30 km",
198        orig_alt_km - altitude_km.last().unwrap()
199    );
200
201    // Build the station keeping data frame.
202    let mut sk_df = df!(
203        "Offset (s)" => offset_s.clone(),
204        "Epoch (UTC)" => epoch_str.clone(),
205        "Longitude E-W (deg)" => longitude_deg,
206        "Latitude N-S (deg)" => latitude_deg,
207        "Altitude (km)" => altitude_km,
208
209    )?;
210
211    // Create a file to write the Parquet to
212    let file = File::create("./03_geo_lla.parquet").expect("Could not create file");
213
214    // Create a ParquetWriter and write the DataFrame to the file
215    ParquetWriter::new(file).finish(&mut sk_df)?;
216
217    Ok(())
218}
More examples
Hide additional examples
examples/01_orbit_prop/main.rs (line 75)
30fn main() -> Result<(), Box<dyn Error>> {
31    pel::init();
32    // Dynamics models require planetary constants and ephemerides to be defined.
33    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
34    // This will automatically download the DE440s planetary ephemeris,
35    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
36    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
37    // planetary constants kernels.
38    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
39    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
40    // references to many functions.
41    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
42    // Define the orbit epoch
43    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45    // Define the orbit.
46    // First we need to fetch the Earth J2000 from information from the Almanac.
47    // This allows the frame to include the gravitational parameters and the shape of the Earth,
48    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
49    // by loading a different set of planetary constants.
50    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
51
52    let orbit =
53        Orbit::try_keplerian_altitude(300.0, 0.015, 68.5, 65.2, 75.0, 0.0, epoch, earth_j2000)?;
54    // Print in in Keplerian form.
55    println!("{orbit:x}");
56
57    // There are two ways to propagate an orbit. We can make a quick approximation assuming only two-body
58    // motion. This is a useful first order approximation but it isn't used in real-world applications.
59
60    // This approach is a feature of ANISE.
61    let future_orbit_tb = orbit.at_epoch(epoch + Unit::Day * 3)?;
62    println!("{future_orbit_tb:x}");
63
64    // Two body propagation relies solely on Kepler's laws, so only the true anomaly will change.
65    println!(
66        "SMA changed by {:.3e} km",
67        orbit.sma_km()? - future_orbit_tb.sma_km()?
68    );
69    println!(
70        "ECC changed by {:.3e}",
71        orbit.ecc()? - future_orbit_tb.ecc()?
72    );
73    println!(
74        "INC changed by {:.3e} deg",
75        orbit.inc_deg()? - future_orbit_tb.inc_deg()?
76    );
77    println!(
78        "RAAN changed by {:.3e} deg",
79        orbit.raan_deg()? - future_orbit_tb.raan_deg()?
80    );
81    println!(
82        "AOP changed by {:.3e} deg",
83        orbit.aop_deg()? - future_orbit_tb.aop_deg()?
84    );
85    println!(
86        "TA changed by {:.3} deg",
87        orbit.ta_deg()? - future_orbit_tb.ta_deg()?
88    );
89
90    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
91    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
92    // models such as solar radiation pressure.
93
94    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
95    let sc = Spacecraft::builder()
96        .orbit(orbit)
97        .mass(Mass::from_dry_mass(9.60))
98        .srp(SRPData {
99            area_m2: 10e-4,
100            coeff_reflectivity: 1.1,
101        })
102        .build();
103    println!("{sc:x}");
104
105    // Set up the spacecraft dynamics.
106
107    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
108    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
109    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
110
111    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
112    // We're using the JGM3 model here, which is the default in GMAT.
113    let mut jgm3_meta = MetaFile {
114        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
115        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
116    };
117    // And let's download it if we don't have it yet.
118    jgm3_meta.process(true)?;
119
120    // Build the spherical harmonics.
121    // The harmonics must be computed in the body fixed frame.
122    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
123    let harmonics_21x21 = Harmonics::from_stor(
124        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
125        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
126    );
127
128    // Include the spherical harmonics into the orbital dynamics.
129    orbital_dyn.accel_models.push(harmonics_21x21);
130
131    // We define the solar radiation pressure, using the default solar flux and accounting only
132    // for the eclipsing caused by the Earth.
133    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
134
135    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
136    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
137    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
138
139    println!("{dynamics}");
140
141    // Finally, let's propagate this orbit to the same epoch as above.
142    // The first returned value is the spacecraft state at the final epoch.
143    // The second value is the full trajectory where the step size is variable step used by the propagator.
144    let (future_sc, trajectory) = Propagator::default(dynamics)
145        .with(sc, almanac.clone())
146        .until_epoch_with_traj(future_orbit_tb.epoch)?;
147
148    println!("=== High fidelity propagation ===");
149    println!(
150        "SMA changed by {:.3} km",
151        orbit.sma_km()? - future_sc.orbit.sma_km()?
152    );
153    println!(
154        "ECC changed by {:.6}",
155        orbit.ecc()? - future_sc.orbit.ecc()?
156    );
157    println!(
158        "INC changed by {:.3e} deg",
159        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
160    );
161    println!(
162        "RAAN changed by {:.3} deg",
163        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
164    );
165    println!(
166        "AOP changed by {:.3} deg",
167        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
168    );
169    println!(
170        "TA changed by {:.3} deg",
171        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
172    );
173
174    // We also have access to the full trajectory throughout the propagation.
175    println!("{trajectory}");
176
177    // With the trajectory, let's build a few data products.
178
179    // 1. Export the trajectory as a CCSDS OEM version 2.0 file and as a parquet file, which includes the Keplerian orbital elements.
180
181    trajectory.to_oem_file(
182        "./01_cubesat_hf_prop.oem",
183        ExportCfg::builder().step(Unit::Minute * 2).build(),
184    )?;
185
186    trajectory.to_parquet_with_cfg(
187        "./01_cubesat_hf_prop.parquet",
188        ExportCfg::builder().step(Unit::Minute * 2).build(),
189        almanac.clone(),
190    )?;
191
192    // 2. Compare the difference in the radial-intrack-crosstrack frame between the high fidelity
193    // and Keplerian propagation. The RIC frame is commonly used to compute the difference in position
194    // and velocity of different spacecraft.
195    // 3. Compute the azimuth, elevation, range, and range-rate data of that spacecraft as seen from Boulder, CO, USA.
196
197    let boulder_station = GroundStation::from_point(
198        "Boulder, CO, USA".to_string(),
199        40.014984,   // latitude in degrees
200        -105.270546, // longitude in degrees
201        1.6550,      // altitude in kilometers
202        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
203    );
204
205    // We iterate over the trajectory, grabbing a state every two minutes.
206    let mut offset_s = vec![];
207    let mut epoch_str = vec![];
208    let mut ric_x_km = vec![];
209    let mut ric_y_km = vec![];
210    let mut ric_z_km = vec![];
211    let mut ric_vx_km_s = vec![];
212    let mut ric_vy_km_s = vec![];
213    let mut ric_vz_km_s = vec![];
214
215    let mut azimuth_deg = vec![];
216    let mut elevation_deg = vec![];
217    let mut range_km = vec![];
218    let mut range_rate_km_s = vec![];
219    for state in trajectory.every(Unit::Minute * 2) {
220        // Try to compute the Keplerian/two body state just in time.
221        // This method occasionally fails to converge on an appropriate true anomaly
222        // from the mean anomaly. If that happens, we just skip this state.
223        // The high fidelity and Keplerian states diverge continuously, and we're curious
224        // about the divergence in this quick analysis.
225        let this_epoch = state.epoch();
226        match orbit.at_epoch(this_epoch) {
227            Ok(tb_then) => {
228                offset_s.push((this_epoch - orbit.epoch).to_seconds());
229                epoch_str.push(format!("{this_epoch}"));
230                // Compute the two body state just in time.
231                let ric = state.orbit.ric_difference(&tb_then)?;
232                ric_x_km.push(ric.radius_km.x);
233                ric_y_km.push(ric.radius_km.y);
234                ric_z_km.push(ric.radius_km.z);
235                ric_vx_km_s.push(ric.velocity_km_s.x);
236                ric_vy_km_s.push(ric.velocity_km_s.y);
237                ric_vz_km_s.push(ric.velocity_km_s.z);
238
239                // Compute the AER data for each state.
240                let aer = almanac.azimuth_elevation_range_sez(
241                    state.orbit,
242                    boulder_station.to_orbit(this_epoch, &almanac)?,
243                    None,
244                    None,
245                )?;
246                azimuth_deg.push(aer.azimuth_deg);
247                elevation_deg.push(aer.elevation_deg);
248                range_km.push(aer.range_km);
249                range_rate_km_s.push(aer.range_rate_km_s);
250            }
251            Err(e) => warn!("{} {e}", state.epoch()),
252        };
253    }
254
255    // Build the data frames.
256    let ric_df = df!(
257        "Offset (s)" => offset_s.clone(),
258        "Epoch" => epoch_str.clone(),
259        "RIC X (km)" => ric_x_km,
260        "RIC Y (km)" => ric_y_km,
261        "RIC Z (km)" => ric_z_km,
262        "RIC VX (km/s)" => ric_vx_km_s,
263        "RIC VY (km/s)" => ric_vy_km_s,
264        "RIC VZ (km/s)" => ric_vz_km_s,
265    )?;
266
267    println!("RIC difference at start\n{}", ric_df.head(Some(10)));
268    println!("RIC difference at end\n{}", ric_df.tail(Some(10)));
269
270    let aer_df = df!(
271        "Offset (s)" => offset_s.clone(),
272        "Epoch" => epoch_str.clone(),
273        "azimuth (deg)" => azimuth_deg,
274        "elevation (deg)" => elevation_deg,
275        "range (km)" => range_km,
276        "range rate (km/s)" => range_rate_km_s,
277    )?;
278
279    // Finally, let's see when the spacecraft is visible, assuming 15 degrees minimum elevation.
280    let mask = aer_df
281        .column("elevation (deg)")?
282        .gt(&Column::Scalar(ScalarColumn::new(
283            "elevation mask (deg)".into(),
284            Scalar::new(DataType::Float64, AnyValue::Float64(15.0)),
285            offset_s.len(),
286        )))?;
287    let cubesat_visible = aer_df.filter(&mask)?;
288
289    println!("{cubesat_visible}");
290
291    Ok(())
292}

pub fn set_inc_deg(&mut self, new_inc_deg: f64) -> Result<(), PhysicsError>

Mutates this orbit to change the INC

:type new_inc_deg: float :rtype: None

pub fn with_inc_deg( &self, new_inc_deg: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a new INC

:type new_inc_deg: float :rtype: Orbit

pub fn add_inc_deg( &self, delta_inc_deg: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a provided INC added to the current one

:type delta_inc_deg: float :rtype: None

pub fn aop_deg(&self) -> Result<f64, PhysicsError>

Returns the argument of periapsis in degrees

:rtype: float

Examples found in repository?
examples/03_geo_analysis/drift.rs (line 134)
26fn main() -> Result<(), Box<dyn Error>> {
27    pel::init();
28    // Dynamics models require planetary constants and ephemerides to be defined.
29    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30    // This will automatically download the DE440s planetary ephemeris,
31    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
32    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
33    // planetary constants kernels.
34    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
35    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
36    // references to many functions.
37    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
38    // Define the orbit epoch
39    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
40
41    // Define the orbit.
42    // First we need to fetch the Earth J2000 from information from the Almanac.
43    // This allows the frame to include the gravitational parameters and the shape of the Earth,
44    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
45    // by loading a different set of planetary constants.
46    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
47
48    // Placing this GEO bird just above Colorado.
49    // In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
50    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
51    // Print in in Keplerian form.
52    println!("{orbit:x}");
53
54    let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
55    let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;
56
57    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
58    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
59    // models such as solar radiation pressure.
60
61    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
62    let sc = Spacecraft::builder()
63        .orbit(orbit)
64        .mass(Mass::from_dry_mass(9.60))
65        .srp(SRPData {
66            area_m2: 10e-4,
67            coeff_reflectivity: 1.1,
68        })
69        .build();
70    println!("{sc:x}");
71
72    // Set up the spacecraft dynamics.
73
74    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
75    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
76    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
77
78    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
79    // We're using the JGM3 model here, which is the default in GMAT.
80    let mut jgm3_meta = MetaFile {
81        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
82        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
83    };
84    // And let's download it if we don't have it yet.
85    jgm3_meta.process(true)?;
86
87    // Build the spherical harmonics.
88    // The harmonics must be computed in the body fixed frame.
89    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
90    let harmonics_21x21 = Harmonics::from_stor(
91        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
92        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
93    );
94
95    // Include the spherical harmonics into the orbital dynamics.
96    orbital_dyn.accel_models.push(harmonics_21x21);
97
98    // We define the solar radiation pressure, using the default solar flux and accounting only
99    // for the eclipsing caused by the Earth and Moon.
100    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
103    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
104    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
105
106    println!("{dynamics}");
107
108    // Finally, let's propagate this orbit to the same epoch as above.
109    // The first returned value is the spacecraft state at the final epoch.
110    // The second value is the full trajectory where the step size is variable step used by the propagator.
111    let (future_sc, trajectory) = Propagator::default(dynamics)
112        .with(sc, almanac.clone())
113        .until_epoch_with_traj(epoch + Unit::Century * 0.03)?;
114
115    println!("=== High fidelity propagation ===");
116    println!(
117        "SMA changed by {:.3} km",
118        orbit.sma_km()? - future_sc.orbit.sma_km()?
119    );
120    println!(
121        "ECC changed by {:.6}",
122        orbit.ecc()? - future_sc.orbit.ecc()?
123    );
124    println!(
125        "INC changed by {:.3e} deg",
126        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
127    );
128    println!(
129        "RAAN changed by {:.3} deg",
130        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
131    );
132    println!(
133        "AOP changed by {:.3} deg",
134        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
135    );
136    println!(
137        "TA changed by {:.3} deg",
138        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
139    );
140
141    // We also have access to the full trajectory throughout the propagation.
142    println!("{trajectory}");
143
144    println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");
145
146    // With the trajectory, let's build a few data products.
147
148    // 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.
149
150    let analysis_step = Unit::Minute * 5;
151
152    trajectory.to_parquet(
153        "./03_geo_hf_prop.parquet",
154        Some(vec![
155            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
156        ]),
157        ExportCfg::builder().step(analysis_step).build(),
158        almanac.clone(),
159    )?;
160
161    // 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.
162
163    // We iterate over the trajectory, grabbing a state every two minutes.
164    let mut offset_s = vec![];
165    let mut epoch_str = vec![];
166    let mut longitude_deg = vec![];
167    let mut latitude_deg = vec![];
168    let mut altitude_km = vec![];
169
170    for state in trajectory.every(analysis_step) {
171        // Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
172        // These define the GEO stationkeeping box.
173
174        let this_epoch = state.epoch();
175
176        offset_s.push((this_epoch - orbit.epoch).to_seconds());
177        epoch_str.push(this_epoch.to_isoformat());
178
179        let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
180        let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
181        longitude_deg.push(long_deg);
182        latitude_deg.push(lat_deg);
183        altitude_km.push(alt_km);
184    }
185
186    println!(
187        "Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
188        orig_long_deg - longitude_deg.last().unwrap()
189    );
190
191    println!(
192        "Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
193        orig_lat_deg - latitude_deg.last().unwrap()
194    );
195
196    println!(
197        "Altitude changed by {:.3} km -- Box is 30 km",
198        orig_alt_km - altitude_km.last().unwrap()
199    );
200
201    // Build the station keeping data frame.
202    let mut sk_df = df!(
203        "Offset (s)" => offset_s.clone(),
204        "Epoch (UTC)" => epoch_str.clone(),
205        "Longitude E-W (deg)" => longitude_deg,
206        "Latitude N-S (deg)" => latitude_deg,
207        "Altitude (km)" => altitude_km,
208
209    )?;
210
211    // Create a file to write the Parquet to
212    let file = File::create("./03_geo_lla.parquet").expect("Could not create file");
213
214    // Create a ParquetWriter and write the DataFrame to the file
215    ParquetWriter::new(file).finish(&mut sk_df)?;
216
217    Ok(())
218}
More examples
Hide additional examples
examples/01_orbit_prop/main.rs (line 83)
30fn main() -> Result<(), Box<dyn Error>> {
31    pel::init();
32    // Dynamics models require planetary constants and ephemerides to be defined.
33    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
34    // This will automatically download the DE440s planetary ephemeris,
35    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
36    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
37    // planetary constants kernels.
38    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
39    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
40    // references to many functions.
41    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
42    // Define the orbit epoch
43    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45    // Define the orbit.
46    // First we need to fetch the Earth J2000 from information from the Almanac.
47    // This allows the frame to include the gravitational parameters and the shape of the Earth,
48    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
49    // by loading a different set of planetary constants.
50    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
51
52    let orbit =
53        Orbit::try_keplerian_altitude(300.0, 0.015, 68.5, 65.2, 75.0, 0.0, epoch, earth_j2000)?;
54    // Print in in Keplerian form.
55    println!("{orbit:x}");
56
57    // There are two ways to propagate an orbit. We can make a quick approximation assuming only two-body
58    // motion. This is a useful first order approximation but it isn't used in real-world applications.
59
60    // This approach is a feature of ANISE.
61    let future_orbit_tb = orbit.at_epoch(epoch + Unit::Day * 3)?;
62    println!("{future_orbit_tb:x}");
63
64    // Two body propagation relies solely on Kepler's laws, so only the true anomaly will change.
65    println!(
66        "SMA changed by {:.3e} km",
67        orbit.sma_km()? - future_orbit_tb.sma_km()?
68    );
69    println!(
70        "ECC changed by {:.3e}",
71        orbit.ecc()? - future_orbit_tb.ecc()?
72    );
73    println!(
74        "INC changed by {:.3e} deg",
75        orbit.inc_deg()? - future_orbit_tb.inc_deg()?
76    );
77    println!(
78        "RAAN changed by {:.3e} deg",
79        orbit.raan_deg()? - future_orbit_tb.raan_deg()?
80    );
81    println!(
82        "AOP changed by {:.3e} deg",
83        orbit.aop_deg()? - future_orbit_tb.aop_deg()?
84    );
85    println!(
86        "TA changed by {:.3} deg",
87        orbit.ta_deg()? - future_orbit_tb.ta_deg()?
88    );
89
90    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
91    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
92    // models such as solar radiation pressure.
93
94    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
95    let sc = Spacecraft::builder()
96        .orbit(orbit)
97        .mass(Mass::from_dry_mass(9.60))
98        .srp(SRPData {
99            area_m2: 10e-4,
100            coeff_reflectivity: 1.1,
101        })
102        .build();
103    println!("{sc:x}");
104
105    // Set up the spacecraft dynamics.
106
107    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
108    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
109    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
110
111    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
112    // We're using the JGM3 model here, which is the default in GMAT.
113    let mut jgm3_meta = MetaFile {
114        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
115        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
116    };
117    // And let's download it if we don't have it yet.
118    jgm3_meta.process(true)?;
119
120    // Build the spherical harmonics.
121    // The harmonics must be computed in the body fixed frame.
122    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
123    let harmonics_21x21 = Harmonics::from_stor(
124        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
125        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
126    );
127
128    // Include the spherical harmonics into the orbital dynamics.
129    orbital_dyn.accel_models.push(harmonics_21x21);
130
131    // We define the solar radiation pressure, using the default solar flux and accounting only
132    // for the eclipsing caused by the Earth.
133    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
134
135    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
136    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
137    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
138
139    println!("{dynamics}");
140
141    // Finally, let's propagate this orbit to the same epoch as above.
142    // The first returned value is the spacecraft state at the final epoch.
143    // The second value is the full trajectory where the step size is variable step used by the propagator.
144    let (future_sc, trajectory) = Propagator::default(dynamics)
145        .with(sc, almanac.clone())
146        .until_epoch_with_traj(future_orbit_tb.epoch)?;
147
148    println!("=== High fidelity propagation ===");
149    println!(
150        "SMA changed by {:.3} km",
151        orbit.sma_km()? - future_sc.orbit.sma_km()?
152    );
153    println!(
154        "ECC changed by {:.6}",
155        orbit.ecc()? - future_sc.orbit.ecc()?
156    );
157    println!(
158        "INC changed by {:.3e} deg",
159        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
160    );
161    println!(
162        "RAAN changed by {:.3} deg",
163        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
164    );
165    println!(
166        "AOP changed by {:.3} deg",
167        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
168    );
169    println!(
170        "TA changed by {:.3} deg",
171        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
172    );
173
174    // We also have access to the full trajectory throughout the propagation.
175    println!("{trajectory}");
176
177    // With the trajectory, let's build a few data products.
178
179    // 1. Export the trajectory as a CCSDS OEM version 2.0 file and as a parquet file, which includes the Keplerian orbital elements.
180
181    trajectory.to_oem_file(
182        "./01_cubesat_hf_prop.oem",
183        ExportCfg::builder().step(Unit::Minute * 2).build(),
184    )?;
185
186    trajectory.to_parquet_with_cfg(
187        "./01_cubesat_hf_prop.parquet",
188        ExportCfg::builder().step(Unit::Minute * 2).build(),
189        almanac.clone(),
190    )?;
191
192    // 2. Compare the difference in the radial-intrack-crosstrack frame between the high fidelity
193    // and Keplerian propagation. The RIC frame is commonly used to compute the difference in position
194    // and velocity of different spacecraft.
195    // 3. Compute the azimuth, elevation, range, and range-rate data of that spacecraft as seen from Boulder, CO, USA.
196
197    let boulder_station = GroundStation::from_point(
198        "Boulder, CO, USA".to_string(),
199        40.014984,   // latitude in degrees
200        -105.270546, // longitude in degrees
201        1.6550,      // altitude in kilometers
202        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
203    );
204
205    // We iterate over the trajectory, grabbing a state every two minutes.
206    let mut offset_s = vec![];
207    let mut epoch_str = vec![];
208    let mut ric_x_km = vec![];
209    let mut ric_y_km = vec![];
210    let mut ric_z_km = vec![];
211    let mut ric_vx_km_s = vec![];
212    let mut ric_vy_km_s = vec![];
213    let mut ric_vz_km_s = vec![];
214
215    let mut azimuth_deg = vec![];
216    let mut elevation_deg = vec![];
217    let mut range_km = vec![];
218    let mut range_rate_km_s = vec![];
219    for state in trajectory.every(Unit::Minute * 2) {
220        // Try to compute the Keplerian/two body state just in time.
221        // This method occasionally fails to converge on an appropriate true anomaly
222        // from the mean anomaly. If that happens, we just skip this state.
223        // The high fidelity and Keplerian states diverge continuously, and we're curious
224        // about the divergence in this quick analysis.
225        let this_epoch = state.epoch();
226        match orbit.at_epoch(this_epoch) {
227            Ok(tb_then) => {
228                offset_s.push((this_epoch - orbit.epoch).to_seconds());
229                epoch_str.push(format!("{this_epoch}"));
230                // Compute the two body state just in time.
231                let ric = state.orbit.ric_difference(&tb_then)?;
232                ric_x_km.push(ric.radius_km.x);
233                ric_y_km.push(ric.radius_km.y);
234                ric_z_km.push(ric.radius_km.z);
235                ric_vx_km_s.push(ric.velocity_km_s.x);
236                ric_vy_km_s.push(ric.velocity_km_s.y);
237                ric_vz_km_s.push(ric.velocity_km_s.z);
238
239                // Compute the AER data for each state.
240                let aer = almanac.azimuth_elevation_range_sez(
241                    state.orbit,
242                    boulder_station.to_orbit(this_epoch, &almanac)?,
243                    None,
244                    None,
245                )?;
246                azimuth_deg.push(aer.azimuth_deg);
247                elevation_deg.push(aer.elevation_deg);
248                range_km.push(aer.range_km);
249                range_rate_km_s.push(aer.range_rate_km_s);
250            }
251            Err(e) => warn!("{} {e}", state.epoch()),
252        };
253    }
254
255    // Build the data frames.
256    let ric_df = df!(
257        "Offset (s)" => offset_s.clone(),
258        "Epoch" => epoch_str.clone(),
259        "RIC X (km)" => ric_x_km,
260        "RIC Y (km)" => ric_y_km,
261        "RIC Z (km)" => ric_z_km,
262        "RIC VX (km/s)" => ric_vx_km_s,
263        "RIC VY (km/s)" => ric_vy_km_s,
264        "RIC VZ (km/s)" => ric_vz_km_s,
265    )?;
266
267    println!("RIC difference at start\n{}", ric_df.head(Some(10)));
268    println!("RIC difference at end\n{}", ric_df.tail(Some(10)));
269
270    let aer_df = df!(
271        "Offset (s)" => offset_s.clone(),
272        "Epoch" => epoch_str.clone(),
273        "azimuth (deg)" => azimuth_deg,
274        "elevation (deg)" => elevation_deg,
275        "range (km)" => range_km,
276        "range rate (km/s)" => range_rate_km_s,
277    )?;
278
279    // Finally, let's see when the spacecraft is visible, assuming 15 degrees minimum elevation.
280    let mask = aer_df
281        .column("elevation (deg)")?
282        .gt(&Column::Scalar(ScalarColumn::new(
283            "elevation mask (deg)".into(),
284            Scalar::new(DataType::Float64, AnyValue::Float64(15.0)),
285            offset_s.len(),
286        )))?;
287    let cubesat_visible = aer_df.filter(&mask)?;
288
289    println!("{cubesat_visible}");
290
291    Ok(())
292}

pub fn set_aop_deg(&mut self, new_aop_deg: f64) -> Result<(), PhysicsError>

Mutates this orbit to change the AOP

:type new_aop_deg: float :rtype: None

pub fn with_aop_deg( &self, new_aop_deg: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a new AOP

:type new_aop_deg: float :rtype: Orbit

pub fn add_aop_deg( &self, delta_aop_deg: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a provided AOP added to the current one

:type delta_aop_deg: float :rtype: Orbit

pub fn raan_deg(&self) -> Result<f64, PhysicsError>

Returns the right ascension of the ascending node in degrees

:rtype: float

Examples found in repository?
examples/03_geo_analysis/drift.rs (line 130)
26fn main() -> Result<(), Box<dyn Error>> {
27    pel::init();
28    // Dynamics models require planetary constants and ephemerides to be defined.
29    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30    // This will automatically download the DE440s planetary ephemeris,
31    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
32    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
33    // planetary constants kernels.
34    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
35    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
36    // references to many functions.
37    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
38    // Define the orbit epoch
39    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
40
41    // Define the orbit.
42    // First we need to fetch the Earth J2000 from information from the Almanac.
43    // This allows the frame to include the gravitational parameters and the shape of the Earth,
44    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
45    // by loading a different set of planetary constants.
46    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
47
48    // Placing this GEO bird just above Colorado.
49    // In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
50    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
51    // Print in in Keplerian form.
52    println!("{orbit:x}");
53
54    let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
55    let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;
56
57    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
58    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
59    // models such as solar radiation pressure.
60
61    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
62    let sc = Spacecraft::builder()
63        .orbit(orbit)
64        .mass(Mass::from_dry_mass(9.60))
65        .srp(SRPData {
66            area_m2: 10e-4,
67            coeff_reflectivity: 1.1,
68        })
69        .build();
70    println!("{sc:x}");
71
72    // Set up the spacecraft dynamics.
73
74    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
75    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
76    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
77
78    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
79    // We're using the JGM3 model here, which is the default in GMAT.
80    let mut jgm3_meta = MetaFile {
81        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
82        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
83    };
84    // And let's download it if we don't have it yet.
85    jgm3_meta.process(true)?;
86
87    // Build the spherical harmonics.
88    // The harmonics must be computed in the body fixed frame.
89    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
90    let harmonics_21x21 = Harmonics::from_stor(
91        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
92        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
93    );
94
95    // Include the spherical harmonics into the orbital dynamics.
96    orbital_dyn.accel_models.push(harmonics_21x21);
97
98    // We define the solar radiation pressure, using the default solar flux and accounting only
99    // for the eclipsing caused by the Earth and Moon.
100    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
103    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
104    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
105
106    println!("{dynamics}");
107
108    // Finally, let's propagate this orbit to the same epoch as above.
109    // The first returned value is the spacecraft state at the final epoch.
110    // The second value is the full trajectory where the step size is variable step used by the propagator.
111    let (future_sc, trajectory) = Propagator::default(dynamics)
112        .with(sc, almanac.clone())
113        .until_epoch_with_traj(epoch + Unit::Century * 0.03)?;
114
115    println!("=== High fidelity propagation ===");
116    println!(
117        "SMA changed by {:.3} km",
118        orbit.sma_km()? - future_sc.orbit.sma_km()?
119    );
120    println!(
121        "ECC changed by {:.6}",
122        orbit.ecc()? - future_sc.orbit.ecc()?
123    );
124    println!(
125        "INC changed by {:.3e} deg",
126        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
127    );
128    println!(
129        "RAAN changed by {:.3} deg",
130        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
131    );
132    println!(
133        "AOP changed by {:.3} deg",
134        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
135    );
136    println!(
137        "TA changed by {:.3} deg",
138        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
139    );
140
141    // We also have access to the full trajectory throughout the propagation.
142    println!("{trajectory}");
143
144    println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");
145
146    // With the trajectory, let's build a few data products.
147
148    // 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.
149
150    let analysis_step = Unit::Minute * 5;
151
152    trajectory.to_parquet(
153        "./03_geo_hf_prop.parquet",
154        Some(vec![
155            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
156        ]),
157        ExportCfg::builder().step(analysis_step).build(),
158        almanac.clone(),
159    )?;
160
161    // 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.
162
163    // We iterate over the trajectory, grabbing a state every two minutes.
164    let mut offset_s = vec![];
165    let mut epoch_str = vec![];
166    let mut longitude_deg = vec![];
167    let mut latitude_deg = vec![];
168    let mut altitude_km = vec![];
169
170    for state in trajectory.every(analysis_step) {
171        // Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
172        // These define the GEO stationkeeping box.
173
174        let this_epoch = state.epoch();
175
176        offset_s.push((this_epoch - orbit.epoch).to_seconds());
177        epoch_str.push(this_epoch.to_isoformat());
178
179        let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
180        let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
181        longitude_deg.push(long_deg);
182        latitude_deg.push(lat_deg);
183        altitude_km.push(alt_km);
184    }
185
186    println!(
187        "Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
188        orig_long_deg - longitude_deg.last().unwrap()
189    );
190
191    println!(
192        "Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
193        orig_lat_deg - latitude_deg.last().unwrap()
194    );
195
196    println!(
197        "Altitude changed by {:.3} km -- Box is 30 km",
198        orig_alt_km - altitude_km.last().unwrap()
199    );
200
201    // Build the station keeping data frame.
202    let mut sk_df = df!(
203        "Offset (s)" => offset_s.clone(),
204        "Epoch (UTC)" => epoch_str.clone(),
205        "Longitude E-W (deg)" => longitude_deg,
206        "Latitude N-S (deg)" => latitude_deg,
207        "Altitude (km)" => altitude_km,
208
209    )?;
210
211    // Create a file to write the Parquet to
212    let file = File::create("./03_geo_lla.parquet").expect("Could not create file");
213
214    // Create a ParquetWriter and write the DataFrame to the file
215    ParquetWriter::new(file).finish(&mut sk_df)?;
216
217    Ok(())
218}
More examples
Hide additional examples
examples/01_orbit_prop/main.rs (line 79)
30fn main() -> Result<(), Box<dyn Error>> {
31    pel::init();
32    // Dynamics models require planetary constants and ephemerides to be defined.
33    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
34    // This will automatically download the DE440s planetary ephemeris,
35    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
36    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
37    // planetary constants kernels.
38    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
39    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
40    // references to many functions.
41    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
42    // Define the orbit epoch
43    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45    // Define the orbit.
46    // First we need to fetch the Earth J2000 from information from the Almanac.
47    // This allows the frame to include the gravitational parameters and the shape of the Earth,
48    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
49    // by loading a different set of planetary constants.
50    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
51
52    let orbit =
53        Orbit::try_keplerian_altitude(300.0, 0.015, 68.5, 65.2, 75.0, 0.0, epoch, earth_j2000)?;
54    // Print in in Keplerian form.
55    println!("{orbit:x}");
56
57    // There are two ways to propagate an orbit. We can make a quick approximation assuming only two-body
58    // motion. This is a useful first order approximation but it isn't used in real-world applications.
59
60    // This approach is a feature of ANISE.
61    let future_orbit_tb = orbit.at_epoch(epoch + Unit::Day * 3)?;
62    println!("{future_orbit_tb:x}");
63
64    // Two body propagation relies solely on Kepler's laws, so only the true anomaly will change.
65    println!(
66        "SMA changed by {:.3e} km",
67        orbit.sma_km()? - future_orbit_tb.sma_km()?
68    );
69    println!(
70        "ECC changed by {:.3e}",
71        orbit.ecc()? - future_orbit_tb.ecc()?
72    );
73    println!(
74        "INC changed by {:.3e} deg",
75        orbit.inc_deg()? - future_orbit_tb.inc_deg()?
76    );
77    println!(
78        "RAAN changed by {:.3e} deg",
79        orbit.raan_deg()? - future_orbit_tb.raan_deg()?
80    );
81    println!(
82        "AOP changed by {:.3e} deg",
83        orbit.aop_deg()? - future_orbit_tb.aop_deg()?
84    );
85    println!(
86        "TA changed by {:.3} deg",
87        orbit.ta_deg()? - future_orbit_tb.ta_deg()?
88    );
89
90    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
91    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
92    // models such as solar radiation pressure.
93
94    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
95    let sc = Spacecraft::builder()
96        .orbit(orbit)
97        .mass(Mass::from_dry_mass(9.60))
98        .srp(SRPData {
99            area_m2: 10e-4,
100            coeff_reflectivity: 1.1,
101        })
102        .build();
103    println!("{sc:x}");
104
105    // Set up the spacecraft dynamics.
106
107    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
108    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
109    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
110
111    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
112    // We're using the JGM3 model here, which is the default in GMAT.
113    let mut jgm3_meta = MetaFile {
114        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
115        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
116    };
117    // And let's download it if we don't have it yet.
118    jgm3_meta.process(true)?;
119
120    // Build the spherical harmonics.
121    // The harmonics must be computed in the body fixed frame.
122    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
123    let harmonics_21x21 = Harmonics::from_stor(
124        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
125        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
126    );
127
128    // Include the spherical harmonics into the orbital dynamics.
129    orbital_dyn.accel_models.push(harmonics_21x21);
130
131    // We define the solar radiation pressure, using the default solar flux and accounting only
132    // for the eclipsing caused by the Earth.
133    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
134
135    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
136    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
137    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
138
139    println!("{dynamics}");
140
141    // Finally, let's propagate this orbit to the same epoch as above.
142    // The first returned value is the spacecraft state at the final epoch.
143    // The second value is the full trajectory where the step size is variable step used by the propagator.
144    let (future_sc, trajectory) = Propagator::default(dynamics)
145        .with(sc, almanac.clone())
146        .until_epoch_with_traj(future_orbit_tb.epoch)?;
147
148    println!("=== High fidelity propagation ===");
149    println!(
150        "SMA changed by {:.3} km",
151        orbit.sma_km()? - future_sc.orbit.sma_km()?
152    );
153    println!(
154        "ECC changed by {:.6}",
155        orbit.ecc()? - future_sc.orbit.ecc()?
156    );
157    println!(
158        "INC changed by {:.3e} deg",
159        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
160    );
161    println!(
162        "RAAN changed by {:.3} deg",
163        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
164    );
165    println!(
166        "AOP changed by {:.3} deg",
167        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
168    );
169    println!(
170        "TA changed by {:.3} deg",
171        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
172    );
173
174    // We also have access to the full trajectory throughout the propagation.
175    println!("{trajectory}");
176
177    // With the trajectory, let's build a few data products.
178
179    // 1. Export the trajectory as a CCSDS OEM version 2.0 file and as a parquet file, which includes the Keplerian orbital elements.
180
181    trajectory.to_oem_file(
182        "./01_cubesat_hf_prop.oem",
183        ExportCfg::builder().step(Unit::Minute * 2).build(),
184    )?;
185
186    trajectory.to_parquet_with_cfg(
187        "./01_cubesat_hf_prop.parquet",
188        ExportCfg::builder().step(Unit::Minute * 2).build(),
189        almanac.clone(),
190    )?;
191
192    // 2. Compare the difference in the radial-intrack-crosstrack frame between the high fidelity
193    // and Keplerian propagation. The RIC frame is commonly used to compute the difference in position
194    // and velocity of different spacecraft.
195    // 3. Compute the azimuth, elevation, range, and range-rate data of that spacecraft as seen from Boulder, CO, USA.
196
197    let boulder_station = GroundStation::from_point(
198        "Boulder, CO, USA".to_string(),
199        40.014984,   // latitude in degrees
200        -105.270546, // longitude in degrees
201        1.6550,      // altitude in kilometers
202        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
203    );
204
205    // We iterate over the trajectory, grabbing a state every two minutes.
206    let mut offset_s = vec![];
207    let mut epoch_str = vec![];
208    let mut ric_x_km = vec![];
209    let mut ric_y_km = vec![];
210    let mut ric_z_km = vec![];
211    let mut ric_vx_km_s = vec![];
212    let mut ric_vy_km_s = vec![];
213    let mut ric_vz_km_s = vec![];
214
215    let mut azimuth_deg = vec![];
216    let mut elevation_deg = vec![];
217    let mut range_km = vec![];
218    let mut range_rate_km_s = vec![];
219    for state in trajectory.every(Unit::Minute * 2) {
220        // Try to compute the Keplerian/two body state just in time.
221        // This method occasionally fails to converge on an appropriate true anomaly
222        // from the mean anomaly. If that happens, we just skip this state.
223        // The high fidelity and Keplerian states diverge continuously, and we're curious
224        // about the divergence in this quick analysis.
225        let this_epoch = state.epoch();
226        match orbit.at_epoch(this_epoch) {
227            Ok(tb_then) => {
228                offset_s.push((this_epoch - orbit.epoch).to_seconds());
229                epoch_str.push(format!("{this_epoch}"));
230                // Compute the two body state just in time.
231                let ric = state.orbit.ric_difference(&tb_then)?;
232                ric_x_km.push(ric.radius_km.x);
233                ric_y_km.push(ric.radius_km.y);
234                ric_z_km.push(ric.radius_km.z);
235                ric_vx_km_s.push(ric.velocity_km_s.x);
236                ric_vy_km_s.push(ric.velocity_km_s.y);
237                ric_vz_km_s.push(ric.velocity_km_s.z);
238
239                // Compute the AER data for each state.
240                let aer = almanac.azimuth_elevation_range_sez(
241                    state.orbit,
242                    boulder_station.to_orbit(this_epoch, &almanac)?,
243                    None,
244                    None,
245                )?;
246                azimuth_deg.push(aer.azimuth_deg);
247                elevation_deg.push(aer.elevation_deg);
248                range_km.push(aer.range_km);
249                range_rate_km_s.push(aer.range_rate_km_s);
250            }
251            Err(e) => warn!("{} {e}", state.epoch()),
252        };
253    }
254
255    // Build the data frames.
256    let ric_df = df!(
257        "Offset (s)" => offset_s.clone(),
258        "Epoch" => epoch_str.clone(),
259        "RIC X (km)" => ric_x_km,
260        "RIC Y (km)" => ric_y_km,
261        "RIC Z (km)" => ric_z_km,
262        "RIC VX (km/s)" => ric_vx_km_s,
263        "RIC VY (km/s)" => ric_vy_km_s,
264        "RIC VZ (km/s)" => ric_vz_km_s,
265    )?;
266
267    println!("RIC difference at start\n{}", ric_df.head(Some(10)));
268    println!("RIC difference at end\n{}", ric_df.tail(Some(10)));
269
270    let aer_df = df!(
271        "Offset (s)" => offset_s.clone(),
272        "Epoch" => epoch_str.clone(),
273        "azimuth (deg)" => azimuth_deg,
274        "elevation (deg)" => elevation_deg,
275        "range (km)" => range_km,
276        "range rate (km/s)" => range_rate_km_s,
277    )?;
278
279    // Finally, let's see when the spacecraft is visible, assuming 15 degrees minimum elevation.
280    let mask = aer_df
281        .column("elevation (deg)")?
282        .gt(&Column::Scalar(ScalarColumn::new(
283            "elevation mask (deg)".into(),
284            Scalar::new(DataType::Float64, AnyValue::Float64(15.0)),
285            offset_s.len(),
286        )))?;
287    let cubesat_visible = aer_df.filter(&mask)?;
288
289    println!("{cubesat_visible}");
290
291    Ok(())
292}

pub fn set_raan_deg(&mut self, new_raan_deg: f64) -> Result<(), PhysicsError>

Mutates this orbit to change the RAAN

:type new_raan_deg: float :rtype: None

pub fn with_raan_deg( &self, new_raan_deg: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a new RAAN

:type new_raan_deg: float :rtype: Orbit

pub fn add_raan_deg( &self, delta_raan_deg: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a provided RAAN added to the current one

:type delta_raan_deg: float :rtype: Orbit

pub fn ta_deg(&self) -> Result<f64, PhysicsError>

Returns the true anomaly in degrees between 0 and 360.0

NOTE: This function will emit a warning stating that the TA should be avoided if in a very near circular orbit Code from https://github.com/ChristopherRabotin/GMAT/blob/80bde040e12946a61dae90d9fc3538f16df34190/src/gmatutil/util/StateConversionUtil.cpp#L6835

LIMITATION: For an orbit whose true anomaly is (very nearly) 0.0 or 180.0, this function may return either 0.0 or 180.0 with a very small time increment. This is due to the precision of the cosine calculation: if the arccosine calculation is out of bounds, the sign of the cosine of the true anomaly is used to determine whether the true anomaly should be 0.0 or 180.0. In other words, there is an ambiguity in the computation in the true anomaly exactly at 180.0 and 0.0.

:rtype: float

Examples found in repository?
examples/03_geo_analysis/drift.rs (line 138)
26fn main() -> Result<(), Box<dyn Error>> {
27    pel::init();
28    // Dynamics models require planetary constants and ephemerides to be defined.
29    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30    // This will automatically download the DE440s planetary ephemeris,
31    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
32    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
33    // planetary constants kernels.
34    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
35    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
36    // references to many functions.
37    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
38    // Define the orbit epoch
39    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
40
41    // Define the orbit.
42    // First we need to fetch the Earth J2000 from information from the Almanac.
43    // This allows the frame to include the gravitational parameters and the shape of the Earth,
44    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
45    // by loading a different set of planetary constants.
46    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
47
48    // Placing this GEO bird just above Colorado.
49    // In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
50    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
51    // Print in in Keplerian form.
52    println!("{orbit:x}");
53
54    let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
55    let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;
56
57    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
58    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
59    // models such as solar radiation pressure.
60
61    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
62    let sc = Spacecraft::builder()
63        .orbit(orbit)
64        .mass(Mass::from_dry_mass(9.60))
65        .srp(SRPData {
66            area_m2: 10e-4,
67            coeff_reflectivity: 1.1,
68        })
69        .build();
70    println!("{sc:x}");
71
72    // Set up the spacecraft dynamics.
73
74    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
75    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
76    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
77
78    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
79    // We're using the JGM3 model here, which is the default in GMAT.
80    let mut jgm3_meta = MetaFile {
81        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
82        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
83    };
84    // And let's download it if we don't have it yet.
85    jgm3_meta.process(true)?;
86
87    // Build the spherical harmonics.
88    // The harmonics must be computed in the body fixed frame.
89    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
90    let harmonics_21x21 = Harmonics::from_stor(
91        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
92        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
93    );
94
95    // Include the spherical harmonics into the orbital dynamics.
96    orbital_dyn.accel_models.push(harmonics_21x21);
97
98    // We define the solar radiation pressure, using the default solar flux and accounting only
99    // for the eclipsing caused by the Earth and Moon.
100    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
103    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
104    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
105
106    println!("{dynamics}");
107
108    // Finally, let's propagate this orbit to the same epoch as above.
109    // The first returned value is the spacecraft state at the final epoch.
110    // The second value is the full trajectory where the step size is variable step used by the propagator.
111    let (future_sc, trajectory) = Propagator::default(dynamics)
112        .with(sc, almanac.clone())
113        .until_epoch_with_traj(epoch + Unit::Century * 0.03)?;
114
115    println!("=== High fidelity propagation ===");
116    println!(
117        "SMA changed by {:.3} km",
118        orbit.sma_km()? - future_sc.orbit.sma_km()?
119    );
120    println!(
121        "ECC changed by {:.6}",
122        orbit.ecc()? - future_sc.orbit.ecc()?
123    );
124    println!(
125        "INC changed by {:.3e} deg",
126        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
127    );
128    println!(
129        "RAAN changed by {:.3} deg",
130        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
131    );
132    println!(
133        "AOP changed by {:.3} deg",
134        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
135    );
136    println!(
137        "TA changed by {:.3} deg",
138        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
139    );
140
141    // We also have access to the full trajectory throughout the propagation.
142    println!("{trajectory}");
143
144    println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");
145
146    // With the trajectory, let's build a few data products.
147
148    // 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.
149
150    let analysis_step = Unit::Minute * 5;
151
152    trajectory.to_parquet(
153        "./03_geo_hf_prop.parquet",
154        Some(vec![
155            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
156        ]),
157        ExportCfg::builder().step(analysis_step).build(),
158        almanac.clone(),
159    )?;
160
161    // 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.
162
163    // We iterate over the trajectory, grabbing a state every two minutes.
164    let mut offset_s = vec![];
165    let mut epoch_str = vec![];
166    let mut longitude_deg = vec![];
167    let mut latitude_deg = vec![];
168    let mut altitude_km = vec![];
169
170    for state in trajectory.every(analysis_step) {
171        // Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
172        // These define the GEO stationkeeping box.
173
174        let this_epoch = state.epoch();
175
176        offset_s.push((this_epoch - orbit.epoch).to_seconds());
177        epoch_str.push(this_epoch.to_isoformat());
178
179        let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
180        let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
181        longitude_deg.push(long_deg);
182        latitude_deg.push(lat_deg);
183        altitude_km.push(alt_km);
184    }
185
186    println!(
187        "Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
188        orig_long_deg - longitude_deg.last().unwrap()
189    );
190
191    println!(
192        "Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
193        orig_lat_deg - latitude_deg.last().unwrap()
194    );
195
196    println!(
197        "Altitude changed by {:.3} km -- Box is 30 km",
198        orig_alt_km - altitude_km.last().unwrap()
199    );
200
201    // Build the station keeping data frame.
202    let mut sk_df = df!(
203        "Offset (s)" => offset_s.clone(),
204        "Epoch (UTC)" => epoch_str.clone(),
205        "Longitude E-W (deg)" => longitude_deg,
206        "Latitude N-S (deg)" => latitude_deg,
207        "Altitude (km)" => altitude_km,
208
209    )?;
210
211    // Create a file to write the Parquet to
212    let file = File::create("./03_geo_lla.parquet").expect("Could not create file");
213
214    // Create a ParquetWriter and write the DataFrame to the file
215    ParquetWriter::new(file).finish(&mut sk_df)?;
216
217    Ok(())
218}
More examples
Hide additional examples
examples/01_orbit_prop/main.rs (line 87)
30fn main() -> Result<(), Box<dyn Error>> {
31    pel::init();
32    // Dynamics models require planetary constants and ephemerides to be defined.
33    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
34    // This will automatically download the DE440s planetary ephemeris,
35    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
36    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
37    // planetary constants kernels.
38    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
39    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
40    // references to many functions.
41    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
42    // Define the orbit epoch
43    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45    // Define the orbit.
46    // First we need to fetch the Earth J2000 from information from the Almanac.
47    // This allows the frame to include the gravitational parameters and the shape of the Earth,
48    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
49    // by loading a different set of planetary constants.
50    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
51
52    let orbit =
53        Orbit::try_keplerian_altitude(300.0, 0.015, 68.5, 65.2, 75.0, 0.0, epoch, earth_j2000)?;
54    // Print in in Keplerian form.
55    println!("{orbit:x}");
56
57    // There are two ways to propagate an orbit. We can make a quick approximation assuming only two-body
58    // motion. This is a useful first order approximation but it isn't used in real-world applications.
59
60    // This approach is a feature of ANISE.
61    let future_orbit_tb = orbit.at_epoch(epoch + Unit::Day * 3)?;
62    println!("{future_orbit_tb:x}");
63
64    // Two body propagation relies solely on Kepler's laws, so only the true anomaly will change.
65    println!(
66        "SMA changed by {:.3e} km",
67        orbit.sma_km()? - future_orbit_tb.sma_km()?
68    );
69    println!(
70        "ECC changed by {:.3e}",
71        orbit.ecc()? - future_orbit_tb.ecc()?
72    );
73    println!(
74        "INC changed by {:.3e} deg",
75        orbit.inc_deg()? - future_orbit_tb.inc_deg()?
76    );
77    println!(
78        "RAAN changed by {:.3e} deg",
79        orbit.raan_deg()? - future_orbit_tb.raan_deg()?
80    );
81    println!(
82        "AOP changed by {:.3e} deg",
83        orbit.aop_deg()? - future_orbit_tb.aop_deg()?
84    );
85    println!(
86        "TA changed by {:.3} deg",
87        orbit.ta_deg()? - future_orbit_tb.ta_deg()?
88    );
89
90    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
91    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
92    // models such as solar radiation pressure.
93
94    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
95    let sc = Spacecraft::builder()
96        .orbit(orbit)
97        .mass(Mass::from_dry_mass(9.60))
98        .srp(SRPData {
99            area_m2: 10e-4,
100            coeff_reflectivity: 1.1,
101        })
102        .build();
103    println!("{sc:x}");
104
105    // Set up the spacecraft dynamics.
106
107    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
108    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
109    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
110
111    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
112    // We're using the JGM3 model here, which is the default in GMAT.
113    let mut jgm3_meta = MetaFile {
114        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
115        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
116    };
117    // And let's download it if we don't have it yet.
118    jgm3_meta.process(true)?;
119
120    // Build the spherical harmonics.
121    // The harmonics must be computed in the body fixed frame.
122    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
123    let harmonics_21x21 = Harmonics::from_stor(
124        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
125        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
126    );
127
128    // Include the spherical harmonics into the orbital dynamics.
129    orbital_dyn.accel_models.push(harmonics_21x21);
130
131    // We define the solar radiation pressure, using the default solar flux and accounting only
132    // for the eclipsing caused by the Earth.
133    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
134
135    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
136    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
137    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
138
139    println!("{dynamics}");
140
141    // Finally, let's propagate this orbit to the same epoch as above.
142    // The first returned value is the spacecraft state at the final epoch.
143    // The second value is the full trajectory where the step size is variable step used by the propagator.
144    let (future_sc, trajectory) = Propagator::default(dynamics)
145        .with(sc, almanac.clone())
146        .until_epoch_with_traj(future_orbit_tb.epoch)?;
147
148    println!("=== High fidelity propagation ===");
149    println!(
150        "SMA changed by {:.3} km",
151        orbit.sma_km()? - future_sc.orbit.sma_km()?
152    );
153    println!(
154        "ECC changed by {:.6}",
155        orbit.ecc()? - future_sc.orbit.ecc()?
156    );
157    println!(
158        "INC changed by {:.3e} deg",
159        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
160    );
161    println!(
162        "RAAN changed by {:.3} deg",
163        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
164    );
165    println!(
166        "AOP changed by {:.3} deg",
167        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
168    );
169    println!(
170        "TA changed by {:.3} deg",
171        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
172    );
173
174    // We also have access to the full trajectory throughout the propagation.
175    println!("{trajectory}");
176
177    // With the trajectory, let's build a few data products.
178
179    // 1. Export the trajectory as a CCSDS OEM version 2.0 file and as a parquet file, which includes the Keplerian orbital elements.
180
181    trajectory.to_oem_file(
182        "./01_cubesat_hf_prop.oem",
183        ExportCfg::builder().step(Unit::Minute * 2).build(),
184    )?;
185
186    trajectory.to_parquet_with_cfg(
187        "./01_cubesat_hf_prop.parquet",
188        ExportCfg::builder().step(Unit::Minute * 2).build(),
189        almanac.clone(),
190    )?;
191
192    // 2. Compare the difference in the radial-intrack-crosstrack frame between the high fidelity
193    // and Keplerian propagation. The RIC frame is commonly used to compute the difference in position
194    // and velocity of different spacecraft.
195    // 3. Compute the azimuth, elevation, range, and range-rate data of that spacecraft as seen from Boulder, CO, USA.
196
197    let boulder_station = GroundStation::from_point(
198        "Boulder, CO, USA".to_string(),
199        40.014984,   // latitude in degrees
200        -105.270546, // longitude in degrees
201        1.6550,      // altitude in kilometers
202        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
203    );
204
205    // We iterate over the trajectory, grabbing a state every two minutes.
206    let mut offset_s = vec![];
207    let mut epoch_str = vec![];
208    let mut ric_x_km = vec![];
209    let mut ric_y_km = vec![];
210    let mut ric_z_km = vec![];
211    let mut ric_vx_km_s = vec![];
212    let mut ric_vy_km_s = vec![];
213    let mut ric_vz_km_s = vec![];
214
215    let mut azimuth_deg = vec![];
216    let mut elevation_deg = vec![];
217    let mut range_km = vec![];
218    let mut range_rate_km_s = vec![];
219    for state in trajectory.every(Unit::Minute * 2) {
220        // Try to compute the Keplerian/two body state just in time.
221        // This method occasionally fails to converge on an appropriate true anomaly
222        // from the mean anomaly. If that happens, we just skip this state.
223        // The high fidelity and Keplerian states diverge continuously, and we're curious
224        // about the divergence in this quick analysis.
225        let this_epoch = state.epoch();
226        match orbit.at_epoch(this_epoch) {
227            Ok(tb_then) => {
228                offset_s.push((this_epoch - orbit.epoch).to_seconds());
229                epoch_str.push(format!("{this_epoch}"));
230                // Compute the two body state just in time.
231                let ric = state.orbit.ric_difference(&tb_then)?;
232                ric_x_km.push(ric.radius_km.x);
233                ric_y_km.push(ric.radius_km.y);
234                ric_z_km.push(ric.radius_km.z);
235                ric_vx_km_s.push(ric.velocity_km_s.x);
236                ric_vy_km_s.push(ric.velocity_km_s.y);
237                ric_vz_km_s.push(ric.velocity_km_s.z);
238
239                // Compute the AER data for each state.
240                let aer = almanac.azimuth_elevation_range_sez(
241                    state.orbit,
242                    boulder_station.to_orbit(this_epoch, &almanac)?,
243                    None,
244                    None,
245                )?;
246                azimuth_deg.push(aer.azimuth_deg);
247                elevation_deg.push(aer.elevation_deg);
248                range_km.push(aer.range_km);
249                range_rate_km_s.push(aer.range_rate_km_s);
250            }
251            Err(e) => warn!("{} {e}", state.epoch()),
252        };
253    }
254
255    // Build the data frames.
256    let ric_df = df!(
257        "Offset (s)" => offset_s.clone(),
258        "Epoch" => epoch_str.clone(),
259        "RIC X (km)" => ric_x_km,
260        "RIC Y (km)" => ric_y_km,
261        "RIC Z (km)" => ric_z_km,
262        "RIC VX (km/s)" => ric_vx_km_s,
263        "RIC VY (km/s)" => ric_vy_km_s,
264        "RIC VZ (km/s)" => ric_vz_km_s,
265    )?;
266
267    println!("RIC difference at start\n{}", ric_df.head(Some(10)));
268    println!("RIC difference at end\n{}", ric_df.tail(Some(10)));
269
270    let aer_df = df!(
271        "Offset (s)" => offset_s.clone(),
272        "Epoch" => epoch_str.clone(),
273        "azimuth (deg)" => azimuth_deg,
274        "elevation (deg)" => elevation_deg,
275        "range (km)" => range_km,
276        "range rate (km/s)" => range_rate_km_s,
277    )?;
278
279    // Finally, let's see when the spacecraft is visible, assuming 15 degrees minimum elevation.
280    let mask = aer_df
281        .column("elevation (deg)")?
282        .gt(&Column::Scalar(ScalarColumn::new(
283            "elevation mask (deg)".into(),
284            Scalar::new(DataType::Float64, AnyValue::Float64(15.0)),
285            offset_s.len(),
286        )))?;
287    let cubesat_visible = aer_df.filter(&mask)?;
288
289    println!("{cubesat_visible}");
290
291    Ok(())
292}

pub fn set_ta_deg(&mut self, new_ta_deg: f64) -> Result<(), PhysicsError>

Mutates this orbit to change the TA

:type new_ta_deg: float :rtype: None

pub fn ta_dot_deg_s(&self) -> Result<f64, PhysicsError>

Returns the time derivative of the true anomaly computed as the 360.0 degrees divided by the orbital period (in seconds).

:rtype: float

pub fn with_ta_deg( &self, new_ta_deg: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a new TA

:type new_ta_deg: float :rtype: Orbit

pub fn add_ta_deg( &self, delta_ta_deg: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of the state with a provided TA added to the current one

:type delta_ta_deg: float :rtype: Orbit

pub fn with_apoapsis_periapsis_km( &self, new_ra_km: f64, new_rp_km: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of this state with the provided apoasis and periapsis

:type new_ra_km: float :type new_rp_km: float :rtype: Orbit

pub fn add_apoapsis_periapsis_km( &self, delta_ra_km: f64, delta_rp_km: f64, ) -> Result<CartesianState, PhysicsError>

Returns a copy of this state with the provided apoasis and periapsis added to the current values

:type delta_ra_km: float :type delta_rp_km: float :rtype: Orbit

pub fn tlong_deg(&self) -> Result<f64, PhysicsError>

Returns the true longitude in degrees

:rtype: float

pub fn aol_deg(&self) -> Result<f64, PhysicsError>

Returns the argument of latitude in degrees

NOTE: If the orbit is near circular, the AoL will be computed from the true longitude instead of relying on the ill-defined true anomaly.

:rtype: float

pub fn periapsis_km(&self) -> Result<f64, PhysicsError>

Returns the radius of periapsis (or perigee around Earth), in kilometers.

:rtype: float

pub fn apoapsis_km(&self) -> Result<f64, PhysicsError>

Returns the radius of apoapsis (or apogee around Earth), in kilometers.

:rtype: float

pub fn ea_deg(&self) -> Result<f64, PhysicsError>

Returns the eccentric anomaly in degrees

This is a conversion from GMAT’s StateConversionUtil::TrueToEccentricAnomaly

:rtype: float

pub fn fpa_deg(&self) -> Result<f64, PhysicsError>

Returns the flight path angle in degrees

:rtype: float

pub fn ma_deg(&self) -> Result<f64, PhysicsError>

Returns the mean anomaly in degrees

This is a conversion from GMAT’s StateConversionUtil::TrueToMeanAnomaly

:rtype: float

pub fn semi_parameter_km(&self) -> Result<f64, PhysicsError>

Returns the semi parameter (or semilatus rectum)

:rtype: float

pub fn is_brouwer_short_valid(&self) -> Result<bool, PhysicsError>

Returns whether this state satisfies the requirement to compute the Mean Brouwer Short orbital element set.

This is a conversion from GMAT’s StateConversionUtil::CartesianToBrouwerMeanShort. The details are at the log level info. NOTE: Mean Brouwer Short are only defined around Earth. However, nyx does not check the main celestial body around which the state is defined (GMAT does perform this verification).

:rtype: bool

pub fn right_ascension_deg(&self) -> f64

Returns the right ascension of this orbit in degrees

:rtype: float

pub fn declination_deg(&self) -> f64

Returns the declination of this orbit in degrees

:rtype: float

pub fn semi_minor_axis_km(&self) -> Result<f64, PhysicsError>

Returns the semi minor axis in km, includes code for a hyperbolic orbit

:rtype: float

pub fn velocity_declination_deg(&self) -> f64

Returns the velocity declination of this orbit in degrees

:rtype: float

pub fn c3_km2_s2(&self) -> Result<f64, PhysicsError>

Returns the $C_3$ of this orbit in km^2/s^2

:rtype: float

pub fn vinf_periapsis_km( &self, turn_angle_degrees: f64, ) -> Result<f64, PhysicsError>

Returns the radius of periapse in kilometers for the provided turn angle of this hyperbolic orbit. Returns an error if the orbit is not hyperbolic.

:type turn_angle_degrees: float :rtype: float

pub fn vinf_turn_angle_deg( &self, periapsis_km: f64, ) -> Result<f64, PhysicsError>

Returns the turn angle in degrees for the provided radius of periapse passage of this hyperbolic orbit Returns an error if the orbit is not hyperbolic.

:type periapsis_km: float :rtype: float

pub fn hyperbolic_anomaly_deg(&self) -> Result<f64, PhysicsError>

Returns the hyperbolic anomaly in degrees between 0 and 360.0 Returns an error if the orbit is not hyperbolic.

:rtype: float

pub fn at_epoch(&self, new_epoch: Epoch) -> Result<CartesianState, PhysicsError>

Adjusts the true anomaly of this orbit using the mean anomaly.

§Astrodynamics note

This is not a true propagation of the orbit. This is akin to a two body propagation ONLY without any other force models applied. Use Nyx for high fidelity propagation.

:type new_epoch: Epoch :rtype: Orbit

Examples found in repository?
examples/01_orbit_prop/main.rs (line 61)
30fn main() -> Result<(), Box<dyn Error>> {
31    pel::init();
32    // Dynamics models require planetary constants and ephemerides to be defined.
33    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
34    // This will automatically download the DE440s planetary ephemeris,
35    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
36    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
37    // planetary constants kernels.
38    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
39    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
40    // references to many functions.
41    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
42    // Define the orbit epoch
43    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45    // Define the orbit.
46    // First we need to fetch the Earth J2000 from information from the Almanac.
47    // This allows the frame to include the gravitational parameters and the shape of the Earth,
48    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
49    // by loading a different set of planetary constants.
50    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
51
52    let orbit =
53        Orbit::try_keplerian_altitude(300.0, 0.015, 68.5, 65.2, 75.0, 0.0, epoch, earth_j2000)?;
54    // Print in in Keplerian form.
55    println!("{orbit:x}");
56
57    // There are two ways to propagate an orbit. We can make a quick approximation assuming only two-body
58    // motion. This is a useful first order approximation but it isn't used in real-world applications.
59
60    // This approach is a feature of ANISE.
61    let future_orbit_tb = orbit.at_epoch(epoch + Unit::Day * 3)?;
62    println!("{future_orbit_tb:x}");
63
64    // Two body propagation relies solely on Kepler's laws, so only the true anomaly will change.
65    println!(
66        "SMA changed by {:.3e} km",
67        orbit.sma_km()? - future_orbit_tb.sma_km()?
68    );
69    println!(
70        "ECC changed by {:.3e}",
71        orbit.ecc()? - future_orbit_tb.ecc()?
72    );
73    println!(
74        "INC changed by {:.3e} deg",
75        orbit.inc_deg()? - future_orbit_tb.inc_deg()?
76    );
77    println!(
78        "RAAN changed by {:.3e} deg",
79        orbit.raan_deg()? - future_orbit_tb.raan_deg()?
80    );
81    println!(
82        "AOP changed by {:.3e} deg",
83        orbit.aop_deg()? - future_orbit_tb.aop_deg()?
84    );
85    println!(
86        "TA changed by {:.3} deg",
87        orbit.ta_deg()? - future_orbit_tb.ta_deg()?
88    );
89
90    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
91    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
92    // models such as solar radiation pressure.
93
94    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
95    let sc = Spacecraft::builder()
96        .orbit(orbit)
97        .mass(Mass::from_dry_mass(9.60))
98        .srp(SRPData {
99            area_m2: 10e-4,
100            coeff_reflectivity: 1.1,
101        })
102        .build();
103    println!("{sc:x}");
104
105    // Set up the spacecraft dynamics.
106
107    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
108    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
109    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
110
111    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
112    // We're using the JGM3 model here, which is the default in GMAT.
113    let mut jgm3_meta = MetaFile {
114        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
115        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
116    };
117    // And let's download it if we don't have it yet.
118    jgm3_meta.process(true)?;
119
120    // Build the spherical harmonics.
121    // The harmonics must be computed in the body fixed frame.
122    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
123    let harmonics_21x21 = Harmonics::from_stor(
124        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
125        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
126    );
127
128    // Include the spherical harmonics into the orbital dynamics.
129    orbital_dyn.accel_models.push(harmonics_21x21);
130
131    // We define the solar radiation pressure, using the default solar flux and accounting only
132    // for the eclipsing caused by the Earth.
133    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
134
135    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
136    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
137    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
138
139    println!("{dynamics}");
140
141    // Finally, let's propagate this orbit to the same epoch as above.
142    // The first returned value is the spacecraft state at the final epoch.
143    // The second value is the full trajectory where the step size is variable step used by the propagator.
144    let (future_sc, trajectory) = Propagator::default(dynamics)
145        .with(sc, almanac.clone())
146        .until_epoch_with_traj(future_orbit_tb.epoch)?;
147
148    println!("=== High fidelity propagation ===");
149    println!(
150        "SMA changed by {:.3} km",
151        orbit.sma_km()? - future_sc.orbit.sma_km()?
152    );
153    println!(
154        "ECC changed by {:.6}",
155        orbit.ecc()? - future_sc.orbit.ecc()?
156    );
157    println!(
158        "INC changed by {:.3e} deg",
159        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
160    );
161    println!(
162        "RAAN changed by {:.3} deg",
163        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
164    );
165    println!(
166        "AOP changed by {:.3} deg",
167        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
168    );
169    println!(
170        "TA changed by {:.3} deg",
171        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
172    );
173
174    // We also have access to the full trajectory throughout the propagation.
175    println!("{trajectory}");
176
177    // With the trajectory, let's build a few data products.
178
179    // 1. Export the trajectory as a CCSDS OEM version 2.0 file and as a parquet file, which includes the Keplerian orbital elements.
180
181    trajectory.to_oem_file(
182        "./01_cubesat_hf_prop.oem",
183        ExportCfg::builder().step(Unit::Minute * 2).build(),
184    )?;
185
186    trajectory.to_parquet_with_cfg(
187        "./01_cubesat_hf_prop.parquet",
188        ExportCfg::builder().step(Unit::Minute * 2).build(),
189        almanac.clone(),
190    )?;
191
192    // 2. Compare the difference in the radial-intrack-crosstrack frame between the high fidelity
193    // and Keplerian propagation. The RIC frame is commonly used to compute the difference in position
194    // and velocity of different spacecraft.
195    // 3. Compute the azimuth, elevation, range, and range-rate data of that spacecraft as seen from Boulder, CO, USA.
196
197    let boulder_station = GroundStation::from_point(
198        "Boulder, CO, USA".to_string(),
199        40.014984,   // latitude in degrees
200        -105.270546, // longitude in degrees
201        1.6550,      // altitude in kilometers
202        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
203    );
204
205    // We iterate over the trajectory, grabbing a state every two minutes.
206    let mut offset_s = vec![];
207    let mut epoch_str = vec![];
208    let mut ric_x_km = vec![];
209    let mut ric_y_km = vec![];
210    let mut ric_z_km = vec![];
211    let mut ric_vx_km_s = vec![];
212    let mut ric_vy_km_s = vec![];
213    let mut ric_vz_km_s = vec![];
214
215    let mut azimuth_deg = vec![];
216    let mut elevation_deg = vec![];
217    let mut range_km = vec![];
218    let mut range_rate_km_s = vec![];
219    for state in trajectory.every(Unit::Minute * 2) {
220        // Try to compute the Keplerian/two body state just in time.
221        // This method occasionally fails to converge on an appropriate true anomaly
222        // from the mean anomaly. If that happens, we just skip this state.
223        // The high fidelity and Keplerian states diverge continuously, and we're curious
224        // about the divergence in this quick analysis.
225        let this_epoch = state.epoch();
226        match orbit.at_epoch(this_epoch) {
227            Ok(tb_then) => {
228                offset_s.push((this_epoch - orbit.epoch).to_seconds());
229                epoch_str.push(format!("{this_epoch}"));
230                // Compute the two body state just in time.
231                let ric = state.orbit.ric_difference(&tb_then)?;
232                ric_x_km.push(ric.radius_km.x);
233                ric_y_km.push(ric.radius_km.y);
234                ric_z_km.push(ric.radius_km.z);
235                ric_vx_km_s.push(ric.velocity_km_s.x);
236                ric_vy_km_s.push(ric.velocity_km_s.y);
237                ric_vz_km_s.push(ric.velocity_km_s.z);
238
239                // Compute the AER data for each state.
240                let aer = almanac.azimuth_elevation_range_sez(
241                    state.orbit,
242                    boulder_station.to_orbit(this_epoch, &almanac)?,
243                    None,
244                    None,
245                )?;
246                azimuth_deg.push(aer.azimuth_deg);
247                elevation_deg.push(aer.elevation_deg);
248                range_km.push(aer.range_km);
249                range_rate_km_s.push(aer.range_rate_km_s);
250            }
251            Err(e) => warn!("{} {e}", state.epoch()),
252        };
253    }
254
255    // Build the data frames.
256    let ric_df = df!(
257        "Offset (s)" => offset_s.clone(),
258        "Epoch" => epoch_str.clone(),
259        "RIC X (km)" => ric_x_km,
260        "RIC Y (km)" => ric_y_km,
261        "RIC Z (km)" => ric_z_km,
262        "RIC VX (km/s)" => ric_vx_km_s,
263        "RIC VY (km/s)" => ric_vy_km_s,
264        "RIC VZ (km/s)" => ric_vz_km_s,
265    )?;
266
267    println!("RIC difference at start\n{}", ric_df.head(Some(10)));
268    println!("RIC difference at end\n{}", ric_df.tail(Some(10)));
269
270    let aer_df = df!(
271        "Offset (s)" => offset_s.clone(),
272        "Epoch" => epoch_str.clone(),
273        "azimuth (deg)" => azimuth_deg,
274        "elevation (deg)" => elevation_deg,
275        "range (km)" => range_km,
276        "range rate (km/s)" => range_rate_km_s,
277    )?;
278
279    // Finally, let's see when the spacecraft is visible, assuming 15 degrees minimum elevation.
280    let mask = aer_df
281        .column("elevation (deg)")?
282        .gt(&Column::Scalar(ScalarColumn::new(
283            "elevation mask (deg)".into(),
284            Scalar::new(DataType::Float64, AnyValue::Float64(15.0)),
285            offset_s.len(),
286        )))?;
287    let cubesat_visible = aer_df.filter(&mask)?;
288
289    println!("{cubesat_visible}");
290
291    Ok(())
292}

pub fn ric_difference( &self, other: &CartesianState, ) -> Result<CartesianState, PhysicsError>

Returns a Cartesian state representing the RIC difference between self and other, in position and velocity (with transport theorem). Refer to dcm_from_ric_to_inertial for details on the RIC frame.

§Algorithm
  1. Compute the RIC DCM of self
  2. Rotate self into the RIC frame
  3. Rotation other into the RIC frame
  4. Compute the difference between these two states
  5. Strip the astrodynamical information from the frame, enabling only computations from CartesianState

:type other: Orbit :rtype: Orbit

Examples found in repository?
examples/04_lro_od/main.rs (line 159)
33fn main() -> Result<(), Box<dyn Error>> {
34    pel::init();
35
36    // ====================== //
37    // === ALMANAC SET UP === //
38    // ====================== //
39
40    // Dynamics models require planetary constants and ephemerides to be defined.
41    // Let's start by grabbing those by using ANISE's MetaAlmanac.
42
43    let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
44        .iter()
45        .collect();
46
47    let meta = data_folder.join("lro-dynamics.dhall");
48
49    // Load this ephem in the general Almanac we're using for this analysis.
50    let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
51        .map_err(Box::new)?
52        .process(true)
53        .map_err(Box::new)?;
54
55    let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
56    moon_pc.mu_km3_s2 = 4902.74987;
57    almanac.planetary_data.set_by_id(MOON, moon_pc)?;
58
59    let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
60    earth_pc.mu_km3_s2 = 398600.436;
61    almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
62
63    // Save this new kernel for reuse.
64    // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
65    almanac
66        .planetary_data
67        .save_as(&data_folder.join("lro-specific.pca"), true)?;
68
69    // Lock the almanac (an Arc is a read only structure).
70    let almanac = Arc::new(almanac);
71
72    // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
73    // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
74    // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
75    // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
76    let lro_frame = Frame::from_ephem_j2000(-85);
77
78    // To build the trajectory we need to provide a spacecraft template.
79    let sc_template = Spacecraft::builder()
80        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
81        .srp(SRPData {
82            // SRP configuration is arbitrary, but we will be estimating it anyway.
83            area_m2: 3.9 * 2.7,
84            coeff_reflectivity: 0.96,
85        })
86        .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
87        .build();
88    // Now we can build the trajectory from the BSP file.
89    // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
90    let traj_as_flown = Traj::from_bsp(
91        lro_frame,
92        MOON_J2000,
93        almanac.clone(),
94        sc_template,
95        5.seconds(),
96        Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
97        Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
98        Aberration::LT,
99        Some("LRO".to_string()),
100    )?;
101
102    println!("{traj_as_flown}");
103
104    // ====================== //
105    // === MODEL MATCHING === //
106    // ====================== //
107
108    // Set up the spacecraft dynamics.
109
110    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
111    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
112    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
113
114    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
115    // We're using the GRAIL JGGRX model.
116    let mut jggrx_meta = MetaFile {
117        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
118        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
119    };
120    // And let's download it if we don't have it yet.
121    jggrx_meta.process(true)?;
122
123    // Build the spherical harmonics.
124    // The harmonics must be computed in the body fixed frame.
125    // We're using the long term prediction of the Moon principal axes frame.
126    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
127    let sph_harmonics = Harmonics::from_stor(
128        almanac.frame_from_uid(moon_pa_frame)?,
129        HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
130    );
131
132    // Include the spherical harmonics into the orbital dynamics.
133    orbital_dyn.accel_models.push(sph_harmonics);
134
135    // We define the solar radiation pressure, using the default solar flux and accounting only
136    // for the eclipsing caused by the Earth and Moon.
137    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
138    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
139
140    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
141    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
142    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
143
144    println!("{dynamics}");
145
146    // Now we can build the propagator.
147    let setup = Propagator::default_dp78(dynamics.clone());
148
149    // For reference, let's build the trajectory with Nyx's models from that LRO state.
150    let (sim_final, traj_as_sim) = setup
151        .with(*traj_as_flown.first(), almanac.clone())
152        .until_epoch_with_traj(traj_as_flown.last().epoch())?;
153
154    println!("SIM INIT:  {:x}", traj_as_flown.first());
155    println!("SIM FINAL: {sim_final:x}");
156    // Compute RIC difference between SIM and LRO ephem
157    let sim_lro_delta = sim_final
158        .orbit
159        .ric_difference(&traj_as_flown.last().orbit)?;
160    println!("{traj_as_sim}");
161    println!(
162        "SIM v LRO - RIC Position (m): {:.3}",
163        sim_lro_delta.radius_km * 1e3
164    );
165    println!(
166        "SIM v LRO - RIC Velocity (m/s): {:.3}",
167        sim_lro_delta.velocity_km_s * 1e3
168    );
169
170    traj_as_sim.ric_diff_to_parquet(
171        &traj_as_flown,
172        "./04_lro_sim_truth_error.parquet",
173        ExportCfg::default(),
174    )?;
175
176    // ==================== //
177    // === OD SIMULATOR === //
178    // ==================== //
179
180    // After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
181    // and the truth LRO state.
182
183    // Therefore, we will actually run an estimation from a dispersed LRO state.
184    // The sc_seed is the true LRO state from the BSP.
185    let sc_seed = *traj_as_flown.first();
186
187    // Load the Deep Space Network ground stations.
188    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
189    let ground_station_file: PathBuf = [
190        env!("CARGO_MANIFEST_DIR"),
191        "examples",
192        "04_lro_od",
193        "dsn-network.yaml",
194    ]
195    .iter()
196    .collect();
197
198    let devices = GroundStation::load_named(ground_station_file)?;
199
200    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
201    // Nyx can build a tracking schedule for you based on the first station with access.
202    let trkconfg_yaml: PathBuf = [
203        env!("CARGO_MANIFEST_DIR"),
204        "examples",
205        "04_lro_od",
206        "tracking-cfg.yaml",
207    ]
208    .iter()
209    .collect();
210
211    let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
212
213    // Build the tracking arc simulation to generate a "standard measurement".
214    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
215        devices.clone(),
216        traj_as_flown.clone(),
217        configs,
218    )?;
219
220    trk.build_schedule(almanac.clone())?;
221    let arc = trk.generate_measurements(almanac.clone())?;
222    // Save the simulated tracking data
223    arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
224
225    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
226    println!("{arc}");
227
228    // Now that we have simulated measurements, we'll run the orbit determination.
229
230    // ===================== //
231    // === OD ESTIMATION === //
232    // ===================== //
233
234    let sc = SpacecraftUncertainty::builder()
235        .nominal(sc_seed)
236        .frame(LocalFrame::RIC)
237        .x_km(0.5)
238        .y_km(0.5)
239        .z_km(0.5)
240        .vx_km_s(5e-3)
241        .vy_km_s(5e-3)
242        .vz_km_s(5e-3)
243        .build();
244
245    // Build the filter initial estimate, which we will reuse in the filter.
246    let initial_estimate = sc.to_estimate()?;
247
248    println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
249
250    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
251    let process_noise = ProcessNoise3D::from_velocity_km_s(
252        &[1.8e-9, 1.8e-9, 1.8e-9],
253        1 * Unit::Hour,
254        10 * Unit::Minute,
255        None,
256    );
257
258    println!("{process_noise}");
259
260    let kf = KF::new(
261        // Increase the initial covariance to account for larger deviation.
262        initial_estimate,
263        process_noise,
264    );
265
266    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
267    let mut odp = SpacecraftODProcess::ekf(
268        setup.with(initial_estimate.state().with_stm(), almanac.clone()),
269        kf,
270        devices,
271        EkfTrigger::new(0, Unit::Hour * 1),
272        Some(ResidRejectCrit::default()),
273        almanac.clone(),
274    );
275
276    odp.process_arc(&arc)?;
277
278    let ric_err = traj_as_flown
279        .at(odp.estimates.last().unwrap().epoch())?
280        .orbit
281        .ric_difference(&odp.estimates.last().unwrap().orbital_state())?;
282    println!("== RIC at end ==");
283    println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
284    println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
285
286    odp.to_parquet(&arc, "./04_lro_od_results.parquet", ExportCfg::default())?;
287
288    // In our case, we have the truth trajectory from NASA.
289    // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
290    // Export the OD trajectory first.
291    let od_trajectory = odp.to_traj()?;
292    // Build the RIC difference.
293    od_trajectory.ric_diff_to_parquet(
294        &traj_as_flown,
295        "./04_lro_od_truth_error.parquet",
296        ExportCfg::default(),
297    )?;
298
299    Ok(())
300}
More examples
Hide additional examples
examples/01_orbit_prop/main.rs (line 231)
30fn main() -> Result<(), Box<dyn Error>> {
31    pel::init();
32    // Dynamics models require planetary constants and ephemerides to be defined.
33    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
34    // This will automatically download the DE440s planetary ephemeris,
35    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
36    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
37    // planetary constants kernels.
38    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
39    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
40    // references to many functions.
41    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
42    // Define the orbit epoch
43    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45    // Define the orbit.
46    // First we need to fetch the Earth J2000 from information from the Almanac.
47    // This allows the frame to include the gravitational parameters and the shape of the Earth,
48    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
49    // by loading a different set of planetary constants.
50    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
51
52    let orbit =
53        Orbit::try_keplerian_altitude(300.0, 0.015, 68.5, 65.2, 75.0, 0.0, epoch, earth_j2000)?;
54    // Print in in Keplerian form.
55    println!("{orbit:x}");
56
57    // There are two ways to propagate an orbit. We can make a quick approximation assuming only two-body
58    // motion. This is a useful first order approximation but it isn't used in real-world applications.
59
60    // This approach is a feature of ANISE.
61    let future_orbit_tb = orbit.at_epoch(epoch + Unit::Day * 3)?;
62    println!("{future_orbit_tb:x}");
63
64    // Two body propagation relies solely on Kepler's laws, so only the true anomaly will change.
65    println!(
66        "SMA changed by {:.3e} km",
67        orbit.sma_km()? - future_orbit_tb.sma_km()?
68    );
69    println!(
70        "ECC changed by {:.3e}",
71        orbit.ecc()? - future_orbit_tb.ecc()?
72    );
73    println!(
74        "INC changed by {:.3e} deg",
75        orbit.inc_deg()? - future_orbit_tb.inc_deg()?
76    );
77    println!(
78        "RAAN changed by {:.3e} deg",
79        orbit.raan_deg()? - future_orbit_tb.raan_deg()?
80    );
81    println!(
82        "AOP changed by {:.3e} deg",
83        orbit.aop_deg()? - future_orbit_tb.aop_deg()?
84    );
85    println!(
86        "TA changed by {:.3} deg",
87        orbit.ta_deg()? - future_orbit_tb.ta_deg()?
88    );
89
90    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
91    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
92    // models such as solar radiation pressure.
93
94    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
95    let sc = Spacecraft::builder()
96        .orbit(orbit)
97        .mass(Mass::from_dry_mass(9.60))
98        .srp(SRPData {
99            area_m2: 10e-4,
100            coeff_reflectivity: 1.1,
101        })
102        .build();
103    println!("{sc:x}");
104
105    // Set up the spacecraft dynamics.
106
107    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
108    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
109    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
110
111    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
112    // We're using the JGM3 model here, which is the default in GMAT.
113    let mut jgm3_meta = MetaFile {
114        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
115        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
116    };
117    // And let's download it if we don't have it yet.
118    jgm3_meta.process(true)?;
119
120    // Build the spherical harmonics.
121    // The harmonics must be computed in the body fixed frame.
122    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
123    let harmonics_21x21 = Harmonics::from_stor(
124        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
125        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
126    );
127
128    // Include the spherical harmonics into the orbital dynamics.
129    orbital_dyn.accel_models.push(harmonics_21x21);
130
131    // We define the solar radiation pressure, using the default solar flux and accounting only
132    // for the eclipsing caused by the Earth.
133    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
134
135    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
136    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
137    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
138
139    println!("{dynamics}");
140
141    // Finally, let's propagate this orbit to the same epoch as above.
142    // The first returned value is the spacecraft state at the final epoch.
143    // The second value is the full trajectory where the step size is variable step used by the propagator.
144    let (future_sc, trajectory) = Propagator::default(dynamics)
145        .with(sc, almanac.clone())
146        .until_epoch_with_traj(future_orbit_tb.epoch)?;
147
148    println!("=== High fidelity propagation ===");
149    println!(
150        "SMA changed by {:.3} km",
151        orbit.sma_km()? - future_sc.orbit.sma_km()?
152    );
153    println!(
154        "ECC changed by {:.6}",
155        orbit.ecc()? - future_sc.orbit.ecc()?
156    );
157    println!(
158        "INC changed by {:.3e} deg",
159        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
160    );
161    println!(
162        "RAAN changed by {:.3} deg",
163        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
164    );
165    println!(
166        "AOP changed by {:.3} deg",
167        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
168    );
169    println!(
170        "TA changed by {:.3} deg",
171        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
172    );
173
174    // We also have access to the full trajectory throughout the propagation.
175    println!("{trajectory}");
176
177    // With the trajectory, let's build a few data products.
178
179    // 1. Export the trajectory as a CCSDS OEM version 2.0 file and as a parquet file, which includes the Keplerian orbital elements.
180
181    trajectory.to_oem_file(
182        "./01_cubesat_hf_prop.oem",
183        ExportCfg::builder().step(Unit::Minute * 2).build(),
184    )?;
185
186    trajectory.to_parquet_with_cfg(
187        "./01_cubesat_hf_prop.parquet",
188        ExportCfg::builder().step(Unit::Minute * 2).build(),
189        almanac.clone(),
190    )?;
191
192    // 2. Compare the difference in the radial-intrack-crosstrack frame between the high fidelity
193    // and Keplerian propagation. The RIC frame is commonly used to compute the difference in position
194    // and velocity of different spacecraft.
195    // 3. Compute the azimuth, elevation, range, and range-rate data of that spacecraft as seen from Boulder, CO, USA.
196
197    let boulder_station = GroundStation::from_point(
198        "Boulder, CO, USA".to_string(),
199        40.014984,   // latitude in degrees
200        -105.270546, // longitude in degrees
201        1.6550,      // altitude in kilometers
202        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
203    );
204
205    // We iterate over the trajectory, grabbing a state every two minutes.
206    let mut offset_s = vec![];
207    let mut epoch_str = vec![];
208    let mut ric_x_km = vec![];
209    let mut ric_y_km = vec![];
210    let mut ric_z_km = vec![];
211    let mut ric_vx_km_s = vec![];
212    let mut ric_vy_km_s = vec![];
213    let mut ric_vz_km_s = vec![];
214
215    let mut azimuth_deg = vec![];
216    let mut elevation_deg = vec![];
217    let mut range_km = vec![];
218    let mut range_rate_km_s = vec![];
219    for state in trajectory.every(Unit::Minute * 2) {
220        // Try to compute the Keplerian/two body state just in time.
221        // This method occasionally fails to converge on an appropriate true anomaly
222        // from the mean anomaly. If that happens, we just skip this state.
223        // The high fidelity and Keplerian states diverge continuously, and we're curious
224        // about the divergence in this quick analysis.
225        let this_epoch = state.epoch();
226        match orbit.at_epoch(this_epoch) {
227            Ok(tb_then) => {
228                offset_s.push((this_epoch - orbit.epoch).to_seconds());
229                epoch_str.push(format!("{this_epoch}"));
230                // Compute the two body state just in time.
231                let ric = state.orbit.ric_difference(&tb_then)?;
232                ric_x_km.push(ric.radius_km.x);
233                ric_y_km.push(ric.radius_km.y);
234                ric_z_km.push(ric.radius_km.z);
235                ric_vx_km_s.push(ric.velocity_km_s.x);
236                ric_vy_km_s.push(ric.velocity_km_s.y);
237                ric_vz_km_s.push(ric.velocity_km_s.z);
238
239                // Compute the AER data for each state.
240                let aer = almanac.azimuth_elevation_range_sez(
241                    state.orbit,
242                    boulder_station.to_orbit(this_epoch, &almanac)?,
243                    None,
244                    None,
245                )?;
246                azimuth_deg.push(aer.azimuth_deg);
247                elevation_deg.push(aer.elevation_deg);
248                range_km.push(aer.range_km);
249                range_rate_km_s.push(aer.range_rate_km_s);
250            }
251            Err(e) => warn!("{} {e}", state.epoch()),
252        };
253    }
254
255    // Build the data frames.
256    let ric_df = df!(
257        "Offset (s)" => offset_s.clone(),
258        "Epoch" => epoch_str.clone(),
259        "RIC X (km)" => ric_x_km,
260        "RIC Y (km)" => ric_y_km,
261        "RIC Z (km)" => ric_z_km,
262        "RIC VX (km/s)" => ric_vx_km_s,
263        "RIC VY (km/s)" => ric_vy_km_s,
264        "RIC VZ (km/s)" => ric_vz_km_s,
265    )?;
266
267    println!("RIC difference at start\n{}", ric_df.head(Some(10)));
268    println!("RIC difference at end\n{}", ric_df.tail(Some(10)));
269
270    let aer_df = df!(
271        "Offset (s)" => offset_s.clone(),
272        "Epoch" => epoch_str.clone(),
273        "azimuth (deg)" => azimuth_deg,
274        "elevation (deg)" => elevation_deg,
275        "range (km)" => range_km,
276        "range rate (km/s)" => range_rate_km_s,
277    )?;
278
279    // Finally, let's see when the spacecraft is visible, assuming 15 degrees minimum elevation.
280    let mask = aer_df
281        .column("elevation (deg)")?
282        .gt(&Column::Scalar(ScalarColumn::new(
283            "elevation mask (deg)".into(),
284            Scalar::new(DataType::Float64, AnyValue::Float64(15.0)),
285            offset_s.len(),
286        )))?;
287    let cubesat_visible = aer_df.filter(&mask)?;
288
289    println!("{cubesat_visible}");
290
291    Ok(())
292}

pub fn vnc_difference( &self, other: &CartesianState, ) -> Result<CartesianState, PhysicsError>

Returns a Cartesian state representing the VNC difference between self and other, in position and velocity (with transport theorem). Refer to dcm_from_vnc_to_inertial for details on the VNC frame.

§Algorithm
  1. Compute the VNC DCM of self
  2. Rotate self into the VNC frame
  3. Rotation other into the VNC frame
  4. Compute the difference between these two states
  5. Strip the astrodynamical information from the frame, enabling only computations from CartesianState

:type other: Orbit :rtype: Orbit

§

impl CartesianState

pub fn rmag_km(&self) -> f64

Returns the magnitude of the radius vector in km

:rtype: float

pub fn vmag_km_s(&self) -> f64

Returns the magnitude of the velocity vector in km/s

:rtype: float

pub fn distance_to_km( &self, other: &CartesianState, ) -> Result<f64, PhysicsError>

Returns the distance in kilometers between this state and another state, if both frame match (epoch does not need to match).

:type other: Orbit :rtype: float

pub fn rss_radius_km(&self, other: &CartesianState) -> Result<f64, PhysicsError>

Returns the root mean squared (RSS) radius difference between this state and another state, if both frames match (epoch does not need to match)

:type other: Orbit :rtype: float

pub fn rss_velocity_km_s( &self, other: &CartesianState, ) -> Result<f64, PhysicsError>

Returns the root mean squared (RSS) velocity difference between this state and another state, if both frames match (epoch does not need to match)

:type other: Orbit :rtype: float

pub fn rms_radius_km(&self, other: &CartesianState) -> Result<f64, PhysicsError>

Returns the root sum squared (RMS) radius difference between this state and another state, if both frames match (epoch does not need to match)

:type other: Orbit :rtype: float

pub fn rms_velocity_km_s( &self, other: &CartesianState, ) -> Result<f64, PhysicsError>

Returns the root sum squared (RMS) velocity difference between this state and another state, if both frames match (epoch does not need to match)

:type other: Orbit :rtype: float

pub fn eq_within( &self, other: &CartesianState, radial_tol_km: f64, velocity_tol_km_s: f64, ) -> bool

Returns whether this orbit and another are equal within the specified radial and velocity absolute tolerances

:type other: Orbit :type radial_tol_km: float :type velocity_tol_km_s: float :rtype: bool

pub fn light_time(&self) -> Duration

Returns the light time duration between this object and the origin of its reference frame.

:rtype: Duration

pub fn abs_pos_diff_km( &self, other: &CartesianState, ) -> Result<f64, PhysicsError>

Returns the absolute position difference in kilometer between this orbit and another. Raises an error if the frames do not match (epochs do not need to match).

:type other: Orbit :rtype: float

pub fn abs_vel_diff_km_s( &self, other: &CartesianState, ) -> Result<f64, PhysicsError>

Returns the absolute velocity difference in kilometer per second between this orbit and another. Raises an error if the frames do not match (epochs do not need to match).

:type other: Orbit :rtype: float

pub fn abs_difference( &self, other: &CartesianState, ) -> Result<(f64, f64), PhysicsError>

Returns the absolute position and velocity differences in km and km/s between this orbit and another. Raises an error if the frames do not match (epochs do not need to match).

:type other: Orbit :rtype: typing.Tuple

pub fn rel_pos_diff(&self, other: &CartesianState) -> Result<f64, PhysicsError>

Returns the relative position difference (unitless) between this orbit and another. This is computed by dividing the absolute difference by the norm of this object’s radius vector. If the radius is zero, this function raises a math error. Raises an error if the frames do not match or (epochs do not need to match).

:type other: Orbit :rtype: float

pub fn rel_vel_diff(&self, other: &CartesianState) -> Result<f64, PhysicsError>

Returns the absolute velocity difference in kilometer per second between this orbit and another. Raises an error if the frames do not match (epochs do not need to match).

:type other: Orbit :rtype: float

pub fn rel_difference( &self, other: &CartesianState, ) -> Result<(f64, f64), PhysicsError>

Returns the relative difference between this orbit and another for the position and velocity, respectively the first and second return values. Both return values are UNITLESS because the relative difference is computed as the absolute difference divided by the rmag and vmag of this object. Raises an error if the frames do not match, if the position is zero or the velocity is zero.

:type other: Orbit :rtype: typing.Tuple

§

impl CartesianState

pub fn sma_altitude_km(&self) -> Result<f64, PhysicsError>

Returns the SMA altitude in km

:rtype: float

pub fn periapsis_altitude_km(&self) -> Result<f64, PhysicsError>

Returns the altitude of periapsis (or perigee around Earth), in kilometers.

:rtype: float

pub fn apoapsis_altitude_km(&self) -> Result<f64, PhysicsError>

Returns the altitude of apoapsis (or apogee around Earth), in kilometers.

:rtype: float

pub fn latlongalt(&self) -> Result<(f64, f64, f64), PhysicsError>

Returns the geodetic latitude, geodetic longitude, and geodetic height, respectively in degrees, degrees, and kilometers.

§Algorithm

This uses the Heikkinen procedure, which is not iterative. The results match Vallado and GMAT.

:rtype: typing.Tuple

Examples found in repository?
examples/03_geo_analysis/drift.rs (line 55)
26fn main() -> Result<(), Box<dyn Error>> {
27    pel::init();
28    // Dynamics models require planetary constants and ephemerides to be defined.
29    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30    // This will automatically download the DE440s planetary ephemeris,
31    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
32    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
33    // planetary constants kernels.
34    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
35    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
36    // references to many functions.
37    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
38    // Define the orbit epoch
39    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
40
41    // Define the orbit.
42    // First we need to fetch the Earth J2000 from information from the Almanac.
43    // This allows the frame to include the gravitational parameters and the shape of the Earth,
44    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
45    // by loading a different set of planetary constants.
46    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
47
48    // Placing this GEO bird just above Colorado.
49    // In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
50    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
51    // Print in in Keplerian form.
52    println!("{orbit:x}");
53
54    let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
55    let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;
56
57    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
58    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
59    // models such as solar radiation pressure.
60
61    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
62    let sc = Spacecraft::builder()
63        .orbit(orbit)
64        .mass(Mass::from_dry_mass(9.60))
65        .srp(SRPData {
66            area_m2: 10e-4,
67            coeff_reflectivity: 1.1,
68        })
69        .build();
70    println!("{sc:x}");
71
72    // Set up the spacecraft dynamics.
73
74    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
75    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
76    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
77
78    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
79    // We're using the JGM3 model here, which is the default in GMAT.
80    let mut jgm3_meta = MetaFile {
81        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
82        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
83    };
84    // And let's download it if we don't have it yet.
85    jgm3_meta.process(true)?;
86
87    // Build the spherical harmonics.
88    // The harmonics must be computed in the body fixed frame.
89    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
90    let harmonics_21x21 = Harmonics::from_stor(
91        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
92        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
93    );
94
95    // Include the spherical harmonics into the orbital dynamics.
96    orbital_dyn.accel_models.push(harmonics_21x21);
97
98    // We define the solar radiation pressure, using the default solar flux and accounting only
99    // for the eclipsing caused by the Earth and Moon.
100    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
103    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
104    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
105
106    println!("{dynamics}");
107
108    // Finally, let's propagate this orbit to the same epoch as above.
109    // The first returned value is the spacecraft state at the final epoch.
110    // The second value is the full trajectory where the step size is variable step used by the propagator.
111    let (future_sc, trajectory) = Propagator::default(dynamics)
112        .with(sc, almanac.clone())
113        .until_epoch_with_traj(epoch + Unit::Century * 0.03)?;
114
115    println!("=== High fidelity propagation ===");
116    println!(
117        "SMA changed by {:.3} km",
118        orbit.sma_km()? - future_sc.orbit.sma_km()?
119    );
120    println!(
121        "ECC changed by {:.6}",
122        orbit.ecc()? - future_sc.orbit.ecc()?
123    );
124    println!(
125        "INC changed by {:.3e} deg",
126        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
127    );
128    println!(
129        "RAAN changed by {:.3} deg",
130        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
131    );
132    println!(
133        "AOP changed by {:.3} deg",
134        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
135    );
136    println!(
137        "TA changed by {:.3} deg",
138        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
139    );
140
141    // We also have access to the full trajectory throughout the propagation.
142    println!("{trajectory}");
143
144    println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");
145
146    // With the trajectory, let's build a few data products.
147
148    // 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.
149
150    let analysis_step = Unit::Minute * 5;
151
152    trajectory.to_parquet(
153        "./03_geo_hf_prop.parquet",
154        Some(vec![
155            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
156        ]),
157        ExportCfg::builder().step(analysis_step).build(),
158        almanac.clone(),
159    )?;
160
161    // 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.
162
163    // We iterate over the trajectory, grabbing a state every two minutes.
164    let mut offset_s = vec![];
165    let mut epoch_str = vec![];
166    let mut longitude_deg = vec![];
167    let mut latitude_deg = vec![];
168    let mut altitude_km = vec![];
169
170    for state in trajectory.every(analysis_step) {
171        // Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
172        // These define the GEO stationkeeping box.
173
174        let this_epoch = state.epoch();
175
176        offset_s.push((this_epoch - orbit.epoch).to_seconds());
177        epoch_str.push(this_epoch.to_isoformat());
178
179        let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
180        let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
181        longitude_deg.push(long_deg);
182        latitude_deg.push(lat_deg);
183        altitude_km.push(alt_km);
184    }
185
186    println!(
187        "Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
188        orig_long_deg - longitude_deg.last().unwrap()
189    );
190
191    println!(
192        "Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
193        orig_lat_deg - latitude_deg.last().unwrap()
194    );
195
196    println!(
197        "Altitude changed by {:.3} km -- Box is 30 km",
198        orig_alt_km - altitude_km.last().unwrap()
199    );
200
201    // Build the station keeping data frame.
202    let mut sk_df = df!(
203        "Offset (s)" => offset_s.clone(),
204        "Epoch (UTC)" => epoch_str.clone(),
205        "Longitude E-W (deg)" => longitude_deg,
206        "Latitude N-S (deg)" => latitude_deg,
207        "Altitude (km)" => altitude_km,
208
209    )?;
210
211    // Create a file to write the Parquet to
212    let file = File::create("./03_geo_lla.parquet").expect("Could not create file");
213
214    // Create a ParquetWriter and write the DataFrame to the file
215    ParquetWriter::new(file).finish(&mut sk_df)?;
216
217    Ok(())
218}

pub fn longitude_deg(&self) -> f64

Returns the geodetic longitude (λ) in degrees. Value is between -180 and 180 degrees.

§Frame warning

This state MUST be in the body fixed frame (e.g. ITRF93) prior to calling this function, or the computation is invalid.

:rtype: float

pub fn longitude_360_deg(&self) -> f64

Returns the geodetic longitude (λ) in degrees. Value is between 0 and 360 degrees.

§Frame warning

This state MUST be in the body fixed frame (e.g. ITRF93) prior to calling this function, or the computation is invalid.

:rtype: float

pub fn latitude_deg(&self) -> Result<f64, PhysicsError>

Returns the geodetic latitude (φ) in degrees. Value is between -180 and +180 degrees.

§Frame warning

This state MUST be in the body fixed frame (e.g. ITRF93) prior to calling this function, or the computation is invalid.

:rtype: float

pub fn height_km(&self) -> Result<f64, PhysicsError>

Returns the geodetic height in km.

Reference: Vallado, 4th Ed., Algorithm 12 page 172.

:rtype: float

§

impl CartesianState

pub fn try_keplerian( sma_km: f64, ecc: f64, inc_deg: f64, raan_deg: f64, aop_deg: f64, ta_deg: f64, epoch: Epoch, frame: Frame, ) -> Result<CartesianState, PhysicsError>

Attempts to create a new Orbit around the provided Celestial or Geoid frame from the Keplerian orbital elements.

Units: km, none, degrees, degrees, degrees, degrees

WARNING: This function will return an error if the singularities in the conversion are encountered. NOTE: The state is defined in Cartesian coordinates as they are non-singular. This causes rounding errors when creating a state from its Keplerian orbital elements (cf. the state tests). One should expect these errors to be on the order of 1e-12.

Examples found in repository?
examples/03_geo_analysis/stationkeeping.rs (line 36)
28fn main() -> Result<(), Box<dyn Error>> {
29    pel::init();
30    // Set up the dynamics like in the orbit raise.
31    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
32    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
33
34    // Define the GEO orbit, and we're just going to maintain it very tightly.
35    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
36    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
37    println!("{orbit:x}");
38
39    let sc = Spacecraft::builder()
40        .orbit(orbit)
41        .mass(Mass::from_dry_and_prop_masses(1000.0, 1000.0)) // 1000 kg of dry mass and prop, totalling 2.0 tons
42        .srp(SRPData::from_area(3.0 * 6.0)) // Assuming 1 kW/m^2 or 18 kW, giving a margin of 4.35 kW for on-propulsion consumption
43        .thruster(Thruster {
44            // "NEXT-STEP" row in Table 2
45            isp_s: 4435.0,
46            thrust_N: 0.472,
47        })
48        .mode(GuidanceMode::Thrust) // Start thrusting immediately.
49        .build();
50
51    // Set up the spacecraft dynamics like in the orbit raise example.
52
53    let prop_time = 30.0 * Unit::Day;
54
55    // Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
56    let objectives = &[
57        Objective::within_tolerance(StateParameter::SMA, 42_164.0, 5.0), // 5 km
58        Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5),
59        Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2),
60    ];
61
62    let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2)?;
63    println!("{ruggiero_ctrl}");
64
65    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
66
67    let mut jgm3_meta = MetaFile {
68        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
69        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
70    };
71    jgm3_meta.process(true)?;
72
73    let harmonics = Harmonics::from_stor(
74        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
75        HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true)?,
76    );
77    orbital_dyn.accel_models.push(harmonics);
78
79    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
80    let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
81        .with_guidance_law(ruggiero_ctrl.clone());
82
83    println!("{sc_dynamics}");
84
85    // Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.
86
87    // Let's start by defining the dispersion.
88    // The MultivariateNormal structure allows us to define the dispersions in any of the orbital parameters, but these are applied directly in the Cartesian state space.
89    // Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
90    let mc_rv = MvnSpacecraft::new(
91        sc,
92        vec![StateDispersion::zero_mean(StateParameter::SMA, 3.0)],
93    )?;
94
95    let my_mc = MonteCarlo::new(
96        sc, // Nominal state
97        mc_rv,
98        "03_geo_sk".to_string(), // Scenario name
99        None, // No specific seed specified, so one will be drawn from the computer's entropy.
100    );
101
102    // Build the propagator setup.
103    let setup = Propagator::rk89(
104        sc_dynamics.clone(),
105        IntegratorOptions::builder()
106            .min_step(10.0_f64.seconds())
107            .error_ctrl(ErrorControl::RSSCartesianStep)
108            .build(),
109    );
110
111    let num_runs = 25;
112    let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);
113
114    assert_eq!(rslts.runs.len(), num_runs);
115
116    // For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
117
118    rslts.to_parquet(
119        "03_geo_sk.parquet",
120        Some(vec![
121            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
122        ]),
123        ExportCfg::default(),
124        almanac,
125    )?;
126
127    Ok(())
128}
More examples
Hide additional examples
examples/03_geo_analysis/drift.rs (line 50)
26fn main() -> Result<(), Box<dyn Error>> {
27    pel::init();
28    // Dynamics models require planetary constants and ephemerides to be defined.
29    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30    // This will automatically download the DE440s planetary ephemeris,
31    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
32    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
33    // planetary constants kernels.
34    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
35    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
36    // references to many functions.
37    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
38    // Define the orbit epoch
39    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
40
41    // Define the orbit.
42    // First we need to fetch the Earth J2000 from information from the Almanac.
43    // This allows the frame to include the gravitational parameters and the shape of the Earth,
44    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
45    // by loading a different set of planetary constants.
46    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
47
48    // Placing this GEO bird just above Colorado.
49    // In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
50    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
51    // Print in in Keplerian form.
52    println!("{orbit:x}");
53
54    let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
55    let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;
56
57    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
58    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
59    // models such as solar radiation pressure.
60
61    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
62    let sc = Spacecraft::builder()
63        .orbit(orbit)
64        .mass(Mass::from_dry_mass(9.60))
65        .srp(SRPData {
66            area_m2: 10e-4,
67            coeff_reflectivity: 1.1,
68        })
69        .build();
70    println!("{sc:x}");
71
72    // Set up the spacecraft dynamics.
73
74    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
75    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
76    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
77
78    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
79    // We're using the JGM3 model here, which is the default in GMAT.
80    let mut jgm3_meta = MetaFile {
81        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
82        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
83    };
84    // And let's download it if we don't have it yet.
85    jgm3_meta.process(true)?;
86
87    // Build the spherical harmonics.
88    // The harmonics must be computed in the body fixed frame.
89    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
90    let harmonics_21x21 = Harmonics::from_stor(
91        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
92        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
93    );
94
95    // Include the spherical harmonics into the orbital dynamics.
96    orbital_dyn.accel_models.push(harmonics_21x21);
97
98    // We define the solar radiation pressure, using the default solar flux and accounting only
99    // for the eclipsing caused by the Earth and Moon.
100    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
103    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
104    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
105
106    println!("{dynamics}");
107
108    // Finally, let's propagate this orbit to the same epoch as above.
109    // The first returned value is the spacecraft state at the final epoch.
110    // The second value is the full trajectory where the step size is variable step used by the propagator.
111    let (future_sc, trajectory) = Propagator::default(dynamics)
112        .with(sc, almanac.clone())
113        .until_epoch_with_traj(epoch + Unit::Century * 0.03)?;
114
115    println!("=== High fidelity propagation ===");
116    println!(
117        "SMA changed by {:.3} km",
118        orbit.sma_km()? - future_sc.orbit.sma_km()?
119    );
120    println!(
121        "ECC changed by {:.6}",
122        orbit.ecc()? - future_sc.orbit.ecc()?
123    );
124    println!(
125        "INC changed by {:.3e} deg",
126        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
127    );
128    println!(
129        "RAAN changed by {:.3} deg",
130        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
131    );
132    println!(
133        "AOP changed by {:.3} deg",
134        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
135    );
136    println!(
137        "TA changed by {:.3} deg",
138        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
139    );
140
141    // We also have access to the full trajectory throughout the propagation.
142    println!("{trajectory}");
143
144    println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");
145
146    // With the trajectory, let's build a few data products.
147
148    // 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.
149
150    let analysis_step = Unit::Minute * 5;
151
152    trajectory.to_parquet(
153        "./03_geo_hf_prop.parquet",
154        Some(vec![
155            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
156        ]),
157        ExportCfg::builder().step(analysis_step).build(),
158        almanac.clone(),
159    )?;
160
161    // 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.
162
163    // We iterate over the trajectory, grabbing a state every two minutes.
164    let mut offset_s = vec![];
165    let mut epoch_str = vec![];
166    let mut longitude_deg = vec![];
167    let mut latitude_deg = vec![];
168    let mut altitude_km = vec![];
169
170    for state in trajectory.every(analysis_step) {
171        // Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
172        // These define the GEO stationkeeping box.
173
174        let this_epoch = state.epoch();
175
176        offset_s.push((this_epoch - orbit.epoch).to_seconds());
177        epoch_str.push(this_epoch.to_isoformat());
178
179        let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
180        let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
181        longitude_deg.push(long_deg);
182        latitude_deg.push(lat_deg);
183        altitude_km.push(alt_km);
184    }
185
186    println!(
187        "Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
188        orig_long_deg - longitude_deg.last().unwrap()
189    );
190
191    println!(
192        "Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
193        orig_lat_deg - latitude_deg.last().unwrap()
194    );
195
196    println!(
197        "Altitude changed by {:.3} km -- Box is 30 km",
198        orig_alt_km - altitude_km.last().unwrap()
199    );
200
201    // Build the station keeping data frame.
202    let mut sk_df = df!(
203        "Offset (s)" => offset_s.clone(),
204        "Epoch (UTC)" => epoch_str.clone(),
205        "Longitude E-W (deg)" => longitude_deg,
206        "Latitude N-S (deg)" => latitude_deg,
207        "Altitude (km)" => altitude_km,
208
209    )?;
210
211    // Create a file to write the Parquet to
212    let file = File::create("./03_geo_lla.parquet").expect("Could not create file");
213
214    // Create a ParquetWriter and write the DataFrame to the file
215    ParquetWriter::new(file).finish(&mut sk_df)?;
216
217    Ok(())
218}

pub fn try_keplerian_apsis_radii( r_a_km: f64, r_p_km: f64, inc_deg: f64, raan_deg: f64, aop_deg: f64, ta_deg: f64, epoch: Epoch, frame: Frame, ) -> Result<CartesianState, PhysicsError>

Attempts to create a new Orbit from the provided radii of apoapsis and periapsis, in kilometers

pub fn try_keplerian_vec( state: &Matrix<f64, Const<6>, Const<1>, ArrayStorage<f64, 6, 1>>, epoch: Epoch, frame: Frame, ) -> Result<CartesianState, PhysicsError>

Attempts to create a new Orbit around the provided frame from the borrowed state vector

The state vector must be sma, ecc, inc, raan, aop, ta. This function is a shortcut to cartesian and as such it has the same unit requirements.

pub fn keplerian( sma_km: f64, ecc: f64, inc_deg: f64, raan_deg: f64, aop_deg: f64, ta_deg: f64, epoch: Epoch, frame: Frame, ) -> CartesianState

Creates (without error checking) a new Orbit around the provided Celestial or Geoid frame from the Keplerian orbital elements.

Units: km, none, degrees, degrees, degrees, degrees

WARNING: This function will panic if the singularities in the conversion are expected. NOTE: The state is defined in Cartesian coordinates as they are non-singular. This causes rounding errors when creating a state from its Keplerian orbital elements (cf. the state tests). One should expect these errors to be on the order of 1e-12.

Examples found in repository?
examples/03_geo_analysis/raise.rs (line 50)
27fn main() -> Result<(), Box<dyn Error>> {
28    pel::init();
29
30    // Dynamics models require planetary constants and ephemerides to be defined.
31    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
32    // This will automatically download the DE440s planetary ephemeris,
33    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
34    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
35    // planetary constants kernels.
36    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
37    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
38    // references to many functions.
39    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
40    // Fetch the EME2000 frame from the Almabac
41    let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
42    // Define the orbit epoch
43    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45    // Build the spacecraft itself.
46    // Using slide 6 of https://aerospace.org/sites/default/files/2018-11/Davis-Mayberry_HPSEP_11212018.pdf
47    // for the "next gen" SEP characteristics.
48
49    // GTO start
50    let orbit = Orbit::keplerian(24505.9, 0.725, 7.05, 0.0, 0.0, 0.0, epoch, eme2k);
51
52    let sc = Spacecraft::builder()
53        .orbit(orbit)
54        .mass(Mass::from_dry_and_prop_masses(1000.0, 1000.0)) // 1000 kg of dry mass and prop, totalling 2.0 tons
55        .srp(SRPData::from_area(3.0 * 6.0)) // Assuming 1 kW/m^2 or 18 kW, giving a margin of 4.35 kW for on-propulsion consumption
56        .thruster(Thruster {
57            // "NEXT-STEP" row in Table 2
58            isp_s: 4435.0,
59            thrust_N: 0.472,
60        })
61        .mode(GuidanceMode::Thrust) // Start thrusting immediately.
62        .build();
63
64    let prop_time = 180.0 * Unit::Day;
65
66    // Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
67    let objectives = &[
68        Objective::within_tolerance(StateParameter::SMA, 42_165.0, 20.0),
69        Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5),
70        Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2),
71    ];
72
73    // Ensure that we only thrust if we have more than 20% illumination.
74    let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2).unwrap();
75    println!("{ruggiero_ctrl}");
76
77    // Define the high fidelity dynamics
78
79    // Set up the spacecraft dynamics.
80
81    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
82    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
83    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
84
85    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
86    // We're using the JGM3 model here, which is the default in GMAT.
87    let mut jgm3_meta = MetaFile {
88        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
89        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
90    };
91    // And let's download it if we don't have it yet.
92    jgm3_meta.process(true)?;
93
94    // Build the spherical harmonics.
95    // The harmonics must be computed in the body fixed frame.
96    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
97    let harmonics = Harmonics::from_stor(
98        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
99        HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true).unwrap(),
100    );
101
102    // Include the spherical harmonics into the orbital dynamics.
103    orbital_dyn.accel_models.push(harmonics);
104
105    // We define the solar radiation pressure, using the default solar flux and accounting only
106    // for the eclipsing caused by the Earth.
107    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
108
109    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
110    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
111    let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
112        .with_guidance_law(ruggiero_ctrl.clone());
113
114    println!("{:x}", orbit);
115
116    // We specify a minimum step in the propagator because the Ruggiero control would otherwise drive this step very low.
117    let (final_state, traj) = Propagator::rk89(
118        sc_dynamics.clone(),
119        IntegratorOptions::builder()
120            .min_step(10.0_f64.seconds())
121            .error_ctrl(ErrorControl::RSSCartesianStep)
122            .build(),
123    )
124    .with(sc, almanac.clone())
125    .for_duration_with_traj(prop_time)?;
126
127    let prop_usage = sc.mass.prop_mass_kg - final_state.mass.prop_mass_kg;
128    println!("{:x}", final_state.orbit);
129    println!("prop usage: {:.3} kg", prop_usage);
130
131    // Finally, export the results for analysis, including the penumbra percentage throughout the orbit raise.
132    traj.to_parquet(
133        "./03_geo_raise.parquet",
134        Some(vec![
135            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
136        ]),
137        ExportCfg::default(),
138        almanac,
139    )?;
140
141    for status_line in ruggiero_ctrl.status(&final_state) {
142        println!("{status_line}");
143    }
144
145    ruggiero_ctrl
146        .achieved(&final_state)
147        .expect("objective not achieved");
148
149    Ok(())
150}

pub fn keplerian_apsis_radii( r_a_km: f64, r_p_km: f64, inc_deg: f64, raan_deg: f64, aop_deg: f64, ta_deg: f64, epoch: Epoch, frame: Frame, ) -> CartesianState

Creates a new Orbit from the provided radii of apoapsis and periapsis, in kilometers

pub fn try_keplerian_mean_anomaly( sma_km: f64, ecc: f64, inc_deg: f64, raan_deg: f64, aop_deg: f64, ma_deg: f64, epoch: Epoch, frame: Frame, ) -> Result<CartesianState, PhysicsError>

Initializes a new orbit from the Keplerian orbital elements using the mean anomaly instead of the true anomaly.

§Implementation notes

This function starts by converting the mean anomaly to true anomaly, and then it initializes the orbit using the keplerian(..) method. The conversion is from GMAT’s MeanToTrueAnomaly function, transliterated originally by Claude and GPT4 with human adjustments.

pub fn keplerian_vec( state: &Matrix<f64, Const<6>, Const<1>, ArrayStorage<f64, 6, 1>>, epoch: Epoch, frame: Frame, ) -> CartesianState

Creates a new Orbit around the provided frame from the borrowed state vector

The state vector must be sma, ecc, inc, raan, aop, ta. This function is a shortcut to cartesian and as such it has the same unit requirements.

pub fn to_keplerian_vec( self, ) -> Result<Matrix<f64, Const<6>, Const<1>, ArrayStorage<f64, 6, 1>>, PhysicsError>

Returns this state as a Keplerian Vector6 in [km, none, degrees, degrees, degrees, degrees]

Note that the time is not returned in the vector.

pub fn hvec( &self, ) -> Result<Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>, PhysicsError>

Returns the orbital momentum vector

pub fn evec( &self, ) -> Result<Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>, PhysicsError>

Returns the eccentricity vector (no unit)

§

impl CartesianState

pub fn try_keplerian_altitude( sma_altitude_km: f64, ecc: f64, inc_deg: f64, raan_deg: f64, aop_deg: f64, ta_deg: f64, epoch: Epoch, frame: Frame, ) -> Result<CartesianState, PhysicsError>

Creates a new Orbit from the provided semi-major axis altitude in kilometers

Examples found in repository?
examples/01_orbit_prop/main.rs (line 53)
30fn main() -> Result<(), Box<dyn Error>> {
31    pel::init();
32    // Dynamics models require planetary constants and ephemerides to be defined.
33    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
34    // This will automatically download the DE440s planetary ephemeris,
35    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
36    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
37    // planetary constants kernels.
38    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
39    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
40    // references to many functions.
41    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
42    // Define the orbit epoch
43    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45    // Define the orbit.
46    // First we need to fetch the Earth J2000 from information from the Almanac.
47    // This allows the frame to include the gravitational parameters and the shape of the Earth,
48    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
49    // by loading a different set of planetary constants.
50    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
51
52    let orbit =
53        Orbit::try_keplerian_altitude(300.0, 0.015, 68.5, 65.2, 75.0, 0.0, epoch, earth_j2000)?;
54    // Print in in Keplerian form.
55    println!("{orbit:x}");
56
57    // There are two ways to propagate an orbit. We can make a quick approximation assuming only two-body
58    // motion. This is a useful first order approximation but it isn't used in real-world applications.
59
60    // This approach is a feature of ANISE.
61    let future_orbit_tb = orbit.at_epoch(epoch + Unit::Day * 3)?;
62    println!("{future_orbit_tb:x}");
63
64    // Two body propagation relies solely on Kepler's laws, so only the true anomaly will change.
65    println!(
66        "SMA changed by {:.3e} km",
67        orbit.sma_km()? - future_orbit_tb.sma_km()?
68    );
69    println!(
70        "ECC changed by {:.3e}",
71        orbit.ecc()? - future_orbit_tb.ecc()?
72    );
73    println!(
74        "INC changed by {:.3e} deg",
75        orbit.inc_deg()? - future_orbit_tb.inc_deg()?
76    );
77    println!(
78        "RAAN changed by {:.3e} deg",
79        orbit.raan_deg()? - future_orbit_tb.raan_deg()?
80    );
81    println!(
82        "AOP changed by {:.3e} deg",
83        orbit.aop_deg()? - future_orbit_tb.aop_deg()?
84    );
85    println!(
86        "TA changed by {:.3} deg",
87        orbit.ta_deg()? - future_orbit_tb.ta_deg()?
88    );
89
90    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
91    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
92    // models such as solar radiation pressure.
93
94    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
95    let sc = Spacecraft::builder()
96        .orbit(orbit)
97        .mass(Mass::from_dry_mass(9.60))
98        .srp(SRPData {
99            area_m2: 10e-4,
100            coeff_reflectivity: 1.1,
101        })
102        .build();
103    println!("{sc:x}");
104
105    // Set up the spacecraft dynamics.
106
107    // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
108    // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
109    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
110
111    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
112    // We're using the JGM3 model here, which is the default in GMAT.
113    let mut jgm3_meta = MetaFile {
114        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
115        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
116    };
117    // And let's download it if we don't have it yet.
118    jgm3_meta.process(true)?;
119
120    // Build the spherical harmonics.
121    // The harmonics must be computed in the body fixed frame.
122    // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
123    let harmonics_21x21 = Harmonics::from_stor(
124        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
125        HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
126    );
127
128    // Include the spherical harmonics into the orbital dynamics.
129    orbital_dyn.accel_models.push(harmonics_21x21);
130
131    // We define the solar radiation pressure, using the default solar flux and accounting only
132    // for the eclipsing caused by the Earth.
133    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
134
135    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
136    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
137    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
138
139    println!("{dynamics}");
140
141    // Finally, let's propagate this orbit to the same epoch as above.
142    // The first returned value is the spacecraft state at the final epoch.
143    // The second value is the full trajectory where the step size is variable step used by the propagator.
144    let (future_sc, trajectory) = Propagator::default(dynamics)
145        .with(sc, almanac.clone())
146        .until_epoch_with_traj(future_orbit_tb.epoch)?;
147
148    println!("=== High fidelity propagation ===");
149    println!(
150        "SMA changed by {:.3} km",
151        orbit.sma_km()? - future_sc.orbit.sma_km()?
152    );
153    println!(
154        "ECC changed by {:.6}",
155        orbit.ecc()? - future_sc.orbit.ecc()?
156    );
157    println!(
158        "INC changed by {:.3e} deg",
159        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
160    );
161    println!(
162        "RAAN changed by {:.3} deg",
163        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
164    );
165    println!(
166        "AOP changed by {:.3} deg",
167        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
168    );
169    println!(
170        "TA changed by {:.3} deg",
171        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
172    );
173
174    // We also have access to the full trajectory throughout the propagation.
175    println!("{trajectory}");
176
177    // With the trajectory, let's build a few data products.
178
179    // 1. Export the trajectory as a CCSDS OEM version 2.0 file and as a parquet file, which includes the Keplerian orbital elements.
180
181    trajectory.to_oem_file(
182        "./01_cubesat_hf_prop.oem",
183        ExportCfg::builder().step(Unit::Minute * 2).build(),
184    )?;
185
186    trajectory.to_parquet_with_cfg(
187        "./01_cubesat_hf_prop.parquet",
188        ExportCfg::builder().step(Unit::Minute * 2).build(),
189        almanac.clone(),
190    )?;
191
192    // 2. Compare the difference in the radial-intrack-crosstrack frame between the high fidelity
193    // and Keplerian propagation. The RIC frame is commonly used to compute the difference in position
194    // and velocity of different spacecraft.
195    // 3. Compute the azimuth, elevation, range, and range-rate data of that spacecraft as seen from Boulder, CO, USA.
196
197    let boulder_station = GroundStation::from_point(
198        "Boulder, CO, USA".to_string(),
199        40.014984,   // latitude in degrees
200        -105.270546, // longitude in degrees
201        1.6550,      // altitude in kilometers
202        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
203    );
204
205    // We iterate over the trajectory, grabbing a state every two minutes.
206    let mut offset_s = vec![];
207    let mut epoch_str = vec![];
208    let mut ric_x_km = vec![];
209    let mut ric_y_km = vec![];
210    let mut ric_z_km = vec![];
211    let mut ric_vx_km_s = vec![];
212    let mut ric_vy_km_s = vec![];
213    let mut ric_vz_km_s = vec![];
214
215    let mut azimuth_deg = vec![];
216    let mut elevation_deg = vec![];
217    let mut range_km = vec![];
218    let mut range_rate_km_s = vec![];
219    for state in trajectory.every(Unit::Minute * 2) {
220        // Try to compute the Keplerian/two body state just in time.
221        // This method occasionally fails to converge on an appropriate true anomaly
222        // from the mean anomaly. If that happens, we just skip this state.
223        // The high fidelity and Keplerian states diverge continuously, and we're curious
224        // about the divergence in this quick analysis.
225        let this_epoch = state.epoch();
226        match orbit.at_epoch(this_epoch) {
227            Ok(tb_then) => {
228                offset_s.push((this_epoch - orbit.epoch).to_seconds());
229                epoch_str.push(format!("{this_epoch}"));
230                // Compute the two body state just in time.
231                let ric = state.orbit.ric_difference(&tb_then)?;
232                ric_x_km.push(ric.radius_km.x);
233                ric_y_km.push(ric.radius_km.y);
234                ric_z_km.push(ric.radius_km.z);
235                ric_vx_km_s.push(ric.velocity_km_s.x);
236                ric_vy_km_s.push(ric.velocity_km_s.y);
237                ric_vz_km_s.push(ric.velocity_km_s.z);
238
239                // Compute the AER data for each state.
240                let aer = almanac.azimuth_elevation_range_sez(
241                    state.orbit,
242                    boulder_station.to_orbit(this_epoch, &almanac)?,
243                    None,
244                    None,
245                )?;
246                azimuth_deg.push(aer.azimuth_deg);
247                elevation_deg.push(aer.elevation_deg);
248                range_km.push(aer.range_km);
249                range_rate_km_s.push(aer.range_rate_km_s);
250            }
251            Err(e) => warn!("{} {e}", state.epoch()),
252        };
253    }
254
255    // Build the data frames.
256    let ric_df = df!(
257        "Offset (s)" => offset_s.clone(),
258        "Epoch" => epoch_str.clone(),
259        "RIC X (km)" => ric_x_km,
260        "RIC Y (km)" => ric_y_km,
261        "RIC Z (km)" => ric_z_km,
262        "RIC VX (km/s)" => ric_vx_km_s,
263        "RIC VY (km/s)" => ric_vy_km_s,
264        "RIC VZ (km/s)" => ric_vz_km_s,
265    )?;
266
267    println!("RIC difference at start\n{}", ric_df.head(Some(10)));
268    println!("RIC difference at end\n{}", ric_df.tail(Some(10)));
269
270    let aer_df = df!(
271        "Offset (s)" => offset_s.clone(),
272        "Epoch" => epoch_str.clone(),
273        "azimuth (deg)" => azimuth_deg,
274        "elevation (deg)" => elevation_deg,
275        "range (km)" => range_km,
276        "range rate (km/s)" => range_rate_km_s,
277    )?;
278
279    // Finally, let's see when the spacecraft is visible, assuming 15 degrees minimum elevation.
280    let mask = aer_df
281        .column("elevation (deg)")?
282        .gt(&Column::Scalar(ScalarColumn::new(
283            "elevation mask (deg)".into(),
284            Scalar::new(DataType::Float64, AnyValue::Float64(15.0)),
285            offset_s.len(),
286        )))?;
287    let cubesat_visible = aer_df.filter(&mask)?;
288
289    println!("{cubesat_visible}");
290
291    Ok(())
292}

pub fn try_keplerian_apsis_altitude( apo_alt_km: f64, peri_alt_km: f64, inc_deg: f64, raan_deg: f64, aop_deg: f64, ta_deg: f64, epoch: Epoch, frame: Frame, ) -> Result<CartesianState, PhysicsError>

Creates a new Orbit from the provided altitudes of apoapsis and periapsis, in kilometers

pub fn try_latlongalt( latitude_deg: f64, longitude_deg: f64, height_km: f64, angular_velocity_deg_s: f64, epoch: Epoch, frame: Frame, ) -> Result<CartesianState, PhysicsError>

Creates a new Orbit from the latitude (φ), longitude (λ) and height (in km) with respect to the frame’s ellipsoid given the angular velocity.

Note: The mean Earth angular velocity is 0.004178079012116429 deg/s.

NOTE: This computation differs from the spherical coordinates because we consider the flattening of body. Reference: G. Xu and Y. Xu, “GPS”, DOI 10.1007/978-3-662-50367-6_2, 2016

§

impl CartesianState

pub fn zero(frame: Frame) -> CartesianState

Builds a state of zero radius and velocity at zero seconds TDB (01 Jan 2000, midnight TDB) in the provided frame.

Examples found in repository?
examples/04_lro_od/main.rs (line 86)
33fn main() -> Result<(), Box<dyn Error>> {
34    pel::init();
35
36    // ====================== //
37    // === ALMANAC SET UP === //
38    // ====================== //
39
40    // Dynamics models require planetary constants and ephemerides to be defined.
41    // Let's start by grabbing those by using ANISE's MetaAlmanac.
42
43    let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
44        .iter()
45        .collect();
46
47    let meta = data_folder.join("lro-dynamics.dhall");
48
49    // Load this ephem in the general Almanac we're using for this analysis.
50    let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
51        .map_err(Box::new)?
52        .process(true)
53        .map_err(Box::new)?;
54
55    let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
56    moon_pc.mu_km3_s2 = 4902.74987;
57    almanac.planetary_data.set_by_id(MOON, moon_pc)?;
58
59    let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
60    earth_pc.mu_km3_s2 = 398600.436;
61    almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
62
63    // Save this new kernel for reuse.
64    // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
65    almanac
66        .planetary_data
67        .save_as(&data_folder.join("lro-specific.pca"), true)?;
68
69    // Lock the almanac (an Arc is a read only structure).
70    let almanac = Arc::new(almanac);
71
72    // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
73    // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
74    // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
75    // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
76    let lro_frame = Frame::from_ephem_j2000(-85);
77
78    // To build the trajectory we need to provide a spacecraft template.
79    let sc_template = Spacecraft::builder()
80        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
81        .srp(SRPData {
82            // SRP configuration is arbitrary, but we will be estimating it anyway.
83            area_m2: 3.9 * 2.7,
84            coeff_reflectivity: 0.96,
85        })
86        .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
87        .build();
88    // Now we can build the trajectory from the BSP file.
89    // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
90    let traj_as_flown = Traj::from_bsp(
91        lro_frame,
92        MOON_J2000,
93        almanac.clone(),
94        sc_template,
95        5.seconds(),
96        Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
97        Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
98        Aberration::LT,
99        Some("LRO".to_string()),
100    )?;
101
102    println!("{traj_as_flown}");
103
104    // ====================== //
105    // === MODEL MATCHING === //
106    // ====================== //
107
108    // Set up the spacecraft dynamics.
109
110    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
111    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
112    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
113
114    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
115    // We're using the GRAIL JGGRX model.
116    let mut jggrx_meta = MetaFile {
117        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
118        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
119    };
120    // And let's download it if we don't have it yet.
121    jggrx_meta.process(true)?;
122
123    // Build the spherical harmonics.
124    // The harmonics must be computed in the body fixed frame.
125    // We're using the long term prediction of the Moon principal axes frame.
126    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
127    let sph_harmonics = Harmonics::from_stor(
128        almanac.frame_from_uid(moon_pa_frame)?,
129        HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
130    );
131
132    // Include the spherical harmonics into the orbital dynamics.
133    orbital_dyn.accel_models.push(sph_harmonics);
134
135    // We define the solar radiation pressure, using the default solar flux and accounting only
136    // for the eclipsing caused by the Earth and Moon.
137    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
138    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
139
140    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
141    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
142    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
143
144    println!("{dynamics}");
145
146    // Now we can build the propagator.
147    let setup = Propagator::default_dp78(dynamics.clone());
148
149    // For reference, let's build the trajectory with Nyx's models from that LRO state.
150    let (sim_final, traj_as_sim) = setup
151        .with(*traj_as_flown.first(), almanac.clone())
152        .until_epoch_with_traj(traj_as_flown.last().epoch())?;
153
154    println!("SIM INIT:  {:x}", traj_as_flown.first());
155    println!("SIM FINAL: {sim_final:x}");
156    // Compute RIC difference between SIM and LRO ephem
157    let sim_lro_delta = sim_final
158        .orbit
159        .ric_difference(&traj_as_flown.last().orbit)?;
160    println!("{traj_as_sim}");
161    println!(
162        "SIM v LRO - RIC Position (m): {:.3}",
163        sim_lro_delta.radius_km * 1e3
164    );
165    println!(
166        "SIM v LRO - RIC Velocity (m/s): {:.3}",
167        sim_lro_delta.velocity_km_s * 1e3
168    );
169
170    traj_as_sim.ric_diff_to_parquet(
171        &traj_as_flown,
172        "./04_lro_sim_truth_error.parquet",
173        ExportCfg::default(),
174    )?;
175
176    // ==================== //
177    // === OD SIMULATOR === //
178    // ==================== //
179
180    // After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
181    // and the truth LRO state.
182
183    // Therefore, we will actually run an estimation from a dispersed LRO state.
184    // The sc_seed is the true LRO state from the BSP.
185    let sc_seed = *traj_as_flown.first();
186
187    // Load the Deep Space Network ground stations.
188    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
189    let ground_station_file: PathBuf = [
190        env!("CARGO_MANIFEST_DIR"),
191        "examples",
192        "04_lro_od",
193        "dsn-network.yaml",
194    ]
195    .iter()
196    .collect();
197
198    let devices = GroundStation::load_named(ground_station_file)?;
199
200    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
201    // Nyx can build a tracking schedule for you based on the first station with access.
202    let trkconfg_yaml: PathBuf = [
203        env!("CARGO_MANIFEST_DIR"),
204        "examples",
205        "04_lro_od",
206        "tracking-cfg.yaml",
207    ]
208    .iter()
209    .collect();
210
211    let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
212
213    // Build the tracking arc simulation to generate a "standard measurement".
214    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
215        devices.clone(),
216        traj_as_flown.clone(),
217        configs,
218    )?;
219
220    trk.build_schedule(almanac.clone())?;
221    let arc = trk.generate_measurements(almanac.clone())?;
222    // Save the simulated tracking data
223    arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
224
225    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
226    println!("{arc}");
227
228    // Now that we have simulated measurements, we'll run the orbit determination.
229
230    // ===================== //
231    // === OD ESTIMATION === //
232    // ===================== //
233
234    let sc = SpacecraftUncertainty::builder()
235        .nominal(sc_seed)
236        .frame(LocalFrame::RIC)
237        .x_km(0.5)
238        .y_km(0.5)
239        .z_km(0.5)
240        .vx_km_s(5e-3)
241        .vy_km_s(5e-3)
242        .vz_km_s(5e-3)
243        .build();
244
245    // Build the filter initial estimate, which we will reuse in the filter.
246    let initial_estimate = sc.to_estimate()?;
247
248    println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
249
250    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
251    let process_noise = ProcessNoise3D::from_velocity_km_s(
252        &[1.8e-9, 1.8e-9, 1.8e-9],
253        1 * Unit::Hour,
254        10 * Unit::Minute,
255        None,
256    );
257
258    println!("{process_noise}");
259
260    let kf = KF::new(
261        // Increase the initial covariance to account for larger deviation.
262        initial_estimate,
263        process_noise,
264    );
265
266    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
267    let mut odp = SpacecraftODProcess::ekf(
268        setup.with(initial_estimate.state().with_stm(), almanac.clone()),
269        kf,
270        devices,
271        EkfTrigger::new(0, Unit::Hour * 1),
272        Some(ResidRejectCrit::default()),
273        almanac.clone(),
274    );
275
276    odp.process_arc(&arc)?;
277
278    let ric_err = traj_as_flown
279        .at(odp.estimates.last().unwrap().epoch())?
280        .orbit
281        .ric_difference(&odp.estimates.last().unwrap().orbital_state())?;
282    println!("== RIC at end ==");
283    println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
284    println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
285
286    odp.to_parquet(&arc, "./04_lro_od_results.parquet", ExportCfg::default())?;
287
288    // In our case, we have the truth trajectory from NASA.
289    // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
290    // Export the OD trajectory first.
291    let od_trajectory = odp.to_traj()?;
292    // Build the RIC difference.
293    od_trajectory.ric_diff_to_parquet(
294        &traj_as_flown,
295        "./04_lro_od_truth_error.parquet",
296        ExportCfg::default(),
297    )?;
298
299    Ok(())
300}

pub fn zero_at_epoch(epoch: Epoch, frame: Frame) -> CartesianState

Builds a state of zero radius and velocity at the provided epoch in the provided frame.

pub fn new( x_km: f64, y_km: f64, z_km: f64, vx_km_s: f64, vy_km_s: f64, vz_km_s: f64, epoch: Epoch, frame: Frame, ) -> CartesianState

Creates a new Cartesian state in the provided frame at the provided Epoch.

Units: km, km, km, km/s, km/s, km/s

pub fn cartesian( x_km: f64, y_km: f64, z_km: f64, vx_km_s: f64, vy_km_s: f64, vz_km_s: f64, epoch: Epoch, frame: Frame, ) -> CartesianState

Creates a new Cartesian state in the provided frame at the provided Epoch (shortcut to new).

Units: km, km, km, km/s, km/s, km/s

pub fn from_position( x_km: f64, y_km: f64, z_km: f64, epoch: Epoch, frame: Frame, ) -> CartesianState

Creates a new Cartesian in the provided frame at the provided Epoch in time with 0.0 velocity.

Units: km, km, km

pub fn from_cartesian_pos_vel( pos_vel: Matrix<f64, Const<6>, Const<1>, ArrayStorage<f64, 6, 1>>, epoch: Epoch, frame: Frame, ) -> CartesianState

Creates a new Cartesian around in the provided frame from the borrowed state vector

The state vector must be x, y, z, vx, vy, vz. This function is a shortcut to cartesian and as such it has the same unit requirements.

Units: position data must be in kilometers, velocity data must be in kilometers per second.

pub fn with_radius_km( self, new_radius_km: Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>, ) -> CartesianState

Returns a copy of the state with a new radius

pub fn with_velocity_km_s( self, new_velocity_km_s: Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>, ) -> CartesianState

Returns a copy of the state with a new radius

pub fn to_cartesian_pos_vel( self, ) -> Matrix<f64, Const<6>, Const<1>, ArrayStorage<f64, 6, 1>>

Returns this state as a Cartesian Vector6 in [km, km, km, km/s, km/s, km/s]

Note that the time is not returned in the vector.

pub fn with_cartesian_pos_vel( self, pos_vel: Matrix<f64, Const<6>, Const<1>, ArrayStorage<f64, 6, 1>>, ) -> CartesianState

Returns a copy of this state where the position and velocity are set to the input vector whose units must be [km, km, km, km/s, km/s, km/s]

pub fn distance_to_point_km( &self, other_km: &Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>, ) -> f64

Returns the distance in kilometers between this state and a point assumed to be in the same frame.

pub fn r_hat(&self) -> Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>

Returns the unit vector in the direction of the state radius

pub fn v_hat(&self) -> Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>

Returns the unit vector in the direction of the state velocity

pub fn apply_dv_km_s( &mut self, dv_km_s: Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>, )

Adds the provided delta-v (in km/s) to the current velocity vector, mimicking an impulsive maneuver.

pub fn with_dv_km_s( &self, dv_km_s: Matrix<f64, Const<3>, Const<1>, ArrayStorage<f64, 3, 1>>, ) -> CartesianState

Copies this orbit after adding the provided delta-v (in km/s) to the velocity vector, mimicking an impulsive maneuver.

pub fn has_velocity_dynamics(&self) -> bool

Returns True if velocity (dynamics) are defined. Which may not be the case if [Self] was built with [Self::from_position] and you only intend to use partial rotations.

Trait Implementations

§

impl Add for CartesianState

§

fn add(self, other: CartesianState) -> <CartesianState as Add>::Output

Adds one state to another. This will return an error if the epochs or frames are different.

§

type Output = Result<CartesianState, PhysicsError>

The resulting type after applying the + operator.
§

impl Clone for CartesianState

§

fn clone(&self) -> CartesianState

Returns a copy of the value. Read more
1.0.0 · Source§

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source. Read more
§

impl Debug for CartesianState

§

fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error>

Formats the value using the given formatter. Read more
§

impl<'de> Deserialize<'de> for CartesianState

§

fn deserialize<__D>( __deserializer: __D, ) -> Result<CartesianState, <__D as Deserializer<'de>>::Error>
where __D: Deserializer<'de>,

Deserialize this value from the given Serde deserializer. Read more
§

impl Display for CartesianState

§

fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error>

Formats the value using the given formatter. Read more
§

impl LowerExp for CartesianState

§

fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error>

Formats the value using the given formatter. Read more
§

impl LowerHex for CartesianState

§

fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error>

Prints the Keplerian orbital elements in floating point with units if frame is celestial, If frame is geodetic, prints the range, altitude, latitude, and longitude with respect to the planetocentric frame

§

impl Neg for CartesianState

§

type Output = CartesianState

The resulting type after applying the - operator.
§

fn neg(self) -> <CartesianState as Neg>::Output

Performs the unary - operation. Read more
§

impl PartialEq for CartesianState

§

fn eq(&self, other: &CartesianState) -> bool

Two states are equal if their position are equal within one centimeter and their velocities within one centimeter per second.

1.0.0 · Source§

fn ne(&self, other: &Rhs) -> bool

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.
§

impl Serialize for CartesianState

§

fn serialize<__S>( &self, __serializer: __S, ) -> Result<<__S as Serializer>::Ok, <__S as Serializer>::Error>
where __S: Serializer,

Serialize this value into the given Serde serializer. Read more
§

impl Sub for CartesianState

§

fn sub(self, other: CartesianState) -> <CartesianState as Sub>::Output

Adds one state to another. This will return an error if the epochs or frames are different.

§

type Output = Result<CartesianState, PhysicsError>

The resulting type after applying the - operator.
§

impl UpperHex for CartesianState

§

fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error>

Prints the prints the range, altitude, latitude, and longitude with respect to the planetocentric frame in floating point with units if frame is celestial,

§

impl Copy for CartesianState