Skip to content

Commit 692dd1c

Browse files
AndreyAPopovAndreyAPopov
committed
Cleaned up, added docs, and made the NRHO preset better.
1 parent a22b10b commit 692dd1c

File tree

6 files changed

+133
-10
lines changed

6 files changed

+133
-10
lines changed

‎docs/references.bib‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -209,3 +209,11 @@ @article{dSL98
209209
pages={83-100},
210210
year={1998}
211211
}
212+
213+
@phdthesis{Spr21,
214+
title={Trajectory design and targeting for applications to the exploration program in cislunar space},
215+
author={Spreen, Emily M Zimovan},
216+
year={2021},
217+
school={Purdue University}
218+
}
219+

‎toolbox/+otp/+cr3bp/+presets/Canonical.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
classdef Canonical < otp.cr3bp.CR3BPProblem
2+
% A trivial preset with a stable oscillating orbit around a lagrange point.
23
methods
34
function obj = Canonical(varargin)
45
% Create the Canonical CR3BP problem object.
5-
% This is a trivial steady state orbit oscillating around a lagrange point.
66

77
mu = 0.5;
88

‎toolbox/+otp/+cr3bp/+presets/NRHO.m‎

Lines changed: 51 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,31 @@
11
classdef NRHO < otp.cr3bp.CR3BPProblem
2-
%
2+
% This preset builds an $L_2$ halo orbit preset for the CR3BP based on
3+
% a table of reference initial conditions and periods. The table
4+
% contains $20$ entries, and thus there are $20$ different possible
5+
% orbits given by this preset. The table of $L_2$ halo orbits is given
6+
% as Table A.1 on page 218 of :cite:p:`Spr21`.
37

48
methods
5-
function obj = NRHO
9+
function obj = NRHO(varargin)
610
% Create the NRHO CR3BP problem object.
7-
% $L_2$ Halo orbit family member
11+
%
12+
% Parameters
13+
% ----------
14+
% Index : numeric(1, 1)
15+
% An integer index between 1 and 20.
16+
17+
T = getL2halotable();
18+
19+
p = inputParser();
20+
p.addParameter('Index', 1, @(x) mod(x, 1) == 0 & x > 0 & x <= size(T, 1));
21+
p.parse(varargin{:});
22+
results = p.Results;
23+
24+
pic = T(results.Index, :);
825

9-
ic = [1.0192741002 -0.1801324242 -0.0971927950];
26+
orbitalperiod = pic(1);
27+
ic = pic(2:end);
28+
1029
mE = otp.utils.PhysicalConstants.EarthMass;
1130
mL = otp.utils.PhysicalConstants.MoonMass;
1231
G = otp.utils.PhysicalConstants.GravitationalConstant;
@@ -17,9 +36,36 @@
1736
mu = muL/(muE + muL);
1837

1938
y0 = [ic(1); 0; ic(2); 0; ic(3); 0];
20-
tspan = [0, 10];
39+
tspan = [0, orbitalperiod];
2140
params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 0);
2241
obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params);
2342
end
2443
end
2544
end
45+
46+
47+
function T = getL2halotable()
48+
% create internal table of L_2 Halo orbits
49+
% see
50+
T = [ ...
51+
1.3632096570 1.0110350588 -0.1731500000 -0.0780141199; ...
52+
1.4748399512 1.0192741002 -0.1801324242 -0.0971927950; ...
53+
1.5872714606 1.0277926091 -0.1858044184 -0.1154896637; ...
54+
1.7008482705 1.0362652156 -0.1904417454 -0.1322667493; ...
55+
1.8155211042 1.0445681848 -0.1942338538 -0.1473971442; ...
56+
1.9311168544 1.0526805665 -0.1972878310 -0.1609628828; ...
57+
2.0474562565 1.0606277874 -0.1996480091 -0.1731020372; ...
58+
2.1741533495 1.0691059976 -0.2014140887 -0.1847950147; ...
59+
2.2915829886 1.0768767277 -0.2022559057 -0.1943508955; ...
60+
2.4093619266 1.0846726654 -0.2022295078 -0.2027817501; ...
61+
2.5273849254 1.0925906981 -0.2011567058 -0.2101017213; ...
62+
2.6455248145 1.1007585320 -0.1987609769 -0.2162644440; ...
63+
2.7635889805 1.1093498794 -0.1946155759 -0.2211327592; ...
64+
2.8909903824 1.1194130163 -0.1873686594 -0.2246002627; ...
65+
3.0073088423 1.1297344316 -0.1769810336 -0.2254855800; ...
66+
3.1205655022 1.1413664663 -0.1612996515 -0.2229158600; ...
67+
3.2266000495 1.1542349115 -0.1379744940 -0.2147411949; ...
68+
3.3173903769 1.1669663066 -0.1049833863 -0.1984458292; ...
69+
3.3833013605 1.1766385512 -0.0621463948 -0.1748356762; ...
70+
3.4154433338 1.1808881373 -0.0032736457 -0.1559184478];
71+
end

‎toolbox/+otp/+cr3bp/CR3BPParameters.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
% Relative mass of the system.
66
Mu %MATLAB ONLY: (1, 1) {otp.utils.validation.mustBeNumerical}
77

8-
% A factor to avoid singularity
8+
% A factor to avoid singularity when computing distances.
99
SoftFactor %MATLAB ONLY: (1, 1) {otp.utils.validation.mustBeNumerical}
1010
end
1111

‎toolbox/+otp/+cr3bp/CR3BPProblem.m‎

Lines changed: 58 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,41 @@
11
classdef CR3BPProblem < otp.Problem
2-
% The circular restricted three body problem
2+
% The circular restricted three body problem, following the formulation
3+
% presented in :cite:p:`Spr21`.
4+
%
5+
% The system represents the evolution of an object of negligent mass
6+
% near two large-mass objects that orbit each other a constant distance
7+
% apart. The ratio of the mass of the first object to the sum of the
8+
% total mass of the system is represented by the non-dimensional
9+
% constant $\mu$. The reference frame of the system is fixed to the
10+
% rotationg frame of the two objects, meaning that the objects have
11+
% fixed constant postions of $(\mu,0,0)^T$ for the first object, and
12+
% $(1 - \mu,0,0)^T$ for the second object. The evolution of the third
13+
% object of negligent mass is given by the following second-order
14+
% non-dimensionalized differential equation:
15+
%
16+
% $$
17+
% x'' &= σ(y - x),\\
18+
% y'' &= ρx - y - xz,\\
19+
% z'' &= xy - βz,\\
20+
% U &= \frac{1}{2} (x^2 + y^2) + \frac{1 - \mu}{d} + \frac{mu}{r},\\
21+
% d &= \sqrt{(x + \mu)^2 + y^2 + z^2},\\
22+
% r &= \sqrt{(x - 1 + \mu)^2 + y^2 + z^2},
23+
% $$
24+
%
25+
% where the system is converted to a differential equation in six
26+
% variables in the standard fashion. The distances $d$ and $r$ can
27+
% cause numerical instability as they approach zero, thus a softening
28+
% factor of $s^2$ is typically added under both of the square-roots.
29+
%
30+
% The system of equations is energy-preserving, meaning that the Jacobi
31+
% constant of the system,
32+
%
33+
% $$
34+
% J = 2U - x'^2 - y'^2 - z'^2,
35+
% $$
36+
%
37+
% is preserved throughout the evolution of the equations, though this
38+
% is typically not true numerically.
339
%
440
% Notes
541
% -----
@@ -14,7 +50,7 @@
1450
% Example
1551
% -------
1652
% >>> problem = otp.cr3bp.presets.NRHO;
17-
% >>> sol = model.solve('AbsTol', 1e-14, 'RelTol', 100*eps);
53+
% >>> sol = model.solve();
1854
% >>> problem.plotPhaseSpace(sol);
1955
%
2056
% See Also
@@ -37,11 +73,17 @@
3773
obj@otp.Problem('Circular Restricted 3 Body Problem', 6, timeSpan, y0, parameters);
3874
end
3975
end
76+
77+
properties (SetAccess = private)
78+
JacobiConstant
79+
end
4080

4181
methods (Access = protected)
4282
function onSettingsChanged(obj)
4383
mu = obj.Parameters.Mu;
4484
soft = obj.Parameters.SoftFactor;
85+
86+
obj.JacobiConstant = @(y) otp.cr3bp.jacobiconstant(y, mu, soft);
4587

4688
obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ...
4789
'Vectorized', 'on');
@@ -65,8 +107,21 @@ function onSettingsChanged(obj)
65107
end
66108

67109
function sol = internalSolve(obj, varargin)
68-
sol = internalSolve@otp.Problem(obj, 'Solver', otp.utils.Solver.Nonstiff, varargin{:});
110+
sol = internalSolve@otp.Problem(obj, ...
111+
'Solver', otp.utils.Solver.Nonstiff, ...
112+
'AbsTol', 1e-14, ...
113+
'RelTol', 100*eps, ...
114+
varargin{:});
115+
end
116+
117+
function fig = internalPlotPhaseSpace(obj, t, y, varargin)
118+
mu = obj.Parameters.Mu;
119+
fig = internalPlotPhaseSpace@otp.Problem(obj, t, y, ...
120+
'Vars', 1:3, varargin{:});
121+
hold on;
122+
scatter3(fig.Children(1), [mu, 1 - mu], [0, 0], [0, 0]);
69123
end
70124

71125
end
126+
72127
end

‎toolbox/+otp/+cr3bp/jacobiconstant.m‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
function J = jacobiconstant(y, mu, soft)
2+
3+
x = y(1:3, :);
4+
v = y(4:6, :);
5+
6+
d = sqrt((x(1, :) + mu).^2 + x(2, :)^2 + x(3, :)^2 + soft^2);
7+
r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + x(3, :)^2 + soft^2);
8+
9+
10+
U = 0.5*(x(1, :).^2 + x(2, :).^2) + (1 - mu)./d + mu./r;
11+
12+
J = 2*U - v(1, :).^2 - v(2, :).^2 - v(3, :).^2;
13+
14+
end

0 commit comments

Comments
 (0)