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
impl CartesianState
pub fn dcm_from_topocentric_to_body_fixed(
&self,
_from: i32,
) -> Result<DCM, PhysicsError>
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>
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>
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
- Compute the state data one millisecond before and one millisecond assuming two body dynamics
- Compute the DCM for this state, and the pre and post states
- Build the c vector as the normalized orbital momentum vector
- Build the i vector as the cross product of \hat{r} and c
- Build the RIC DCM as a 3x3 of the columns [\hat{r}, \hat{i}, \hat{c}], for the post, post, and current states
- Compute the difference between the DCMs of the pre and post states, to build the DCM time derivative
- 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>
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
- Build the c vector as the normalized orbital momentum vector
- Build the i vector as the cross product of \hat{r} and c
- Build the RIC DCM as a 3x3 of the columns [\hat{r}, \hat{i}, \hat{c}]
- Return the DCM structure without accounting for the transport theorem.
pub fn dcm3x3_from_rcn_to_inertial(&self) -> Result<DCM, PhysicsError>
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
- Compute \hat{r}, \hat{h}, the unit vectors of the radius and orbital momentum.
- Compute the cross product of these
- Build the DCM with these unit vectors
- Return the DCM structure
pub fn dcm_from_rcn_to_inertial(&self) -> Result<DCM, PhysicsError>
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
- Compute \hat{r}, \hat{h}, the unit vectors of the radius and orbital momentum.
- Compute the cross product of these
- Build the DCM with these unit vectors
- 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>
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
- Compute \hat{v}, \hat{h}, the unit vectors of the radius and orbital momentum.
- Compute the cross product of these
- Build the DCM with these unit vectors
- Return the DCM structure.
pub fn dcm_from_vnc_to_inertial(&self) -> Result<DCM, PhysicsError>
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
- Compute \hat{v}, \hat{h}, the unit vectors of the radius and orbital momentum.
- Compute the cross product of these
- Build the DCM with these unit vectors
- Compute the difference between the DCMs of the pre and post states (+/- 1 ms), to build the DCM time derivative
- 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>
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>
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>
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>
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>
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>
pub fn sma_km(&self) -> Result<f64, PhysicsError>
Returns the semi-major axis in km
:rtype: float
Examples found in repository?
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
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>
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>
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>
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>
pub fn period(&self) -> Result<Duration, PhysicsError>
Returns the period in seconds
:rtype: Duration
pub fn ecc(&self) -> Result<f64, PhysicsError>
pub fn ecc(&self) -> Result<f64, PhysicsError>
Returns the eccentricity (no unit)
:rtype: float
Examples found in repository?
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
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>
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>
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>
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>
pub fn inc_deg(&self) -> Result<f64, PhysicsError>
Returns the inclination in degrees
:rtype: float
Examples found in repository?
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
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>
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>
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>
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>
pub fn aop_deg(&self) -> Result<f64, PhysicsError>
Returns the argument of periapsis in degrees
:rtype: float
Examples found in repository?
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
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>
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>
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>
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>
pub fn raan_deg(&self) -> Result<f64, PhysicsError>
Returns the right ascension of the ascending node in degrees
:rtype: float
Examples found in repository?
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
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>
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>
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>
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>
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?
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
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>
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>
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>
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>
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>
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>
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>
pub fn tlong_deg(&self) -> Result<f64, PhysicsError>
Returns the true longitude in degrees
:rtype: float
pub fn aol_deg(&self) -> Result<f64, PhysicsError>
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>
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>
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>
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>
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>
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>
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>
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
pub fn right_ascension_deg(&self) -> f64
Returns the right ascension of this orbit in degrees
:rtype: float
pub fn declination_deg(&self) -> f64
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>
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
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>
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>
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>
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>
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>
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?
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>
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
- Compute the RIC DCM of self
- Rotate self into the RIC frame
- Rotation other into the RIC frame
- Compute the difference between these two states
- Strip the astrodynamical information from the frame, enabling only computations from
CartesianState
:type other: Orbit :rtype: Orbit
Examples found in repository?
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
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>
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
- Compute the VNC DCM of self
- Rotate self into the VNC frame
- Rotation other into the VNC frame
- Compute the difference between these two states
- Strip the astrodynamical information from the frame, enabling only computations from
CartesianState
:type other: Orbit :rtype: Orbit
§impl CartesianState
impl CartesianState
pub fn distance_to_km(
&self,
other: &CartesianState,
) -> Result<f64, PhysicsError>
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>
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>
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>
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>
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
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
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>
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>
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>
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>
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>
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>
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
impl CartesianState
pub fn sma_altitude_km(&self) -> Result<f64, PhysicsError>
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>
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>
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>
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?
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
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
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>
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
§impl CartesianState
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>
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?
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
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>
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>
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
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?
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
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>
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
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>
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.
§impl CartesianState
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>
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?
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>
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>
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
impl CartesianState
pub fn zero(frame: Frame) -> 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?
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
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
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
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
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
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
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
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>>
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
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
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>>
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>>
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>>,
)
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
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
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.