pub struct Results<S: Interpolatable, R>where
DefaultAllocator: Allocator<S::Size> + Allocator<S::Size, S::Size> + Allocator<S::VecLength>,
<DefaultAllocator as Allocator<S::VecLength>>::Buffer<f64>: Send,{
pub runs: Vec<Run<S, R>>,
pub scenario: String,
}
Expand description
A structure of Monte Carlo results
Fields§
§runs: Vec<Run<S, R>>
Raw data from each run, sorted by run index for O(1) access to each run
scenario: String
Name of this scenario
Implementations§
source§impl<S: Interpolatable> Results<S, PropResult<S>>
impl<S: Interpolatable> Results<S, PropResult<S>>
sourcepub fn every_value_of_between(
&self,
param: StateParameter,
step: Duration,
start: Epoch,
end: Epoch,
value_if_run_failed: Option<f64>,
) -> Vec<f64>
pub fn every_value_of_between( &self, param: StateParameter, step: Duration, start: Epoch, end: Epoch, value_if_run_failed: Option<f64>, ) -> Vec<f64>
Returns the value of the requested state parameter for all trajectories from start
to end
every step
and
using the value of value_if_run_failed
if set and skipping that run if the run failed
sourcepub fn every_value_of(
&self,
param: StateParameter,
step: Duration,
value_if_run_failed: Option<f64>,
) -> Vec<f64>
pub fn every_value_of( &self, param: StateParameter, step: Duration, value_if_run_failed: Option<f64>, ) -> Vec<f64>
Returns the value of the requested state parameter for all trajectories from the start to the end of each trajectory and
using the value of value_if_run_failed
if set and skipping that run if the run failed
sourcepub fn first_values_of(
&self,
param: StateParameter,
value_if_run_failed: Option<f64>,
) -> Vec<f64>
pub fn first_values_of( &self, param: StateParameter, value_if_run_failed: Option<f64>, ) -> Vec<f64>
Returns the value of the requested state parameter for the first state and
using the value of value_if_run_failed
if set and skipping that run if the run failed
sourcepub fn last_values_of(
&self,
param: StateParameter,
value_if_run_failed: Option<f64>,
) -> Vec<f64>
pub fn last_values_of( &self, param: StateParameter, value_if_run_failed: Option<f64>, ) -> Vec<f64>
Returns the value of the requested state parameter for the first state and
using the value of value_if_run_failed
if set and skipping that run if the run failed
sourcepub fn dispersion_values_of(
&self,
param: StateParameter,
) -> Result<Vec<f64>, MonteCarloError>
pub fn dispersion_values_of( &self, param: StateParameter, ) -> Result<Vec<f64>, MonteCarloError>
Returns the dispersion values of the requested state parameter
sourcepub fn to_parquet<P: AsRef<Path>>(
&self,
path: P,
events: Option<Vec<&dyn EventEvaluator<S>>>,
cfg: ExportCfg,
almanac: Arc<Almanac>,
) -> Result<PathBuf, Box<dyn Error>>
pub fn to_parquet<P: AsRef<Path>>( &self, path: P, events: Option<Vec<&dyn EventEvaluator<S>>>, cfg: ExportCfg, almanac: Arc<Almanac>, ) -> Result<PathBuf, Box<dyn Error>>
Examples found in repository?
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
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(())
}
Auto Trait Implementations§
impl<S, R> Freeze for Results<S, R>where
DefaultAllocator: Sized,
impl<S, R> !RefUnwindSafe for Results<S, R>
impl<S, R> Send for Results<S, R>
impl<S, R> Sync for Results<S, R>
impl<S, R> Unpin for Results<S, R>
impl<S, R> !UnwindSafe for Results<S, R>
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
§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>
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>
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>
self
from the equivalent element of its
superset. Read more§fn is_in_subset(&self) -> bool
fn is_in_subset(&self) -> bool
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
self.to_subset
but without any property checks. Always succeeds.§fn from_subset(element: &SS) -> SP
fn from_subset(element: &SS) -> SP
self
to the equivalent element of its superset.