pub struct EclipseLocator {
pub light_source: Frame,
pub shadow_bodies: Vec<Frame>,
}
Fields§
§light_source: Frame
§shadow_bodies: Vec<Frame>
Implementations§
Source§impl EclipseLocator
impl EclipseLocator
Sourcepub fn cislunar(almanac: Arc<Almanac>) -> Self
pub fn cislunar(almanac: Arc<Almanac>) -> Self
Creates a new typical eclipse locator. The light source is the Sun, and the shadow bodies are the Earth and the Moon.
Examples found in repository?
examples/03_geo_analysis/stationkeeping.rs (line 122)
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// Set up the dynamics like in the orbit raise.
let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
// Define the GEO orbit, and we're just going to maintain it very tightly.
let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
println!("{orbit:x}");
let sc = Spacecraft::builder()
.orbit(orbit)
.dry_mass_kg(1000.0) // 1000 kg of dry mass
.fuel_mass_kg(1000.0) // 1000 kg of fuel, totalling 2.0 tons
.srp(SrpConfig::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
.thruster(Thruster {
// "NEXT-STEP" row in Table 2
isp_s: 4435.0,
thrust_N: 0.472,
})
.mode(GuidanceMode::Thrust) // Start thrusting immediately.
.build();
// Set up the spacecraft dynamics like in the orbit raise example.
let prop_time = 30.0 * Unit::Day;
// Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
let objectives = &[
Objective::within_tolerance(StateParameter::SMA, 42_164.0, 5.0), // 5 km
Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5),
Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2),
];
let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2)?;
println!("{ruggiero_ctrl}");
let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
let mut jgm3_meta = MetaFile {
uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
};
jgm3_meta.process(true)?;
let harmonics = Harmonics::from_stor(
almanac.frame_from_uid(IAU_EARTH_FRAME)?,
HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true)?,
);
orbital_dyn.accel_models.push(harmonics);
let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
.with_guidance_law(ruggiero_ctrl.clone());
println!("{sc_dynamics}");
// Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.
// Let's start by defining the dispersion.
// 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.
// Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
let mc_rv = MultivariateNormal::new(
sc,
vec![StateDispersion::zero_mean(StateParameter::SMA, 3.0)],
)?;
let my_mc = MonteCarlo::new(
sc, // Nominal state
mc_rv,
"03_geo_sk".to_string(), // Scenario name
None, // No specific seed specified, so one will be drawn from the computer's entropy.
);
// Build the propagator setup.
let setup = Propagator::rk89(
sc_dynamics.clone(),
IntegratorOptions::builder()
.min_step(10.0_f64.seconds())
.error_ctrl(ErrorControl::RSSCartesianStep)
.build(),
);
let num_runs = 25;
let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);
assert_eq!(rslts.runs.len(), num_runs);
// For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
rslts.to_parquet(
"03_geo_sk.parquet",
Some(vec![
&EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
]),
ExportCfg::default(),
almanac,
)?;
Ok(())
}
More examples
examples/03_geo_analysis/raise.rs (line 136)
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's latest MetaAlmanac.
// This will automatically download the DE440s planetary ephemeris,
// the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
// parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
// planetary constants kernels.
// For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
// Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
// references to many functions.
let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
// Fetch the EME2000 frame from the Almabac
let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
// Define the orbit epoch
let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
// Build the spacecraft itself.
// Using slide 6 of https://aerospace.org/sites/default/files/2018-11/Davis-Mayberry_HPSEP_11212018.pdf
// for the "next gen" SEP characteristics.
// GTO start
let orbit = Orbit::keplerian(24505.9, 0.725, 7.05, 0.0, 0.0, 0.0, epoch, eme2k);
let sc = Spacecraft::builder()
.orbit(orbit)
.dry_mass_kg(1000.0) // 1000 kg of dry mass
.fuel_mass_kg(1000.0) // 1000 kg of fuel, totalling 2.0 tons
.srp(SrpConfig::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
.thruster(Thruster {
// "NEXT-STEP" row in Table 2
isp_s: 4435.0,
thrust_N: 0.472,
})
.mode(GuidanceMode::Thrust) // Start thrusting immediately.
.build();
let prop_time = 180.0 * Unit::Day;
// Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
let objectives = &[
Objective::within_tolerance(StateParameter::SMA, 42_165.0, 20.0),
Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5),
Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2),
];
// Ensure that we only thrust if we have more than 20% illumination.
let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2).unwrap();
println!("{ruggiero_ctrl}");
// Define the high fidelity dynamics
// Set up the spacecraft dynamics.
// Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
// The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
// We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
// We're using the JGM3 model here, which is the default in GMAT.
let mut jgm3_meta = MetaFile {
uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
};
// And let's download it if we don't have it yet.
jgm3_meta.process(true)?;
// Build the spherical harmonics.
// The harmonics must be computed in the body fixed frame.
// We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
let harmonics = Harmonics::from_stor(
almanac.frame_from_uid(IAU_EARTH_FRAME)?,
HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true).unwrap(),
);
// Include the spherical harmonics into the orbital dynamics.
orbital_dyn.accel_models.push(harmonics);
// We define the solar radiation pressure, using the default solar flux and accounting only
// for the eclipsing caused by the Earth.
let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
// Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
// acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
.with_guidance_law(ruggiero_ctrl.clone());
println!("{:x}", orbit);
// We specify a minimum step in the propagator because the Ruggiero control would otherwise drive this step very low.
let (final_state, traj) = Propagator::rk89(
sc_dynamics.clone(),
IntegratorOptions::builder()
.min_step(10.0_f64.seconds())
.error_ctrl(ErrorControl::RSSCartesianStep)
.build(),
)
.with(sc, almanac.clone())
.for_duration_with_traj(prop_time)?;
let fuel_usage = sc.fuel_mass_kg - final_state.fuel_mass_kg;
println!("{:x}", final_state.orbit);
println!("fuel usage: {:.3} kg", fuel_usage);
// Finally, export the results for analysis, including the penumbra percentage throughout the orbit raise.
traj.to_parquet(
"./03_geo_raise.parquet",
Some(vec![
&EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
]),
ExportCfg::default(),
almanac,
)?;
for status_line in ruggiero_ctrl.status(&final_state) {
println!("{status_line}");
}
ruggiero_ctrl
.achieved(&final_state)
.expect("objective not achieved");
Ok(())
}
examples/02_jwst_covar_monte_carlo/main.rs (line 149)
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's latest MetaAlmanac.
// For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
// Download the regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
// Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
let mut latest_jwst_ephem = MetaFile {
uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
crc32: None,
};
latest_jwst_ephem.process(true)?;
// Load this ephem in the general Almanac we're using for this analysis.
let almanac = Arc::new(
MetaAlmanac::latest()
.map_err(Box::new)?
.load_from_metafile(latest_jwst_ephem, true)?,
);
// By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
// in the BSP. We need this ID in order to query the ephemeris.
const JWST_NAIF_ID: i32 = -170;
// Let's build a frame in the J2000 orientation centered on the JWST.
const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
// Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
println!("JWST defined from {earliest_epoch} to {latest_epoch}");
// Fetch the state, printing it in the Earth J2000 frame.
let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
println!("{jwst_orbit:x}");
// Build the spacecraft
// SRP area assumed to be the full sunshield and mass if 6200.0 kg, c.f. https://webb.nasa.gov/content/about/faqs/facts.html
// SRP Coefficient of reflectivity assumed to be that of Kapton, i.e. 2 - 0.44 = 1.56, table 1 from https://amostech.com/TechnicalPapers/2018/Poster/Bengtson.pdf
let jwst = Spacecraft::builder()
.orbit(jwst_orbit)
.srp(SrpConfig {
area_m2: 21.197 * 14.162,
cr: 1.56,
})
.dry_mass_kg(6200.0)
.build();
// Build up the spacecraft uncertainty builder.
// We can use the spacecraft uncertainty structure to build this up.
// We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
// in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
// Nyx can also estimate the deviation of the spacecraft parameters.
let jwst_uncertainty = SpacecraftUncertainty::builder()
.nominal(jwst)
.frame(LocalFrame::RIC)
.x_km(0.5)
.y_km(0.3)
.z_km(1.5)
.vx_km_s(1e-4)
.vy_km_s(0.6e-3)
.vz_km_s(3e-3)
.build();
println!("{jwst_uncertainty}");
// Build the Kalman filter estimate.
// Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
// but this approach requires quite a bit more boilerplate code.
let jwst_estimate = jwst_uncertainty.to_estimate()?;
// Set up the spacecraft dynamics.
// We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
// We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
// Build the propagator set up to use for the whole analysis.
let setup = Propagator::default(dynamics);
// All of the analysis will use this duration.
let prediction_duration = 6.5 * Unit::Day;
// === Covariance mapping ===
// For the covariance mapping / prediction, we'll use the common orbit determination approach.
// This is done by setting up a spacecraft OD process, and predicting for the analysis duration.
let ckf = KF::no_snc(jwst_estimate);
// Build the propagation instance for the OD process.
let prop = setup.with(jwst.with_stm(), almanac.clone());
let mut odp = SpacecraftODProcess::ckf(prop, ckf, None, almanac.clone());
// Define the prediction step, i.e. how often we want to know the covariance.
let step = 1_i64.minutes();
// Finally, predict, and export the trajectory with covariance to a parquet file.
odp.predict_for(step, prediction_duration)?;
odp.to_parquet("./02_jwst_covar_map.parquet", ExportCfg::default())?;
// === Monte Carlo framework ===
// Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
let my_mc = MonteCarlo::new(
jwst, // Nominal state
jwst_estimate.to_random_variable()?,
"02_jwst".to_string(), // Scenario name
None, // No specific seed specified, so one will be drawn from the computer's entropy.
);
let num_runs = 5_000;
let rslts = my_mc.run_until_epoch(
setup,
almanac.clone(),
jwst.epoch() + prediction_duration,
num_runs,
);
assert_eq!(rslts.runs.len(), num_runs);
// Finally, export these results, computing the eclipse percentage for all of these results.
// For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
let umbra_event = eclipse_loc.to_umbra_event();
let penumbra_event = eclipse_loc.to_penumbra_event();
rslts.to_parquet(
"02_jwst_monte_carlo.parquet",
Some(vec![&umbra_event, &penumbra_event]),
ExportCfg::default(),
almanac,
)?;
Ok(())
}
examples/03_geo_analysis/drift.rs (line 155)
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's latest MetaAlmanac.
// This will automatically download the DE440s planetary ephemeris,
// the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
// parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
// planetary constants kernels.
// For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
// Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
// references to many functions.
let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
// Define the orbit epoch
let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
// Define the orbit.
// First we need to fetch the Earth J2000 from information from the Almanac.
// This allows the frame to include the gravitational parameters and the shape of the Earth,
// defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
// by loading a different set of planetary constants.
let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
// Placing this GEO bird just above Colorado.
// In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
// Print in in Keplerian form.
println!("{orbit:x}");
let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;
// Nyx is used for high fidelity propagation, not Keplerian propagation as above.
// Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
// models such as solar radiation pressure.
// Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
let sc = Spacecraft::builder()
.orbit(orbit)
.dry_mass_kg(9.60)
.srp(SrpConfig {
area_m2: 10e-4,
cr: 1.1,
})
.build();
println!("{sc:x}");
// Set up the spacecraft dynamics.
// Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
// The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
// We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
// We're using the JGM3 model here, which is the default in GMAT.
let mut jgm3_meta = MetaFile {
uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
};
// And let's download it if we don't have it yet.
jgm3_meta.process(true)?;
// Build the spherical harmonics.
// The harmonics must be computed in the body fixed frame.
// We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
let harmonics_21x21 = Harmonics::from_stor(
almanac.frame_from_uid(IAU_EARTH_FRAME)?,
HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
);
// Include the spherical harmonics into the orbital dynamics.
orbital_dyn.accel_models.push(harmonics_21x21);
// We define the solar radiation pressure, using the default solar flux and accounting only
// for the eclipsing caused by the Earth and Moon.
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
// acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
println!("{dynamics}");
// Finally, let's propagate this orbit to the same epoch as above.
// The first returned value is the spacecraft state at the final epoch.
// The second value is the full trajectory where the step size is variable step used by the propagator.
let (future_sc, trajectory) = Propagator::default(dynamics)
.with(sc, almanac.clone())
.until_epoch_with_traj(epoch + Unit::Century * 0.03)?;
println!("=== High fidelity propagation ===");
println!(
"SMA changed by {:.3} km",
orbit.sma_km()? - future_sc.orbit.sma_km()?
);
println!(
"ECC changed by {:.6}",
orbit.ecc()? - future_sc.orbit.ecc()?
);
println!(
"INC changed by {:.3e} deg",
orbit.inc_deg()? - future_sc.orbit.inc_deg()?
);
println!(
"RAAN changed by {:.3} deg",
orbit.raan_deg()? - future_sc.orbit.raan_deg()?
);
println!(
"AOP changed by {:.3} deg",
orbit.aop_deg()? - future_sc.orbit.aop_deg()?
);
println!(
"TA changed by {:.3} deg",
orbit.ta_deg()? - future_sc.orbit.ta_deg()?
);
// We also have access to the full trajectory throughout the propagation.
println!("{trajectory}");
println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");
// With the trajectory, let's build a few data products.
// 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.
let analysis_step = Unit::Minute * 5;
trajectory.to_parquet(
"./03_geo_hf_prop.parquet",
Some(vec![
&EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
]),
ExportCfg::builder().step(analysis_step).build(),
almanac.clone(),
)?;
// 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.
// We iterate over the trajectory, grabbing a state every two minutes.
let mut offset_s = vec![];
let mut epoch_str = vec![];
let mut longitude_deg = vec![];
let mut latitude_deg = vec![];
let mut altitude_km = vec![];
for state in trajectory.every(analysis_step) {
// Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
// These define the GEO stationkeeping box.
let this_epoch = state.epoch();
offset_s.push((this_epoch - orbit.epoch).to_seconds());
epoch_str.push(this_epoch.to_isoformat());
let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
longitude_deg.push(long_deg);
latitude_deg.push(lat_deg);
altitude_km.push(alt_km);
}
println!(
"Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
orig_long_deg - longitude_deg.last().unwrap()
);
println!(
"Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
orig_lat_deg - latitude_deg.last().unwrap()
);
println!(
"Altitude changed by {:.3} km -- Box is 30 km",
orig_alt_km - altitude_km.last().unwrap()
);
// Build the station keeping data frame.
let mut sk_df = df!(
"Offset (s)" => offset_s.clone(),
"Epoch (UTC)" => epoch_str.clone(),
"Longitude E-W (deg)" => longitude_deg,
"Latitude N-S (deg)" => latitude_deg,
"Altitude (km)" => altitude_km,
)?;
// Create a file to write the Parquet to
let file = File::create("./03_geo_lla.parquet").expect("Could not create file");
// Create a ParquetWriter and write the DataFrame to the file
ParquetWriter::new(file).finish(&mut sk_df)?;
Ok(())
}
Sourcepub fn compute(
&self,
observer: Orbit,
almanac: Arc<Almanac>,
) -> AlmanacResult<Occultation>
pub fn compute( &self, observer: Orbit, almanac: Arc<Almanac>, ) -> AlmanacResult<Occultation>
Compute the visibility/eclipse between an observer and an observed state
Sourcepub fn to_umbra_event(&self) -> UmbraEvent
pub fn to_umbra_event(&self) -> UmbraEvent
Creates an umbra event from this eclipse locator. Evaluation of the event, returns 0.0 for umbra, 1.0 for visibility (no shadow) and some value in between for penumbra
Examples found in repository?
examples/02_jwst_covar_monte_carlo/main.rs (line 150)
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's latest MetaAlmanac.
// For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
// Download the regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
// Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
let mut latest_jwst_ephem = MetaFile {
uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
crc32: None,
};
latest_jwst_ephem.process(true)?;
// Load this ephem in the general Almanac we're using for this analysis.
let almanac = Arc::new(
MetaAlmanac::latest()
.map_err(Box::new)?
.load_from_metafile(latest_jwst_ephem, true)?,
);
// By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
// in the BSP. We need this ID in order to query the ephemeris.
const JWST_NAIF_ID: i32 = -170;
// Let's build a frame in the J2000 orientation centered on the JWST.
const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
// Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
println!("JWST defined from {earliest_epoch} to {latest_epoch}");
// Fetch the state, printing it in the Earth J2000 frame.
let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
println!("{jwst_orbit:x}");
// Build the spacecraft
// SRP area assumed to be the full sunshield and mass if 6200.0 kg, c.f. https://webb.nasa.gov/content/about/faqs/facts.html
// SRP Coefficient of reflectivity assumed to be that of Kapton, i.e. 2 - 0.44 = 1.56, table 1 from https://amostech.com/TechnicalPapers/2018/Poster/Bengtson.pdf
let jwst = Spacecraft::builder()
.orbit(jwst_orbit)
.srp(SrpConfig {
area_m2: 21.197 * 14.162,
cr: 1.56,
})
.dry_mass_kg(6200.0)
.build();
// Build up the spacecraft uncertainty builder.
// We can use the spacecraft uncertainty structure to build this up.
// We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
// in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
// Nyx can also estimate the deviation of the spacecraft parameters.
let jwst_uncertainty = SpacecraftUncertainty::builder()
.nominal(jwst)
.frame(LocalFrame::RIC)
.x_km(0.5)
.y_km(0.3)
.z_km(1.5)
.vx_km_s(1e-4)
.vy_km_s(0.6e-3)
.vz_km_s(3e-3)
.build();
println!("{jwst_uncertainty}");
// Build the Kalman filter estimate.
// Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
// but this approach requires quite a bit more boilerplate code.
let jwst_estimate = jwst_uncertainty.to_estimate()?;
// Set up the spacecraft dynamics.
// We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
// We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
// Build the propagator set up to use for the whole analysis.
let setup = Propagator::default(dynamics);
// All of the analysis will use this duration.
let prediction_duration = 6.5 * Unit::Day;
// === Covariance mapping ===
// For the covariance mapping / prediction, we'll use the common orbit determination approach.
// This is done by setting up a spacecraft OD process, and predicting for the analysis duration.
let ckf = KF::no_snc(jwst_estimate);
// Build the propagation instance for the OD process.
let prop = setup.with(jwst.with_stm(), almanac.clone());
let mut odp = SpacecraftODProcess::ckf(prop, ckf, None, almanac.clone());
// Define the prediction step, i.e. how often we want to know the covariance.
let step = 1_i64.minutes();
// Finally, predict, and export the trajectory with covariance to a parquet file.
odp.predict_for(step, prediction_duration)?;
odp.to_parquet("./02_jwst_covar_map.parquet", ExportCfg::default())?;
// === Monte Carlo framework ===
// Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
let my_mc = MonteCarlo::new(
jwst, // Nominal state
jwst_estimate.to_random_variable()?,
"02_jwst".to_string(), // Scenario name
None, // No specific seed specified, so one will be drawn from the computer's entropy.
);
let num_runs = 5_000;
let rslts = my_mc.run_until_epoch(
setup,
almanac.clone(),
jwst.epoch() + prediction_duration,
num_runs,
);
assert_eq!(rslts.runs.len(), num_runs);
// Finally, export these results, computing the eclipse percentage for all of these results.
// For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
let umbra_event = eclipse_loc.to_umbra_event();
let penumbra_event = eclipse_loc.to_penumbra_event();
rslts.to_parquet(
"02_jwst_monte_carlo.parquet",
Some(vec![&umbra_event, &penumbra_event]),
ExportCfg::default(),
almanac,
)?;
Ok(())
}
Sourcepub fn to_penumbra_event(&self) -> PenumbraEvent
pub fn to_penumbra_event(&self) -> PenumbraEvent
Creates a penumbra event from this eclipse locator
Examples found in repository?
examples/03_geo_analysis/stationkeeping.rs (line 122)
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// Set up the dynamics like in the orbit raise.
let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
// Define the GEO orbit, and we're just going to maintain it very tightly.
let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
println!("{orbit:x}");
let sc = Spacecraft::builder()
.orbit(orbit)
.dry_mass_kg(1000.0) // 1000 kg of dry mass
.fuel_mass_kg(1000.0) // 1000 kg of fuel, totalling 2.0 tons
.srp(SrpConfig::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
.thruster(Thruster {
// "NEXT-STEP" row in Table 2
isp_s: 4435.0,
thrust_N: 0.472,
})
.mode(GuidanceMode::Thrust) // Start thrusting immediately.
.build();
// Set up the spacecraft dynamics like in the orbit raise example.
let prop_time = 30.0 * Unit::Day;
// Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
let objectives = &[
Objective::within_tolerance(StateParameter::SMA, 42_164.0, 5.0), // 5 km
Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5),
Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2),
];
let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2)?;
println!("{ruggiero_ctrl}");
let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
let mut jgm3_meta = MetaFile {
uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
};
jgm3_meta.process(true)?;
let harmonics = Harmonics::from_stor(
almanac.frame_from_uid(IAU_EARTH_FRAME)?,
HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true)?,
);
orbital_dyn.accel_models.push(harmonics);
let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
.with_guidance_law(ruggiero_ctrl.clone());
println!("{sc_dynamics}");
// Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.
// Let's start by defining the dispersion.
// 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.
// Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
let mc_rv = MultivariateNormal::new(
sc,
vec![StateDispersion::zero_mean(StateParameter::SMA, 3.0)],
)?;
let my_mc = MonteCarlo::new(
sc, // Nominal state
mc_rv,
"03_geo_sk".to_string(), // Scenario name
None, // No specific seed specified, so one will be drawn from the computer's entropy.
);
// Build the propagator setup.
let setup = Propagator::rk89(
sc_dynamics.clone(),
IntegratorOptions::builder()
.min_step(10.0_f64.seconds())
.error_ctrl(ErrorControl::RSSCartesianStep)
.build(),
);
let num_runs = 25;
let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);
assert_eq!(rslts.runs.len(), num_runs);
// For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
rslts.to_parquet(
"03_geo_sk.parquet",
Some(vec![
&EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
]),
ExportCfg::default(),
almanac,
)?;
Ok(())
}
More examples
examples/03_geo_analysis/raise.rs (line 136)
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's latest MetaAlmanac.
// This will automatically download the DE440s planetary ephemeris,
// the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
// parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
// planetary constants kernels.
// For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
// Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
// references to many functions.
let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
// Fetch the EME2000 frame from the Almabac
let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
// Define the orbit epoch
let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
// Build the spacecraft itself.
// Using slide 6 of https://aerospace.org/sites/default/files/2018-11/Davis-Mayberry_HPSEP_11212018.pdf
// for the "next gen" SEP characteristics.
// GTO start
let orbit = Orbit::keplerian(24505.9, 0.725, 7.05, 0.0, 0.0, 0.0, epoch, eme2k);
let sc = Spacecraft::builder()
.orbit(orbit)
.dry_mass_kg(1000.0) // 1000 kg of dry mass
.fuel_mass_kg(1000.0) // 1000 kg of fuel, totalling 2.0 tons
.srp(SrpConfig::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
.thruster(Thruster {
// "NEXT-STEP" row in Table 2
isp_s: 4435.0,
thrust_N: 0.472,
})
.mode(GuidanceMode::Thrust) // Start thrusting immediately.
.build();
let prop_time = 180.0 * Unit::Day;
// Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
let objectives = &[
Objective::within_tolerance(StateParameter::SMA, 42_165.0, 20.0),
Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5),
Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2),
];
// Ensure that we only thrust if we have more than 20% illumination.
let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2).unwrap();
println!("{ruggiero_ctrl}");
// Define the high fidelity dynamics
// Set up the spacecraft dynamics.
// Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
// The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
// We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
// We're using the JGM3 model here, which is the default in GMAT.
let mut jgm3_meta = MetaFile {
uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
};
// And let's download it if we don't have it yet.
jgm3_meta.process(true)?;
// Build the spherical harmonics.
// The harmonics must be computed in the body fixed frame.
// We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
let harmonics = Harmonics::from_stor(
almanac.frame_from_uid(IAU_EARTH_FRAME)?,
HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true).unwrap(),
);
// Include the spherical harmonics into the orbital dynamics.
orbital_dyn.accel_models.push(harmonics);
// We define the solar radiation pressure, using the default solar flux and accounting only
// for the eclipsing caused by the Earth.
let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
// Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
// acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
.with_guidance_law(ruggiero_ctrl.clone());
println!("{:x}", orbit);
// We specify a minimum step in the propagator because the Ruggiero control would otherwise drive this step very low.
let (final_state, traj) = Propagator::rk89(
sc_dynamics.clone(),
IntegratorOptions::builder()
.min_step(10.0_f64.seconds())
.error_ctrl(ErrorControl::RSSCartesianStep)
.build(),
)
.with(sc, almanac.clone())
.for_duration_with_traj(prop_time)?;
let fuel_usage = sc.fuel_mass_kg - final_state.fuel_mass_kg;
println!("{:x}", final_state.orbit);
println!("fuel usage: {:.3} kg", fuel_usage);
// Finally, export the results for analysis, including the penumbra percentage throughout the orbit raise.
traj.to_parquet(
"./03_geo_raise.parquet",
Some(vec![
&EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
]),
ExportCfg::default(),
almanac,
)?;
for status_line in ruggiero_ctrl.status(&final_state) {
println!("{status_line}");
}
ruggiero_ctrl
.achieved(&final_state)
.expect("objective not achieved");
Ok(())
}
examples/02_jwst_covar_monte_carlo/main.rs (line 151)
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's latest MetaAlmanac.
// For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
// Download the regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
// Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
let mut latest_jwst_ephem = MetaFile {
uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
crc32: None,
};
latest_jwst_ephem.process(true)?;
// Load this ephem in the general Almanac we're using for this analysis.
let almanac = Arc::new(
MetaAlmanac::latest()
.map_err(Box::new)?
.load_from_metafile(latest_jwst_ephem, true)?,
);
// By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
// in the BSP. We need this ID in order to query the ephemeris.
const JWST_NAIF_ID: i32 = -170;
// Let's build a frame in the J2000 orientation centered on the JWST.
const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
// Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
println!("JWST defined from {earliest_epoch} to {latest_epoch}");
// Fetch the state, printing it in the Earth J2000 frame.
let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
println!("{jwst_orbit:x}");
// Build the spacecraft
// SRP area assumed to be the full sunshield and mass if 6200.0 kg, c.f. https://webb.nasa.gov/content/about/faqs/facts.html
// SRP Coefficient of reflectivity assumed to be that of Kapton, i.e. 2 - 0.44 = 1.56, table 1 from https://amostech.com/TechnicalPapers/2018/Poster/Bengtson.pdf
let jwst = Spacecraft::builder()
.orbit(jwst_orbit)
.srp(SrpConfig {
area_m2: 21.197 * 14.162,
cr: 1.56,
})
.dry_mass_kg(6200.0)
.build();
// Build up the spacecraft uncertainty builder.
// We can use the spacecraft uncertainty structure to build this up.
// We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
// in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
// Nyx can also estimate the deviation of the spacecraft parameters.
let jwst_uncertainty = SpacecraftUncertainty::builder()
.nominal(jwst)
.frame(LocalFrame::RIC)
.x_km(0.5)
.y_km(0.3)
.z_km(1.5)
.vx_km_s(1e-4)
.vy_km_s(0.6e-3)
.vz_km_s(3e-3)
.build();
println!("{jwst_uncertainty}");
// Build the Kalman filter estimate.
// Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
// but this approach requires quite a bit more boilerplate code.
let jwst_estimate = jwst_uncertainty.to_estimate()?;
// Set up the spacecraft dynamics.
// We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
// We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
// Build the propagator set up to use for the whole analysis.
let setup = Propagator::default(dynamics);
// All of the analysis will use this duration.
let prediction_duration = 6.5 * Unit::Day;
// === Covariance mapping ===
// For the covariance mapping / prediction, we'll use the common orbit determination approach.
// This is done by setting up a spacecraft OD process, and predicting for the analysis duration.
let ckf = KF::no_snc(jwst_estimate);
// Build the propagation instance for the OD process.
let prop = setup.with(jwst.with_stm(), almanac.clone());
let mut odp = SpacecraftODProcess::ckf(prop, ckf, None, almanac.clone());
// Define the prediction step, i.e. how often we want to know the covariance.
let step = 1_i64.minutes();
// Finally, predict, and export the trajectory with covariance to a parquet file.
odp.predict_for(step, prediction_duration)?;
odp.to_parquet("./02_jwst_covar_map.parquet", ExportCfg::default())?;
// === Monte Carlo framework ===
// Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
let my_mc = MonteCarlo::new(
jwst, // Nominal state
jwst_estimate.to_random_variable()?,
"02_jwst".to_string(), // Scenario name
None, // No specific seed specified, so one will be drawn from the computer's entropy.
);
let num_runs = 5_000;
let rslts = my_mc.run_until_epoch(
setup,
almanac.clone(),
jwst.epoch() + prediction_duration,
num_runs,
);
assert_eq!(rslts.runs.len(), num_runs);
// Finally, export these results, computing the eclipse percentage for all of these results.
// For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
let umbra_event = eclipse_loc.to_umbra_event();
let penumbra_event = eclipse_loc.to_penumbra_event();
rslts.to_parquet(
"02_jwst_monte_carlo.parquet",
Some(vec![&umbra_event, &penumbra_event]),
ExportCfg::default(),
almanac,
)?;
Ok(())
}
examples/03_geo_analysis/drift.rs (line 155)
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's latest MetaAlmanac.
// This will automatically download the DE440s planetary ephemeris,
// the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
// parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
// planetary constants kernels.
// For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
// Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
// references to many functions.
let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
// Define the orbit epoch
let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
// Define the orbit.
// First we need to fetch the Earth J2000 from information from the Almanac.
// This allows the frame to include the gravitational parameters and the shape of the Earth,
// defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
// by loading a different set of planetary constants.
let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
// Placing this GEO bird just above Colorado.
// In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
// Print in in Keplerian form.
println!("{orbit:x}");
let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;
// Nyx is used for high fidelity propagation, not Keplerian propagation as above.
// Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
// models such as solar radiation pressure.
// Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
let sc = Spacecraft::builder()
.orbit(orbit)
.dry_mass_kg(9.60)
.srp(SrpConfig {
area_m2: 10e-4,
cr: 1.1,
})
.build();
println!("{sc:x}");
// Set up the spacecraft dynamics.
// Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
// The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
// We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
// We're using the JGM3 model here, which is the default in GMAT.
let mut jgm3_meta = MetaFile {
uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
};
// And let's download it if we don't have it yet.
jgm3_meta.process(true)?;
// Build the spherical harmonics.
// The harmonics must be computed in the body fixed frame.
// We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
let harmonics_21x21 = Harmonics::from_stor(
almanac.frame_from_uid(IAU_EARTH_FRAME)?,
HarmonicsMem::from_cof(&jgm3_meta.uri, 21, 21, true).unwrap(),
);
// Include the spherical harmonics into the orbital dynamics.
orbital_dyn.accel_models.push(harmonics_21x21);
// We define the solar radiation pressure, using the default solar flux and accounting only
// for the eclipsing caused by the Earth and Moon.
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
// acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
println!("{dynamics}");
// Finally, let's propagate this orbit to the same epoch as above.
// The first returned value is the spacecraft state at the final epoch.
// The second value is the full trajectory where the step size is variable step used by the propagator.
let (future_sc, trajectory) = Propagator::default(dynamics)
.with(sc, almanac.clone())
.until_epoch_with_traj(epoch + Unit::Century * 0.03)?;
println!("=== High fidelity propagation ===");
println!(
"SMA changed by {:.3} km",
orbit.sma_km()? - future_sc.orbit.sma_km()?
);
println!(
"ECC changed by {:.6}",
orbit.ecc()? - future_sc.orbit.ecc()?
);
println!(
"INC changed by {:.3e} deg",
orbit.inc_deg()? - future_sc.orbit.inc_deg()?
);
println!(
"RAAN changed by {:.3} deg",
orbit.raan_deg()? - future_sc.orbit.raan_deg()?
);
println!(
"AOP changed by {:.3} deg",
orbit.aop_deg()? - future_sc.orbit.aop_deg()?
);
println!(
"TA changed by {:.3} deg",
orbit.ta_deg()? - future_sc.orbit.ta_deg()?
);
// We also have access to the full trajectory throughout the propagation.
println!("{trajectory}");
println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");
// With the trajectory, let's build a few data products.
// 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.
let analysis_step = Unit::Minute * 5;
trajectory.to_parquet(
"./03_geo_hf_prop.parquet",
Some(vec![
&EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
]),
ExportCfg::builder().step(analysis_step).build(),
almanac.clone(),
)?;
// 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.
// We iterate over the trajectory, grabbing a state every two minutes.
let mut offset_s = vec![];
let mut epoch_str = vec![];
let mut longitude_deg = vec![];
let mut latitude_deg = vec![];
let mut altitude_km = vec![];
for state in trajectory.every(analysis_step) {
// Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
// These define the GEO stationkeeping box.
let this_epoch = state.epoch();
offset_s.push((this_epoch - orbit.epoch).to_seconds());
epoch_str.push(this_epoch.to_isoformat());
let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
longitude_deg.push(long_deg);
latitude_deg.push(lat_deg);
altitude_km.push(alt_km);
}
println!(
"Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
orig_long_deg - longitude_deg.last().unwrap()
);
println!(
"Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
orig_lat_deg - latitude_deg.last().unwrap()
);
println!(
"Altitude changed by {:.3} km -- Box is 30 km",
orig_alt_km - altitude_km.last().unwrap()
);
// Build the station keeping data frame.
let mut sk_df = df!(
"Offset (s)" => offset_s.clone(),
"Epoch (UTC)" => epoch_str.clone(),
"Longitude E-W (deg)" => longitude_deg,
"Latitude N-S (deg)" => latitude_deg,
"Altitude (km)" => altitude_km,
)?;
// Create a file to write the Parquet to
let file = File::create("./03_geo_lla.parquet").expect("Could not create file");
// Create a ParquetWriter and write the DataFrame to the file
ParquetWriter::new(file).finish(&mut sk_df)?;
Ok(())
}
Trait Implementations§
Source§impl Clone for EclipseLocator
impl Clone for EclipseLocator
Source§fn clone(&self) -> EclipseLocator
fn clone(&self) -> EclipseLocator
Returns a copy of the value. Read more
1.0.0 · Source§fn clone_from(&mut self, source: &Self)
fn clone_from(&mut self, source: &Self)
Performs copy-assignment from
source
. Read moreAuto Trait Implementations§
impl Freeze for EclipseLocator
impl RefUnwindSafe for EclipseLocator
impl Send for EclipseLocator
impl Sync for EclipseLocator
impl Unpin for EclipseLocator
impl UnwindSafe for EclipseLocator
Blanket Implementations§
Source§impl<T> BorrowMut<T> for Twhere
T: ?Sized,
impl<T> BorrowMut<T> for Twhere
T: ?Sized,
Source§fn borrow_mut(&mut self) -> &mut T
fn borrow_mut(&mut self) -> &mut T
Mutably borrows from an owned value. Read more
Source§impl<T> CloneToUninit for Twhere
T: Clone,
impl<T> CloneToUninit for Twhere
T: Clone,
§impl<T> Instrument for T
impl<T> Instrument for T
§fn instrument(self, span: Span) -> Instrumented<Self>
fn instrument(self, span: Span) -> Instrumented<Self>
§fn in_current_span(self) -> Instrumented<Self>
fn in_current_span(self) -> Instrumented<Self>
Source§impl<T> IntoEither for T
impl<T> IntoEither for T
Source§fn into_either(self, into_left: bool) -> Either<Self, Self>
fn into_either(self, into_left: bool) -> Either<Self, Self>
Converts
self
into a Left
variant of Either<Self, Self>
if into_left
is true
.
Converts self
into a Right
variant of Either<Self, Self>
otherwise. Read moreSource§fn into_either_with<F>(self, into_left: F) -> Either<Self, Self>
fn into_either_with<F>(self, into_left: F) -> Either<Self, Self>
Converts
self
into a Left
variant of Either<Self, Self>
if into_left(&self)
returns true
.
Converts self
into a Right
variant of Either<Self, Self>
otherwise. Read more§impl<T> Pointable for T
impl<T> Pointable for T
§impl<SS, SP> SupersetOf<SS> for SPwhere
SS: SubsetOf<SP>,
impl<SS, SP> SupersetOf<SS> for SPwhere
SS: SubsetOf<SP>,
§fn to_subset(&self) -> Option<SS>
fn to_subset(&self) -> Option<SS>
The inverse inclusion map: attempts to construct
self
from the equivalent element of its
superset. Read more§fn is_in_subset(&self) -> bool
fn is_in_subset(&self) -> bool
Checks if
self
is actually part of its subset T
(and can be converted to it).§fn to_subset_unchecked(&self) -> SS
fn to_subset_unchecked(&self) -> SS
Use with care! Same as
self.to_subset
but without any property checks. Always succeeds.§fn from_subset(element: &SS) -> SP
fn from_subset(element: &SS) -> SP
The inclusion map: converts
self
to the equivalent element of its superset.