pub struct MonteCarlo<S: Interpolatable, Distr: Distribution<DispersedState<S>>>{
pub seed: Option<u128>,
pub random_state: Distr,
pub scenario: String,
pub nominal_state: S,
}
Expand description
A Monte Carlo framework, automatically running on all threads via a thread pool. This framework is targeted toward analysis of time-continuous variables. One caveat of the design is that the trajectory is used for post processing, not each individual state. This may prevent some event switching from being shown in GNC simulations.
Fields§
§seed: Option<u128>
Seed of the 64bit PCG random number generator
random_state: Distr
Generator of states for the Monte Carlo run
scenario: String
Name of this run, will be reflected in the progress bar and in the output structure
nominal_state: S
Implementations§
Source§impl<S: Interpolatable, Distr: Distribution<DispersedState<S>>> MonteCarlo<S, Distr>
impl<S: Interpolatable, Distr: Distribution<DispersedState<S>>> MonteCarlo<S, Distr>
Sourcepub fn new(
nominal_state: S,
random_variable: Distr,
scenario: String,
seed: Option<u128>,
) -> Self
pub fn new( nominal_state: S, random_variable: Distr, scenario: String, seed: Option<u128>, ) -> Self
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 // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
31
32 // Download the regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
33 // Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
34 let mut latest_jwst_ephem = MetaFile {
35 uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
36 crc32: None,
37 };
38 latest_jwst_ephem.process(true)?;
39
40 // Load this ephem in the general Almanac we're using for this analysis.
41 let almanac = Arc::new(
42 MetaAlmanac::latest()
43 .map_err(Box::new)?
44 .load_from_metafile(latest_jwst_ephem, true)?,
45 );
46
47 // By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
48 // in the BSP. We need this ID in order to query the ephemeris.
49 const JWST_NAIF_ID: i32 = -170;
50 // Let's build a frame in the J2000 orientation centered on the JWST.
51 const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
52
53 // Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
54 let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
55 println!("JWST defined from {earliest_epoch} to {latest_epoch}");
56 // Fetch the state, printing it in the Earth J2000 frame.
57 let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
58 println!("{jwst_orbit:x}");
59
60 // Build the spacecraft
61 // 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
62 // 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
63 let jwst = Spacecraft::builder()
64 .orbit(jwst_orbit)
65 .srp(SRPData {
66 area_m2: 21.197 * 14.162,
67 coeff_reflectivity: 1.56,
68 })
69 .mass(Mass::from_dry_mass(6200.0))
70 .build();
71
72 // Build up the spacecraft uncertainty builder.
73 // We can use the spacecraft uncertainty structure to build this up.
74 // We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
75 // in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
76 // Nyx can also estimate the deviation of the spacecraft parameters.
77 let jwst_uncertainty = SpacecraftUncertainty::builder()
78 .nominal(jwst)
79 .frame(LocalFrame::RIC)
80 .x_km(0.5)
81 .y_km(0.3)
82 .z_km(1.5)
83 .vx_km_s(1e-4)
84 .vy_km_s(0.6e-3)
85 .vz_km_s(3e-3)
86 .build();
87
88 println!("{jwst_uncertainty}");
89
90 // Build the Kalman filter estimate.
91 // Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
92 // but this approach requires quite a bit more boilerplate code.
93 let jwst_estimate = jwst_uncertainty.to_estimate()?;
94
95 // Set up the spacecraft dynamics.
96 // We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
97 // We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
98
99 let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
100 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102 // Finalize setting up the dynamics.
103 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
104
105 // Build the propagator set up to use for the whole analysis.
106 let setup = Propagator::default(dynamics);
107
108 // All of the analysis will use this duration.
109 let prediction_duration = 6.5 * Unit::Day;
110
111 // === Covariance mapping ===
112 // For the covariance mapping / prediction, we'll use the common orbit determination approach.
113 // This is done by setting up a spacecraft OD process, and predicting for the analysis duration.
114
115 let ckf = KF::no_snc(jwst_estimate);
116
117 // Build the propagation instance for the OD process.
118 let prop = setup.with(jwst.with_stm(), almanac.clone());
119 let mut odp = SpacecraftODProcess::ckf(prop, ckf, BTreeMap::new(), None, almanac.clone());
120
121 // Define the prediction step, i.e. how often we want to know the covariance.
122 let step = 1_i64.minutes();
123 // Finally, predict, and export the trajectory with covariance to a parquet file.
124 odp.predict_for(step, prediction_duration)?;
125 odp.to_parquet(
126 &TrackingDataArc::default(),
127 "./02_jwst_covar_map.parquet",
128 ExportCfg::default(),
129 )?;
130
131 // === Monte Carlo framework ===
132 // Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
133
134 let my_mc = MonteCarlo::new(
135 jwst, // Nominal state
136 jwst_estimate.to_random_variable()?,
137 "02_jwst".to_string(), // Scenario name
138 None, // No specific seed specified, so one will be drawn from the computer's entropy.
139 );
140
141 let num_runs = 5_000;
142 let rslts = my_mc.run_until_epoch(
143 setup,
144 almanac.clone(),
145 jwst.epoch() + prediction_duration,
146 num_runs,
147 );
148
149 assert_eq!(rslts.runs.len(), num_runs);
150 // Finally, export these results, computing the eclipse percentage for all of these results.
151
152 // For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
153 let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
154 let umbra_event = eclipse_loc.to_umbra_event();
155 let penumbra_event = eclipse_loc.to_penumbra_event();
156
157 rslts.to_parquet(
158 "02_jwst_monte_carlo.parquet",
159 Some(vec![&umbra_event, &penumbra_event]),
160 ExportCfg::default(),
161 almanac,
162 )?;
163
164 Ok(())
165}
Sourcepub fn run_until_nth_event<D, F>(
self,
prop: Propagator<D>,
almanac: Arc<Almanac>,
max_duration: Duration,
event: &F,
trigger: usize,
num_runs: usize,
) -> Results<S, PropResult<S>>where
D: Dynamics<StateType = S>,
F: EventEvaluator<S>,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
pub fn run_until_nth_event<D, F>(
self,
prop: Propagator<D>,
almanac: Arc<Almanac>,
max_duration: Duration,
event: &F,
trigger: usize,
num_runs: usize,
) -> Results<S, PropResult<S>>where
D: Dynamics<StateType = S>,
F: EventEvaluator<S>,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
Generate states and propagate each independently until a specific event is found trigger
times.
Sourcepub fn resume_run_until_nth_event<D, F>(
&self,
prop: Propagator<D>,
almanac: Arc<Almanac>,
skip: usize,
max_duration: Duration,
event: &F,
trigger: usize,
num_runs: usize,
) -> Results<S, PropResult<S>>where
D: Dynamics<StateType = S>,
F: EventEvaluator<S>,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
pub fn resume_run_until_nth_event<D, F>(
&self,
prop: Propagator<D>,
almanac: Arc<Almanac>,
skip: usize,
max_duration: Duration,
event: &F,
trigger: usize,
num_runs: usize,
) -> Results<S, PropResult<S>>where
D: Dynamics<StateType = S>,
F: EventEvaluator<S>,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
Generate states and propagate each independently until a specific event is found trigger
times.
Sourcepub fn run_until_epoch<D>(
self,
prop: Propagator<D>,
almanac: Arc<Almanac>,
end_epoch: Epoch,
num_runs: usize,
) -> Results<S, PropResult<S>>where
D: Dynamics<StateType = S>,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
pub fn run_until_epoch<D>(
self,
prop: Propagator<D>,
almanac: Arc<Almanac>,
end_epoch: Epoch,
num_runs: usize,
) -> Results<S, PropResult<S>>where
D: Dynamics<StateType = S>,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
Generate states and propagate each independently until a specific event is found trigger
times.
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 // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
31
32 // Download the regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
33 // Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
34 let mut latest_jwst_ephem = MetaFile {
35 uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
36 crc32: None,
37 };
38 latest_jwst_ephem.process(true)?;
39
40 // Load this ephem in the general Almanac we're using for this analysis.
41 let almanac = Arc::new(
42 MetaAlmanac::latest()
43 .map_err(Box::new)?
44 .load_from_metafile(latest_jwst_ephem, true)?,
45 );
46
47 // By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
48 // in the BSP. We need this ID in order to query the ephemeris.
49 const JWST_NAIF_ID: i32 = -170;
50 // Let's build a frame in the J2000 orientation centered on the JWST.
51 const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
52
53 // Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
54 let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
55 println!("JWST defined from {earliest_epoch} to {latest_epoch}");
56 // Fetch the state, printing it in the Earth J2000 frame.
57 let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
58 println!("{jwst_orbit:x}");
59
60 // Build the spacecraft
61 // 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
62 // 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
63 let jwst = Spacecraft::builder()
64 .orbit(jwst_orbit)
65 .srp(SRPData {
66 area_m2: 21.197 * 14.162,
67 coeff_reflectivity: 1.56,
68 })
69 .mass(Mass::from_dry_mass(6200.0))
70 .build();
71
72 // Build up the spacecraft uncertainty builder.
73 // We can use the spacecraft uncertainty structure to build this up.
74 // We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
75 // in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
76 // Nyx can also estimate the deviation of the spacecraft parameters.
77 let jwst_uncertainty = SpacecraftUncertainty::builder()
78 .nominal(jwst)
79 .frame(LocalFrame::RIC)
80 .x_km(0.5)
81 .y_km(0.3)
82 .z_km(1.5)
83 .vx_km_s(1e-4)
84 .vy_km_s(0.6e-3)
85 .vz_km_s(3e-3)
86 .build();
87
88 println!("{jwst_uncertainty}");
89
90 // Build the Kalman filter estimate.
91 // Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
92 // but this approach requires quite a bit more boilerplate code.
93 let jwst_estimate = jwst_uncertainty.to_estimate()?;
94
95 // Set up the spacecraft dynamics.
96 // We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
97 // We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
98
99 let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
100 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102 // Finalize setting up the dynamics.
103 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
104
105 // Build the propagator set up to use for the whole analysis.
106 let setup = Propagator::default(dynamics);
107
108 // All of the analysis will use this duration.
109 let prediction_duration = 6.5 * Unit::Day;
110
111 // === Covariance mapping ===
112 // For the covariance mapping / prediction, we'll use the common orbit determination approach.
113 // This is done by setting up a spacecraft OD process, and predicting for the analysis duration.
114
115 let ckf = KF::no_snc(jwst_estimate);
116
117 // Build the propagation instance for the OD process.
118 let prop = setup.with(jwst.with_stm(), almanac.clone());
119 let mut odp = SpacecraftODProcess::ckf(prop, ckf, BTreeMap::new(), None, almanac.clone());
120
121 // Define the prediction step, i.e. how often we want to know the covariance.
122 let step = 1_i64.minutes();
123 // Finally, predict, and export the trajectory with covariance to a parquet file.
124 odp.predict_for(step, prediction_duration)?;
125 odp.to_parquet(
126 &TrackingDataArc::default(),
127 "./02_jwst_covar_map.parquet",
128 ExportCfg::default(),
129 )?;
130
131 // === Monte Carlo framework ===
132 // Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
133
134 let my_mc = MonteCarlo::new(
135 jwst, // Nominal state
136 jwst_estimate.to_random_variable()?,
137 "02_jwst".to_string(), // Scenario name
138 None, // No specific seed specified, so one will be drawn from the computer's entropy.
139 );
140
141 let num_runs = 5_000;
142 let rslts = my_mc.run_until_epoch(
143 setup,
144 almanac.clone(),
145 jwst.epoch() + prediction_duration,
146 num_runs,
147 );
148
149 assert_eq!(rslts.runs.len(), num_runs);
150 // Finally, export these results, computing the eclipse percentage for all of these results.
151
152 // For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
153 let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
154 let umbra_event = eclipse_loc.to_umbra_event();
155 let penumbra_event = eclipse_loc.to_penumbra_event();
156
157 rslts.to_parquet(
158 "02_jwst_monte_carlo.parquet",
159 Some(vec![&umbra_event, &penumbra_event]),
160 ExportCfg::default(),
161 almanac,
162 )?;
163
164 Ok(())
165}
Sourcepub fn resume_run_until_epoch<D>(
&self,
prop: Propagator<D>,
almanac: Arc<Almanac>,
skip: usize,
end_epoch: Epoch,
num_runs: usize,
) -> Results<S, PropResult<S>>where
D: Dynamics<StateType = S>,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
pub fn resume_run_until_epoch<D>(
&self,
prop: Propagator<D>,
almanac: Arc<Almanac>,
skip: usize,
end_epoch: Epoch,
num_runs: usize,
) -> Results<S, PropResult<S>>where
D: Dynamics<StateType = S>,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
Resumes a Monte Carlo run by skipping the first skip
items, generating states only after that, and propagate each independently until the specified epoch.
Sourcepub fn generate_states(
&self,
skip: usize,
num_runs: usize,
seed: Option<u128>,
) -> Vec<(usize, DispersedState<S>)>
pub fn generate_states( &self, skip: usize, num_runs: usize, seed: Option<u128>, ) -> Vec<(usize, DispersedState<S>)>
Set up the seed and generate the states. This is useful for checking the generated states before running a large scale Monte Carlo.
Trait Implementations§
Source§impl<S: Interpolatable, Distr: Distribution<DispersedState<S>>> Display for MonteCarlo<S, Distr>
impl<S: Interpolatable, Distr: Distribution<DispersedState<S>>> Display for MonteCarlo<S, Distr>
Source§impl<S: Interpolatable, Distr: Distribution<DispersedState<S>>> LowerHex for MonteCarlo<S, Distr>
impl<S: Interpolatable, Distr: Distribution<DispersedState<S>>> LowerHex for MonteCarlo<S, Distr>
Auto Trait Implementations§
impl<S, Distr> Freeze for MonteCarlo<S, Distr>
impl<S, Distr> RefUnwindSafe for MonteCarlo<S, Distr>
impl<S, Distr> Send for MonteCarlo<S, Distr>
impl<S, Distr> Sync for MonteCarlo<S, Distr>
impl<S, Distr> Unpin for MonteCarlo<S, Distr>
impl<S, Distr> UnwindSafe for MonteCarlo<S, Distr>
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.