From 98c3b8ce7741c614401b0da212486271f50577a5 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Fri, 26 Sep 2025 11:12:19 -1000 Subject: [PATCH 01/33] initial commit of CR3BP --- toolbox/+otp/+cr3bp/+presets/Canonical.m | 3 + toolbox/+otp/+cr3bp/+presets/NRHO.m | 25 ++++++++ toolbox/+otp/+cr3bp/CR3BPParameters.m | 25 ++++++++ toolbox/+otp/+cr3bp/CR3BPProblem.m | 72 ++++++++++++++++++++++++ toolbox/+otp/+cr3bp/f.m | 19 +++++++ 5 files changed, 144 insertions(+) create mode 100644 toolbox/+otp/+cr3bp/+presets/Canonical.m create mode 100644 toolbox/+otp/+cr3bp/+presets/NRHO.m create mode 100644 toolbox/+otp/+cr3bp/CR3BPParameters.m create mode 100644 toolbox/+otp/+cr3bp/CR3BPProblem.m create mode 100644 toolbox/+otp/+cr3bp/f.m diff --git a/toolbox/+otp/+cr3bp/+presets/Canonical.m b/toolbox/+otp/+cr3bp/+presets/Canonical.m new file mode 100644 index 00000000..2949ce7b --- /dev/null +++ b/toolbox/+otp/+cr3bp/+presets/Canonical.m @@ -0,0 +1,3 @@ +classdef Canonical < otp.cr3bp.CR3BPProblem + +end diff --git a/toolbox/+otp/+cr3bp/+presets/NRHO.m b/toolbox/+otp/+cr3bp/+presets/NRHO.m new file mode 100644 index 00000000..cf51bf8e --- /dev/null +++ b/toolbox/+otp/+cr3bp/+presets/NRHO.m @@ -0,0 +1,25 @@ +classdef NRHO < otp.cr3bp.CR3BPProblem + % + + methods + function obj = NRHO + % Create the NRHO CR3BP problem object. + % $L_2$ Halo orbit family member + + ic = [1.0192741002 -0.1801324242 -0.0971927950]; + mE = otp.utils.PhysicalConstants.EarthMass; + mL = otp.utils.PhysicalConstants.MoonMass; + G = otp.utils.PhysicalConstants.GravitationalConstant; + + % derive mu + muE = G*mE; + muL = G*mL; + mu = muL/(muE + muL); + + y0 = [ic(1); 0; ic(2); 0; ic(3); 0]; + tspan = [0, 10]; + params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 0); + obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); + end + end +end diff --git a/toolbox/+otp/+cr3bp/CR3BPParameters.m b/toolbox/+otp/+cr3bp/CR3BPParameters.m new file mode 100644 index 00000000..17540b74 --- /dev/null +++ b/toolbox/+otp/+cr3bp/CR3BPParameters.m @@ -0,0 +1,25 @@ +classdef CR3BPParameters < otp.Parameters + % Parameters for the CR3BP problem. + + properties + % Relative mass of the system. + Mu %MATLAB ONLY: (1, 1) {otp.utils.validation.mustBeNumerical} + + % A factor to avoid singularity + SoftFactor %MATLAB ONLY: (1, 1) {otp.utils.validation.mustBeNumerical} + end + + methods + function obj = CR3BPParameters(varargin) + % Create a CR3BP parameters object. + % + % Parameters + % ---------- + % varargin + % A variable number of name-value pairs. A name can be any property of this class, and the subsequent + % value initializes that property. + + obj = obj@otp.Parameters(varargin{:}); + end + end +end diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m new file mode 100644 index 00000000..ccc1362c --- /dev/null +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -0,0 +1,72 @@ +classdef CR3BPProblem < otp.Problem + % The circular restricted three body problem + % + % Notes + % ----- + % +---------------------+-----------------------------------------------+ + % | Type | ODE | + % +---------------------+-----------------------------------------------+ + % | Number of Variables | 6 | + % +---------------------+-----------------------------------------------+ + % | Stiff | no | + % +---------------------+-----------------------------------------------+ + % + % Example + % ------- + % >>> problem = otp.cr3bp.presets.Canonical; + % >>> sol = model.solve('AbsTol', 1e-14, 'RelTol', 100*eps); + % >>> problem.plotPhaseSpace(sol); + % + % See Also + % -------- + % otp.nbody.NBodyProblem + + methods + function obj = CR3BPProblem(timeSpan, y0, parameters) + % Create a CR3BP problem object. + % + % Parameters + % ---------- + % timeSpan : numeric(1, 2) + % The start and final time. + % y0 : numeric(6, 1) + % The initial conditions. + % parameters : otp.cr3bp.CR3BPParameters + % The parameters. + + obj@otp.Problem('Circular Restricted 3 Body Problem', 6, timeSpan, y0, parameters); + end + end + + methods (Access = protected) + function onSettingsChanged(obj) + mu = obj.Parameters.Mu; + soft = obj.Parameters.SoftFactor; + + obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ... + 'Vectorized', 'on'); + end + + function label = internalIndex2label(~, index) + switch index + case 1 + label = 'x'; + case 2 + label = 'y'; + case 3 + label = 'z'; + case 4 + label = 'dx'; + case 5 + label = 'dy'; + case 6 + label = 'dz'; + end + end + + function sol = internalSolve(obj, varargin) + sol = internalSolve@otp.Problem(obj, 'Solver', otp.utils.Solver.Nonstiff, varargin{:}); + end + + end +end diff --git a/toolbox/+otp/+cr3bp/f.m b/toolbox/+otp/+cr3bp/f.m new file mode 100644 index 00000000..14ad0427 --- /dev/null +++ b/toolbox/+otp/+cr3bp/f.m @@ -0,0 +1,19 @@ +function dy = f(~, y, mu, soft) + +x = y(1:3, :); +v = y(4:6, :); + +d = sqrt((x(1, :) + mu).^2 + x(2, :)^2 + x(3, :)^2 + soft^2); +r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + x(3, :)^2 + soft^2); +cd = (1 - mu)./d.^3; +cr = mu./r.^3; + +dx = v; +dv = [... + 2*v(2, :) + (1 - cd - cr).*x(1, :) - (cd + cr).*mu + cr; ... + -2*v(1, :) + (1 - cd - cr).*x(2, :); ... + -(cd + cr).*x(3, :)]; + +dy = [dx; dv]; + +end From a384befb442b659eb7f992b98b7826f16e990344 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Fri, 26 Sep 2025 11:16:11 -1000 Subject: [PATCH 02/33] Fixed the example --- toolbox/+otp/+cr3bp/CR3BPProblem.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index ccc1362c..23ad44bd 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -13,7 +13,7 @@ % % Example % ------- - % >>> problem = otp.cr3bp.presets.Canonical; + % >>> problem = otp.cr3bp.presets.NRHO; % >>> sol = model.solve('AbsTol', 1e-14, 'RelTol', 100*eps); % >>> problem.plotPhaseSpace(sol); % From 72108dfdee3bc0e8f61113d52a5a3a89d74d30ab Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Fri, 26 Sep 2025 11:23:38 -1000 Subject: [PATCH 03/33] Added trivial orbit as canonical --- toolbox/+otp/+cr3bp/+presets/Canonical.m | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/toolbox/+otp/+cr3bp/+presets/Canonical.m b/toolbox/+otp/+cr3bp/+presets/Canonical.m index 2949ce7b..f2f04118 100644 --- a/toolbox/+otp/+cr3bp/+presets/Canonical.m +++ b/toolbox/+otp/+cr3bp/+presets/Canonical.m @@ -1,3 +1,15 @@ classdef Canonical < otp.cr3bp.CR3BPProblem - + methods + function obj = Canonical(varargin) + % Create the NRHO CR3BP problem object. + % This is a trivial steady state orbit. + + mu = 0.5; + + y0 = zeros(6, 1); + tspan = [0, 10]; + params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 1e-3, varargin{:}); + obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); + end + end end From a22b10b31bfcde039dadff8eb390e26e41ce02ab Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Fri, 26 Sep 2025 12:27:30 -1000 Subject: [PATCH 04/33] Updated Canonical --- toolbox/+otp/+cr3bp/+presets/Canonical.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/toolbox/+otp/+cr3bp/+presets/Canonical.m b/toolbox/+otp/+cr3bp/+presets/Canonical.m index f2f04118..22f31f54 100644 --- a/toolbox/+otp/+cr3bp/+presets/Canonical.m +++ b/toolbox/+otp/+cr3bp/+presets/Canonical.m @@ -1,12 +1,12 @@ classdef Canonical < otp.cr3bp.CR3BPProblem methods function obj = Canonical(varargin) - % Create the NRHO CR3BP problem object. - % This is a trivial steady state orbit. + % Create the Canonical CR3BP problem object. + % This is a trivial steady state orbit oscillating around a lagrange point. mu = 0.5; - y0 = zeros(6, 1); + y0 = [0; 0; 0; 0; 0; 1]; tspan = [0, 10]; params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 1e-3, varargin{:}); obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); From 692dd1cf6ad7e343006b414c80e5c7bbd5696a2b Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Sat, 27 Sep 2025 09:09:33 -1000 Subject: [PATCH 05/33] Cleaned up, added docs, and made the NRHO preset better. --- docs/references.bib | 8 ++++ toolbox/+otp/+cr3bp/+presets/Canonical.m | 2 +- toolbox/+otp/+cr3bp/+presets/NRHO.m | 56 ++++++++++++++++++++-- toolbox/+otp/+cr3bp/CR3BPParameters.m | 2 +- toolbox/+otp/+cr3bp/CR3BPProblem.m | 61 ++++++++++++++++++++++-- toolbox/+otp/+cr3bp/jacobiconstant.m | 14 ++++++ 6 files changed, 133 insertions(+), 10 deletions(-) create mode 100644 toolbox/+otp/+cr3bp/jacobiconstant.m diff --git a/docs/references.bib b/docs/references.bib index 15a4c373..f460e9a0 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -209,3 +209,11 @@ @article{dSL98 pages={83-100}, year={1998} } + +@phdthesis{Spr21, + title={Trajectory design and targeting for applications to the exploration program in cislunar space}, + author={Spreen, Emily M Zimovan}, + year={2021}, + school={Purdue University} +} + diff --git a/toolbox/+otp/+cr3bp/+presets/Canonical.m b/toolbox/+otp/+cr3bp/+presets/Canonical.m index 22f31f54..b374efc7 100644 --- a/toolbox/+otp/+cr3bp/+presets/Canonical.m +++ b/toolbox/+otp/+cr3bp/+presets/Canonical.m @@ -1,8 +1,8 @@ classdef Canonical < otp.cr3bp.CR3BPProblem + % A trivial preset with a stable oscillating orbit around a lagrange point. methods function obj = Canonical(varargin) % Create the Canonical CR3BP problem object. - % This is a trivial steady state orbit oscillating around a lagrange point. mu = 0.5; diff --git a/toolbox/+otp/+cr3bp/+presets/NRHO.m b/toolbox/+otp/+cr3bp/+presets/NRHO.m index cf51bf8e..ac31da8f 100644 --- a/toolbox/+otp/+cr3bp/+presets/NRHO.m +++ b/toolbox/+otp/+cr3bp/+presets/NRHO.m @@ -1,12 +1,31 @@ classdef NRHO < otp.cr3bp.CR3BPProblem - % + % This preset builds an $L_2$ halo orbit preset for the CR3BP based on + % a table of reference initial conditions and periods. The table + % contains $20$ entries, and thus there are $20$ different possible + % orbits given by this preset. The table of $L_2$ halo orbits is given + % as Table A.1 on page 218 of :cite:p:`Spr21`. methods - function obj = NRHO + function obj = NRHO(varargin) % Create the NRHO CR3BP problem object. - % $L_2$ Halo orbit family member + % + % Parameters + % ---------- + % Index : numeric(1, 1) + % An integer index between 1 and 20. + + T = getL2halotable(); + + p = inputParser(); + p.addParameter('Index', 1, @(x) mod(x, 1) == 0 & x > 0 & x <= size(T, 1)); + p.parse(varargin{:}); + results = p.Results; + + pic = T(results.Index, :); - ic = [1.0192741002 -0.1801324242 -0.0971927950]; + orbitalperiod = pic(1); + ic = pic(2:end); + mE = otp.utils.PhysicalConstants.EarthMass; mL = otp.utils.PhysicalConstants.MoonMass; G = otp.utils.PhysicalConstants.GravitationalConstant; @@ -17,9 +36,36 @@ mu = muL/(muE + muL); y0 = [ic(1); 0; ic(2); 0; ic(3); 0]; - tspan = [0, 10]; + tspan = [0, orbitalperiod]; params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 0); obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); end end end + + +function T = getL2halotable() +% create internal table of L_2 Halo orbits +% see +T = [ ... + 1.3632096570 1.0110350588 -0.1731500000 -0.0780141199; ... + 1.4748399512 1.0192741002 -0.1801324242 -0.0971927950; ... + 1.5872714606 1.0277926091 -0.1858044184 -0.1154896637; ... + 1.7008482705 1.0362652156 -0.1904417454 -0.1322667493; ... + 1.8155211042 1.0445681848 -0.1942338538 -0.1473971442; ... + 1.9311168544 1.0526805665 -0.1972878310 -0.1609628828; ... + 2.0474562565 1.0606277874 -0.1996480091 -0.1731020372; ... + 2.1741533495 1.0691059976 -0.2014140887 -0.1847950147; ... + 2.2915829886 1.0768767277 -0.2022559057 -0.1943508955; ... + 2.4093619266 1.0846726654 -0.2022295078 -0.2027817501; ... + 2.5273849254 1.0925906981 -0.2011567058 -0.2101017213; ... + 2.6455248145 1.1007585320 -0.1987609769 -0.2162644440; ... + 2.7635889805 1.1093498794 -0.1946155759 -0.2211327592; ... + 2.8909903824 1.1194130163 -0.1873686594 -0.2246002627; ... + 3.0073088423 1.1297344316 -0.1769810336 -0.2254855800; ... + 3.1205655022 1.1413664663 -0.1612996515 -0.2229158600; ... + 3.2266000495 1.1542349115 -0.1379744940 -0.2147411949; ... + 3.3173903769 1.1669663066 -0.1049833863 -0.1984458292; ... + 3.3833013605 1.1766385512 -0.0621463948 -0.1748356762; ... + 3.4154433338 1.1808881373 -0.0032736457 -0.1559184478]; +end diff --git a/toolbox/+otp/+cr3bp/CR3BPParameters.m b/toolbox/+otp/+cr3bp/CR3BPParameters.m index 17540b74..77232e81 100644 --- a/toolbox/+otp/+cr3bp/CR3BPParameters.m +++ b/toolbox/+otp/+cr3bp/CR3BPParameters.m @@ -5,7 +5,7 @@ % Relative mass of the system. Mu %MATLAB ONLY: (1, 1) {otp.utils.validation.mustBeNumerical} - % A factor to avoid singularity + % A factor to avoid singularity when computing distances. SoftFactor %MATLAB ONLY: (1, 1) {otp.utils.validation.mustBeNumerical} end diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index 23ad44bd..e5fa1fd2 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -1,5 +1,41 @@ classdef CR3BPProblem < otp.Problem - % The circular restricted three body problem + % The circular restricted three body problem, following the formulation + % presented in :cite:p:`Spr21`. + % + % The system represents the evolution of an object of negligent mass + % near two large-mass objects that orbit each other a constant distance + % apart. The ratio of the mass of the first object to the sum of the + % total mass of the system is represented by the non-dimensional + % constant $\mu$. The reference frame of the system is fixed to the + % rotationg frame of the two objects, meaning that the objects have + % fixed constant postions of $(\mu,0,0)^T$ for the first object, and + % $(1 - \mu,0,0)^T$ for the second object. The evolution of the third + % object of negligent mass is given by the following second-order + % non-dimensionalized differential equation: + % + % $$ + % x'' &= σ(y - x),\\ + % y'' &= ρx - y - xz,\\ + % z'' &= xy - βz,\\ + % U &= \frac{1}{2} (x^2 + y^2) + \frac{1 - \mu}{d} + \frac{mu}{r},\\ + % d &= \sqrt{(x + \mu)^2 + y^2 + z^2},\\ + % r &= \sqrt{(x - 1 + \mu)^2 + y^2 + z^2}, + % $$ + % + % where the system is converted to a differential equation in six + % variables in the standard fashion. The distances $d$ and $r$ can + % cause numerical instability as they approach zero, thus a softening + % factor of $s^2$ is typically added under both of the square-roots. + % + % The system of equations is energy-preserving, meaning that the Jacobi + % constant of the system, + % + % $$ + % J = 2U - x'^2 - y'^2 - z'^2, + % $$ + % + % is preserved throughout the evolution of the equations, though this + % is typically not true numerically. % % Notes % ----- @@ -14,7 +50,7 @@ % Example % ------- % >>> problem = otp.cr3bp.presets.NRHO; - % >>> sol = model.solve('AbsTol', 1e-14, 'RelTol', 100*eps); + % >>> sol = model.solve(); % >>> problem.plotPhaseSpace(sol); % % See Also @@ -37,11 +73,17 @@ obj@otp.Problem('Circular Restricted 3 Body Problem', 6, timeSpan, y0, parameters); end end + + properties (SetAccess = private) + JacobiConstant + end methods (Access = protected) function onSettingsChanged(obj) mu = obj.Parameters.Mu; soft = obj.Parameters.SoftFactor; + + obj.JacobiConstant = @(y) otp.cr3bp.jacobiconstant(y, mu, soft); obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ... 'Vectorized', 'on'); @@ -65,8 +107,21 @@ function onSettingsChanged(obj) end function sol = internalSolve(obj, varargin) - sol = internalSolve@otp.Problem(obj, 'Solver', otp.utils.Solver.Nonstiff, varargin{:}); + sol = internalSolve@otp.Problem(obj, ... + 'Solver', otp.utils.Solver.Nonstiff, ... + 'AbsTol', 1e-14, ... + 'RelTol', 100*eps, ... + varargin{:}); + end + + function fig = internalPlotPhaseSpace(obj, t, y, varargin) + mu = obj.Parameters.Mu; + fig = internalPlotPhaseSpace@otp.Problem(obj, t, y, ... + 'Vars', 1:3, varargin{:}); + hold on; + scatter3(fig.Children(1), [mu, 1 - mu], [0, 0], [0, 0]); end end + end diff --git a/toolbox/+otp/+cr3bp/jacobiconstant.m b/toolbox/+otp/+cr3bp/jacobiconstant.m new file mode 100644 index 00000000..066496b5 --- /dev/null +++ b/toolbox/+otp/+cr3bp/jacobiconstant.m @@ -0,0 +1,14 @@ +function J = jacobiconstant(y, mu, soft) + +x = y(1:3, :); +v = y(4:6, :); + +d = sqrt((x(1, :) + mu).^2 + x(2, :)^2 + x(3, :)^2 + soft^2); +r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + x(3, :)^2 + soft^2); + + +U = 0.5*(x(1, :).^2 + x(2, :).^2) + (1 - mu)./d + mu./r; + +J = 2*U - v(1, :).^2 - v(2, :).^2 - v(3, :).^2; + +end From 1cfb3683aa63e193520f3ab2c2f2c346e9f908bb Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Sat, 27 Sep 2025 09:12:45 -1000 Subject: [PATCH 06/33] Renamed NRHO preset --- toolbox/+otp/+cr3bp/+presets/{NRHO.m => L2NRHO.m} | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) rename toolbox/+otp/+cr3bp/+presets/{NRHO.m => L2NRHO.m} (97%) diff --git a/toolbox/+otp/+cr3bp/+presets/NRHO.m b/toolbox/+otp/+cr3bp/+presets/L2NRHO.m similarity index 97% rename from toolbox/+otp/+cr3bp/+presets/NRHO.m rename to toolbox/+otp/+cr3bp/+presets/L2NRHO.m index ac31da8f..44fa3935 100644 --- a/toolbox/+otp/+cr3bp/+presets/NRHO.m +++ b/toolbox/+otp/+cr3bp/+presets/L2NRHO.m @@ -1,4 +1,4 @@ -classdef NRHO < otp.cr3bp.CR3BPProblem +classdef L2NRHO < otp.cr3bp.CR3BPProblem % This preset builds an $L_2$ halo orbit preset for the CR3BP based on % a table of reference initial conditions and periods. The table % contains $20$ entries, and thus there are $20$ different possible @@ -6,7 +6,7 @@ % as Table A.1 on page 218 of :cite:p:`Spr21`. methods - function obj = NRHO(varargin) + function obj = L2NRHO(varargin) % Create the NRHO CR3BP problem object. % % Parameters From fe52ebd7a85729587cac5c61fa6d4d474a060492 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Sat, 27 Sep 2025 09:14:27 -1000 Subject: [PATCH 07/33] fixed a typo --- toolbox/+otp/+cr3bp/CR3BPProblem.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index e5fa1fd2..caead95b 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -8,7 +8,7 @@ % total mass of the system is represented by the non-dimensional % constant $\mu$. The reference frame of the system is fixed to the % rotationg frame of the two objects, meaning that the objects have - % fixed constant postions of $(\mu,0,0)^T$ for the first object, and + % fixed constant positions of $(\mu,0,0)^T$ for the first object, and % $(1 - \mu,0,0)^T$ for the second object. The evolution of the third % object of negligent mass is given by the following second-order % non-dimensionalized differential equation: From e05cacc5abb0bf8f60a78390d0a3748b98de9bc5 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Sat, 27 Sep 2025 09:30:27 -1000 Subject: [PATCH 08/33] fix for octave --- toolbox/+otp/+cr3bp/CR3BPProblem.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index caead95b..b81577e5 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -119,7 +119,7 @@ function onSettingsChanged(obj) fig = internalPlotPhaseSpace@otp.Problem(obj, t, y, ... 'Vars', 1:3, varargin{:}); hold on; - scatter3(fig.Children(1), [mu, 1 - mu], [0, 0], [0, 0]); + scatter3([mu, 1 - mu], [0, 0], [0, 0]); end end From c6bca5c245add3c8988d3d5af4d66fd29dee5f61 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Sat, 27 Sep 2025 13:09:11 -1000 Subject: [PATCH 09/33] Added more Halo orbit presets --- .gitignore | 1 + toolbox/+otp/+cr3bp/+presets/HaloOrbit.m | 180 +++++++++++++++++++++++ toolbox/+otp/+cr3bp/+presets/L2NRHO.m | 71 --------- 3 files changed, 181 insertions(+), 71 deletions(-) create mode 100644 toolbox/+otp/+cr3bp/+presets/HaloOrbit.m delete mode 100644 toolbox/+otp/+cr3bp/+presets/L2NRHO.m diff --git a/.gitignore b/.gitignore index 18c6c30a..0a083382 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,4 @@ docs/build/ docs/problems/ .DS_Store +tests/octave-workspace diff --git a/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m b/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m new file mode 100644 index 00000000..d04f870f --- /dev/null +++ b/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m @@ -0,0 +1,180 @@ +classdef HaloOrbit < otp.cr3bp.CR3BPProblem + % This preset builds a halo orbit preset for the CR3BP based on + % a table of reference initial conditions and periods. There are four + % possible different orbits given by this preset, $L_2$, $P2HO_1$, + % $P2HO_2$, $P4HO_1$, and $P4HO_2$. The tables contain $20$ entries, and + % thus there are $100$ different possible orbits given by this preset. + % The tables of $L_2$ halo orbits are given as Table A.1--A.5 on pages + % 218--222 of :cite:p:`Spr21`. + + methods + function obj = HaloOrbit(varargin) + % Create the NRHO CR3BP problem object. + % + % Parameters + % ---------- + % OrbitType : string + % One of the following: + % a. 'L2' for the L2 Orbit family (default) + % b. 'P2HO1' for the period doubling halo orbit of type 1 + % c. 'P2HO2' for the period doubling halo orbit of type 2 + % d. 'P4HO1' for the period quadrupling halo orbit of type 1 + % e. 'P4HO2' for the period quadrupling halo orbit of type 2 + % Index : numeric(1, 1) + % An integer index between 1 and 20. + + p = inputParser(); + p.addParameter('OrbitType', 'L2', ... + @(x) strcmp(x, 'L2') || strcmp(x, 'P2HO1') || strcmp(x, 'P2HO2') ... + || strcmp(x, 'P4HO1') || strcmp(x, 'P4HO2')); + p.addParameter('Index', 1, @(x) mod(x, 1) == 0 && x > 0 && x <= 20); + p.parse(varargin{:}); + results = p.Results; + + T = gethalotable(results.OrbitType); + + pic = T(results.Index, :); + + orbitalperiod = pic(1); + ic = pic(2:end); + + mE = otp.utils.PhysicalConstants.EarthMass; + mL = otp.utils.PhysicalConstants.MoonMass; + G = otp.utils.PhysicalConstants.GravitationalConstant; + + % derive mu + muE = G*mE; + muL = G*mL; + mu = muL/(muE + muL); + + y0 = [ic(1); 0; ic(2); 0; ic(3); 0]; + tspan = [0, orbitalperiod]; + params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 0); + obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); + end + end +end + + +function T = gethalotable(orbittype) +% create internal table of halo orbit initial conditions +% The first entry of each row is the orbital period, the second entry is +% the $x$ initial condition, the third entry is the $z$ initial condition +% and the fourth entry is the $y'$ intial condition. +switch orbittype + case 'L2' + T = [ ... + 1.3632096570 1.0110350588 -0.1731500000 -0.0780141199; ... + 1.4748399512 1.0192741002 -0.1801324242 -0.0971927950; ... + 1.5872714606 1.0277926091 -0.1858044184 -0.1154896637; ... + 1.7008482705 1.0362652156 -0.1904417454 -0.1322667493; ... + 1.8155211042 1.0445681848 -0.1942338538 -0.1473971442; ... + 1.9311168544 1.0526805665 -0.1972878310 -0.1609628828; ... + 2.0474562565 1.0606277874 -0.1996480091 -0.1731020372; ... + 2.1741533495 1.0691059976 -0.2014140887 -0.1847950147; ... + 2.2915829886 1.0768767277 -0.2022559057 -0.1943508955; ... + 2.4093619266 1.0846726654 -0.2022295078 -0.2027817501; ... + 2.5273849254 1.0925906981 -0.2011567058 -0.2101017213; ... + 2.6455248145 1.1007585320 -0.1987609769 -0.2162644440; ... + 2.7635889805 1.1093498794 -0.1946155759 -0.2211327592; ... + 2.8909903824 1.1194130163 -0.1873686594 -0.2246002627; ... + 3.0073088423 1.1297344316 -0.1769810336 -0.2254855800; ... + 3.1205655022 1.1413664663 -0.1612996515 -0.2229158600; ... + 3.2266000495 1.1542349115 -0.1379744940 -0.2147411949; ... + 3.3173903769 1.1669663066 -0.1049833863 -0.1984458292; ... + 3.3833013605 1.1766385512 -0.0621463948 -0.1748356762; ... + 3.4154433338 1.1808881373 -0.0032736457 -0.1559184478]; + case 'P2HO1' + T = [ ... + 2.7486723000 1.0124213729 -0.1739463747 -0.0801103670; ... + 2.8235719903 0.9686100061 -0.1684646845 -0.0555868296; ... + 2.9470915455 0.9428902766 -0.1621459176 -0.0367463465; ... + 3.0718373938 0.9280339523 -0.1570215893 -0.0278535009; ... + 3.2098634414 0.9177532010 -0.1525284780 -0.0253346758; ... + 3.3390119532 0.9116052635 -0.1492702168 -0.0275809790; ... + 3.4781621508 0.9074818076 -0.1466452811 -0.0334316688; ... + 3.6071172693 0.9053614223 -0.1449312614 -0.0411597104; ... + 3.7456075513 0.9044951611 -0.1437576475 -0.0512886367; ... + 3.8738428281 0.9047472035 -0.1431979716 -0.0619366228; ... + 4.0115944839 0.9059473646 -0.1430690015 -0.0743674595; ... + 4.1392474276 0.9077780758 -0.1433117179 -0.0865211750; ... + 4.2765261366 0.9103912789 -0.1438936424 -0.0999906823; ... + 4.4039019522 0.9133170339 -0.1446839313 -0.1126257931; ... + 4.5410586110 0.9169169198 -0.1457671148 -0.1261483024; ... + 4.6684708483 0.9206100089 -0.1469654829 -0.1384542728; ... + 4.8057987675 0.9249032213 -0.1484447194 -0.1512842987; ... + 4.9334592177 0.9291366414 -0.1499785567 -0.1627035911; ... + 5.0711130707 0.9339173119 -0.1517856322 -0.1743951579; ... + 5.2089533970 0.9388868141 -0.1537336561 -0.1854168651]; + case 'P2HO2' + T = [ ... + 4.5243776465 0.9902279661 0.0392479429 0.7431504374; ... + 4.3364693699 0.9913034068 0.0386425658 0.7476198379; ... + 4.3344588357 0.9925765901 0.0404028452 0.7275812069; ... + 4.3322427233 0.9936999422 0.0418211289 0.7119297395; ... + 4.3274956856 0.9956999819 0.0440891033 0.6876036122; ... + 4.3154119806 0.9997963248 0.0479397591 0.6474122281; ... + 4.2966706889 1.0050420842 0.0517347036 0.6074272113; ... + 4.2610464740 1.0135348381 0.0559827198 0.5575089518; ... + 4.2236953788 1.0214960279 0.0583763748 0.5199524066; ... + 4.2114971687 1.0239852380 0.0588603231 0.5093886860; ... + 4.2066140086 1.0249716121 0.0590196996 0.5053259810; ... + 4.2037569181 1.0255464522 0.0591042680 0.5029882761; ... + 4.2014068389 1.0260181267 0.0591691271 0.5010860112; ... + 4.1989300286 1.0265141820 0.0592329594 0.4991004650; ... + 4.1953346327 1.0272324868 0.0593174846 0.4962518855; ... + 4.1879862293 1.0286948426 0.0594609346 0.4905454113; ... + 4.1645925208 1.0333161410 0.0596666029 0.4732448565; ... + 4.0891753953 1.0483318888 0.0577756525 0.4227218693; ... + 3.9680003798 1.0758753256 0.0417283774 0.3430241217; ... + 3.9102015266 1.0926153718 0.0005492180 0.2990655449]; + case 'P4HO1' + T = [ ... + 8.2676291788 1.0619393673 -0.1999741894 -0.1749961447; ... + 8.3305139834 1.0681129545 -0.1976100673 -0.1868722998; ... + 8.3392497293 1.0738246899 -0.1946216607 -0.1984738200; ... + 8.3519917053 1.0786654098 -0.1919984087 -0.2083581121; ... + 8.3693678945 1.0833295002 -0.1893927253 -0.2179310067; ... + 8.3878867600 1.0871892065 -0.1871812140 -0.2258888455; ... + 8.4096026243 1.0909131908 -0.1850039744 -0.2335944979; ... + 8.4309304649 1.0940376079 -0.1831478491 -0.2400766927; ... + 8.4547978250 1.0971075676 -0.1813014490 -0.2464568317; ... + 8.4775801091 1.0997316965 -0.1797082140 -0.2519152440; ... + 8.5026193458 1.1023550439 -0.1781043165 -0.2573726538; ... + 8.5262347341 1.1046324969 -0.1767049403 -0.2621078568; ... + 8.5519768095 1.1069405119 -0.1752821493 -0.2669009993; ... + 8.5761129683 1.1089683638 -0.1740297442 -0.2711051374; ... + 8.6023105873 1.1110450669 -0.1727463692 -0.2754011052; ... + 8.6267952587 1.1128866230 -0.1716087730 -0.2792006165; ... + 8.6533063618 1.1147879601 -0.1704357869 -0.2831115888; ... + 8.6780368548 1.1164862743 -0.1693902247 -0.2865931191; ... + 8.7047739081 1.1182510772 -0.1683067141 -0.2901975105; ... + 8.7316043012 1.1199566179 -0.1672632043 -0.2936663490]; + case 'P4HO2' + T = [ ... + 9.0487552928 1.0749349573 -0.2021215571 -0.1920690379; ... + 9.0396341760 1.0578442875 -0.2047569352 -0.1678532831; ... + 9.0542791209 1.0524484578 -0.2047233749 -0.1606194256; ... + 9.0601028655 1.0505487899 -0.2045610040 -0.1581773718; ... + 9.0646190673 1.0491368344 -0.2043821958 -0.1564064537; ... + 9.0707205473 1.0472961590 -0.2040685317 -0.1541623227; ... + 9.0786237393 1.0450028327 -0.2035390634 -0.1514838937; ... + 9.0888613168 1.0421535217 -0.2026443039 -0.1483728351; ... + 9.1037491359 1.0382266494 -0.2009270031 -0.1445796433; ... + 9.1306372229 1.0319435819 -0.1968531223 -0.1401578788; ... + 9.2021506431 1.0215863525 -0.1864218292 -0.1399647167; ... + 9.3830140082 1.0126738542 -0.1723454500 -0.1560887631; ... + 9.5440425474 1.0113733655 -0.1648460585 -0.1739113299; ... + 9.6990195652 1.0129644895 -0.1593270671 -0.1918297866; ... + 10.2961547121 1.0286792032 -0.1457428883 -0.2522546927; ... + 11.7792576999 1.0624281188 -0.1343562748 -0.3275531781; ... + 13.7831848966 1.0850769369 -0.1379181777 -0.3653699405; ... + 16.0024825055 1.0955445136 -0.1591111395 -0.3767847569; ... + 18.1058513583 1.0954718438 -0.1988427540 -0.3675757435; ... + 20.0406128078 1.0829659022 -0.2591580427 -0.3388198392]; + otherwise + error('OTP:InvalidInput: %s is not a valid orbit family; try L2, P2HO1, P2HO2, P4HO1, or P4HO2', ... + orbittype) +end + +end diff --git a/toolbox/+otp/+cr3bp/+presets/L2NRHO.m b/toolbox/+otp/+cr3bp/+presets/L2NRHO.m deleted file mode 100644 index 44fa3935..00000000 --- a/toolbox/+otp/+cr3bp/+presets/L2NRHO.m +++ /dev/null @@ -1,71 +0,0 @@ -classdef L2NRHO < otp.cr3bp.CR3BPProblem - % This preset builds an $L_2$ halo orbit preset for the CR3BP based on - % a table of reference initial conditions and periods. The table - % contains $20$ entries, and thus there are $20$ different possible - % orbits given by this preset. The table of $L_2$ halo orbits is given - % as Table A.1 on page 218 of :cite:p:`Spr21`. - - methods - function obj = L2NRHO(varargin) - % Create the NRHO CR3BP problem object. - % - % Parameters - % ---------- - % Index : numeric(1, 1) - % An integer index between 1 and 20. - - T = getL2halotable(); - - p = inputParser(); - p.addParameter('Index', 1, @(x) mod(x, 1) == 0 & x > 0 & x <= size(T, 1)); - p.parse(varargin{:}); - results = p.Results; - - pic = T(results.Index, :); - - orbitalperiod = pic(1); - ic = pic(2:end); - - mE = otp.utils.PhysicalConstants.EarthMass; - mL = otp.utils.PhysicalConstants.MoonMass; - G = otp.utils.PhysicalConstants.GravitationalConstant; - - % derive mu - muE = G*mE; - muL = G*mL; - mu = muL/(muE + muL); - - y0 = [ic(1); 0; ic(2); 0; ic(3); 0]; - tspan = [0, orbitalperiod]; - params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 0); - obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); - end - end -end - - -function T = getL2halotable() -% create internal table of L_2 Halo orbits -% see -T = [ ... - 1.3632096570 1.0110350588 -0.1731500000 -0.0780141199; ... - 1.4748399512 1.0192741002 -0.1801324242 -0.0971927950; ... - 1.5872714606 1.0277926091 -0.1858044184 -0.1154896637; ... - 1.7008482705 1.0362652156 -0.1904417454 -0.1322667493; ... - 1.8155211042 1.0445681848 -0.1942338538 -0.1473971442; ... - 1.9311168544 1.0526805665 -0.1972878310 -0.1609628828; ... - 2.0474562565 1.0606277874 -0.1996480091 -0.1731020372; ... - 2.1741533495 1.0691059976 -0.2014140887 -0.1847950147; ... - 2.2915829886 1.0768767277 -0.2022559057 -0.1943508955; ... - 2.4093619266 1.0846726654 -0.2022295078 -0.2027817501; ... - 2.5273849254 1.0925906981 -0.2011567058 -0.2101017213; ... - 2.6455248145 1.1007585320 -0.1987609769 -0.2162644440; ... - 2.7635889805 1.1093498794 -0.1946155759 -0.2211327592; ... - 2.8909903824 1.1194130163 -0.1873686594 -0.2246002627; ... - 3.0073088423 1.1297344316 -0.1769810336 -0.2254855800; ... - 3.1205655022 1.1413664663 -0.1612996515 -0.2229158600; ... - 3.2266000495 1.1542349115 -0.1379744940 -0.2147411949; ... - 3.3173903769 1.1669663066 -0.1049833863 -0.1984458292; ... - 3.3833013605 1.1766385512 -0.0621463948 -0.1748356762; ... - 3.4154433338 1.1808881373 -0.0032736457 -0.1559184478]; -end From df12d4603c516e6c652fd29adc0f2c67aab74ebe Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Sat, 27 Sep 2025 13:24:41 -1000 Subject: [PATCH 10/33] I see Tupac is still alive. --- toolbox/+otp/+cr3bp/+presets/HaloOrbit.m | 2 +- toolbox/+otp/+cr3bp/CR3BPProblem.m | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m b/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m index d04f870f..1f8fe77a 100644 --- a/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m +++ b/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m @@ -60,7 +60,7 @@ % create internal table of halo orbit initial conditions % The first entry of each row is the orbital period, the second entry is % the $x$ initial condition, the third entry is the $z$ initial condition -% and the fourth entry is the $y'$ intial condition. +% and the fourth entry is the $y'$ initial condition. switch orbittype case 'L2' T = [ ... diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index b81577e5..e559c427 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -14,9 +14,9 @@ % non-dimensionalized differential equation: % % $$ - % x'' &= σ(y - x),\\ - % y'' &= ρx - y - xz,\\ - % z'' &= xy - βz,\\ + % x'' &= 2y' + \frac{\partial U}{\partial x},\\ + % y'' &= -2x' + \frac{\partial U}{\partial y},\\ + % z'' &= \frac{\partial U}{\partial z},\\ % U &= \frac{1}{2} (x^2 + y^2) + \frac{1 - \mu}{d} + \frac{mu}{r},\\ % d &= \sqrt{(x + \mu)^2 + y^2 + z^2},\\ % r &= \sqrt{(x - 1 + \mu)^2 + y^2 + z^2}, @@ -49,7 +49,7 @@ % % Example % ------- - % >>> problem = otp.cr3bp.presets.NRHO; + % >>> problem = otp.cr3bp.presets.HaloOrbit('OrbitType', 'L2', 'Index', 10); % >>> sol = model.solve(); % >>> problem.plotPhaseSpace(sol); % From 1164eb1ab29b50c69a91cc1159a0ff0ec0f848d5 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 1 Oct 2025 07:51:19 -1000 Subject: [PATCH 11/33] Change PhysicalConstants to be more in line with reality --- toolbox/+otp/+utils/PhysicalConstants.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/toolbox/+otp/+utils/PhysicalConstants.m b/toolbox/+otp/+utils/PhysicalConstants.m index 2fbeec9b..faa30fe7 100644 --- a/toolbox/+otp/+utils/PhysicalConstants.m +++ b/toolbox/+otp/+utils/PhysicalConstants.m @@ -5,11 +5,11 @@ TwoD = 2; ThreeD = 3; - GravitationalConstant = 6.67408e-11; % m^3 / (kg s^2) + GravitationalConstant = 6.6743e-11; % m^3 / (kg s^2) SunMass = 1.98855e30; % kg - EarthMass = 5.9722e24; % kg - MoonMass = 7.34767309e22; % kg + EarthMass = 5.972e24; % kg + MoonMass = 7.342e22; % kg EarthSunDistance = 1.496e11; % m EarthMoonDistance = 3.844e8; % m From d716dcf7ed8989da6153a749152d2dffe722f6e9 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 1 Oct 2025 09:04:32 -1000 Subject: [PATCH 12/33] fixed the Jacobi constant and added JCJacobian --- toolbox/+otp/+cr3bp/jacobiconstant.m | 5 ++-- toolbox/+otp/+cr3bp/jacobiconstantjacobian.m | 27 ++++++++++++++++++++ 2 files changed, 29 insertions(+), 3 deletions(-) create mode 100644 toolbox/+otp/+cr3bp/jacobiconstantjacobian.m diff --git a/toolbox/+otp/+cr3bp/jacobiconstant.m b/toolbox/+otp/+cr3bp/jacobiconstant.m index 066496b5..18cc2bbd 100644 --- a/toolbox/+otp/+cr3bp/jacobiconstant.m +++ b/toolbox/+otp/+cr3bp/jacobiconstant.m @@ -3,9 +3,8 @@ x = y(1:3, :); v = y(4:6, :); -d = sqrt((x(1, :) + mu).^2 + x(2, :)^2 + x(3, :)^2 + soft^2); -r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + x(3, :)^2 + soft^2); - +d = sqrt((x(1, :) + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); +r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); U = 0.5*(x(1, :).^2 + x(2, :).^2) + (1 - mu)./d + mu./r; diff --git a/toolbox/+otp/+cr3bp/jacobiconstantjacobian.m b/toolbox/+otp/+cr3bp/jacobiconstantjacobian.m new file mode 100644 index 00000000..eb19fe9b --- /dev/null +++ b/toolbox/+otp/+cr3bp/jacobiconstantjacobian.m @@ -0,0 +1,27 @@ +function dJdy = jacobiconstantjacobian(y, mu, soft) + +x = y(1:3, :); +v = y(4:6, :); + +d = sqrt((x(1, :) + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); +r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); + +dddx = (x(1, :) + mu)./d; +dddy = x(2, :)./d; +dddz = x(3, :)./d; +drdx = (x(1, :) - 1 + mu)./r; +drdy = x(2, :)./r; +drdz = x(3, :)./r; + +dUdx = x(1, :) - ((1 - mu).*dddx)./(d.^2) - (mu.*drdx)./(r.^2); +dUdy = x(2, :) - ((1 - mu).*dddy)./(d.^2) - (mu.*drdy)./(r.^2); +dUdz = - ((1 - mu).*dddz)./(d.^2) - (mu.*drdz)./(r.^2); + +dJdy = [reshape(2*dUdx, 1, 1, []), ... + reshape(2*dUdy, 1, 1, []), ... + reshape(2*dUdz, 1, 1, []), ... + reshape(-2*v(1, :), 1, 1, []), ... + reshape(-2*v(2, :), 1, 1, []), ... + reshape(-2*v(3, :), 1, 1, [])]; + +end From 9a023832762d40115d402b84faf2eb04bd8680eb Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 1 Oct 2025 09:04:57 -1000 Subject: [PATCH 13/33] Fixed f and added a corrected table to better work with ode45 --- .../+cr3bp/+presets/+utils/halocorrector.m | 50 ++++ toolbox/+otp/+cr3bp/+presets/HaloOrbit.m | 246 +++++++++--------- .../+cr3bp/+presets/private/haloorbits.mat | Bin 0 -> 3379 bytes toolbox/+otp/+cr3bp/CR3BPProblem.m | 5 + toolbox/+otp/+cr3bp/f.m | 24 +- 5 files changed, 197 insertions(+), 128 deletions(-) create mode 100644 toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m create mode 100644 toolbox/+otp/+cr3bp/+presets/private/haloorbits.mat diff --git a/toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m b/toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m new file mode 100644 index 00000000..7cdd52c0 --- /dev/null +++ b/toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m @@ -0,0 +1,50 @@ +orbittypes = {'L2', 'P2HO1', 'P2HO2', 'P4HO1', 'P4HO2'}; +nIndices = 20; + +orbits = struct; + +for oti = 1:numel(orbittypes) + ot = orbittypes{oti}; + otable = zeros(20, 4); + for ind = 1:nIndices + + problem = otp.cr3bp.presets.HaloOrbit('OrbitType', ot, 'Index', ind); + + y0 = problem.Y0; + y0 = [y0(1); y0(3); y0(5)]; + period = problem.TimeSpan(2); + f = problem.RHS.F; + + J = @(yc) cost(yc, y0); + nonlc = @(yc) constr(yc, f, period); + + opts = optimoptions("fmincon", "Algorithm","interior-point", ... + "FunctionTolerance", 1e-8, "ConstraintTolerance", 1e-14, 'Display','off'); + [ynew, fev] = fmincon(J, y0, [], [], [], [], [], [], nonlc, opts); + + otable(ind, :) = [period, ynew.']; + + end + orbits.(ot) = otable; +end + +function c = cost(yc, y0) + +c = sum((yc - y0).^2); + +end + +function [c, ceq] = constr(yc, f, period) + +yc = [yc(1); 0; yc(2); 0; yc(3); 0]; + +c = []; + +sol = ode45(f, [0, period], yc, ... + odeset('AbsTol', 1e-14, 'RelTol', 100*eps)); + +ceq = sum((yc - sol.y(:, end)).^2); + +end + + diff --git a/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m b/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m index 1f8fe77a..9583e661 100644 --- a/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m +++ b/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m @@ -57,124 +57,132 @@ function T = gethalotable(orbittype) -% create internal table of halo orbit initial conditions -% The first entry of each row is the orbital period, the second entry is -% the $x$ initial condition, the third entry is the $z$ initial condition -% and the fourth entry is the $y'$ initial condition. -switch orbittype - case 'L2' - T = [ ... - 1.3632096570 1.0110350588 -0.1731500000 -0.0780141199; ... - 1.4748399512 1.0192741002 -0.1801324242 -0.0971927950; ... - 1.5872714606 1.0277926091 -0.1858044184 -0.1154896637; ... - 1.7008482705 1.0362652156 -0.1904417454 -0.1322667493; ... - 1.8155211042 1.0445681848 -0.1942338538 -0.1473971442; ... - 1.9311168544 1.0526805665 -0.1972878310 -0.1609628828; ... - 2.0474562565 1.0606277874 -0.1996480091 -0.1731020372; ... - 2.1741533495 1.0691059976 -0.2014140887 -0.1847950147; ... - 2.2915829886 1.0768767277 -0.2022559057 -0.1943508955; ... - 2.4093619266 1.0846726654 -0.2022295078 -0.2027817501; ... - 2.5273849254 1.0925906981 -0.2011567058 -0.2101017213; ... - 2.6455248145 1.1007585320 -0.1987609769 -0.2162644440; ... - 2.7635889805 1.1093498794 -0.1946155759 -0.2211327592; ... - 2.8909903824 1.1194130163 -0.1873686594 -0.2246002627; ... - 3.0073088423 1.1297344316 -0.1769810336 -0.2254855800; ... - 3.1205655022 1.1413664663 -0.1612996515 -0.2229158600; ... - 3.2266000495 1.1542349115 -0.1379744940 -0.2147411949; ... - 3.3173903769 1.1669663066 -0.1049833863 -0.1984458292; ... - 3.3833013605 1.1766385512 -0.0621463948 -0.1748356762; ... - 3.4154433338 1.1808881373 -0.0032736457 -0.1559184478]; - case 'P2HO1' - T = [ ... - 2.7486723000 1.0124213729 -0.1739463747 -0.0801103670; ... - 2.8235719903 0.9686100061 -0.1684646845 -0.0555868296; ... - 2.9470915455 0.9428902766 -0.1621459176 -0.0367463465; ... - 3.0718373938 0.9280339523 -0.1570215893 -0.0278535009; ... - 3.2098634414 0.9177532010 -0.1525284780 -0.0253346758; ... - 3.3390119532 0.9116052635 -0.1492702168 -0.0275809790; ... - 3.4781621508 0.9074818076 -0.1466452811 -0.0334316688; ... - 3.6071172693 0.9053614223 -0.1449312614 -0.0411597104; ... - 3.7456075513 0.9044951611 -0.1437576475 -0.0512886367; ... - 3.8738428281 0.9047472035 -0.1431979716 -0.0619366228; ... - 4.0115944839 0.9059473646 -0.1430690015 -0.0743674595; ... - 4.1392474276 0.9077780758 -0.1433117179 -0.0865211750; ... - 4.2765261366 0.9103912789 -0.1438936424 -0.0999906823; ... - 4.4039019522 0.9133170339 -0.1446839313 -0.1126257931; ... - 4.5410586110 0.9169169198 -0.1457671148 -0.1261483024; ... - 4.6684708483 0.9206100089 -0.1469654829 -0.1384542728; ... - 4.8057987675 0.9249032213 -0.1484447194 -0.1512842987; ... - 4.9334592177 0.9291366414 -0.1499785567 -0.1627035911; ... - 5.0711130707 0.9339173119 -0.1517856322 -0.1743951579; ... - 5.2089533970 0.9388868141 -0.1537336561 -0.1854168651]; - case 'P2HO2' - T = [ ... - 4.5243776465 0.9902279661 0.0392479429 0.7431504374; ... - 4.3364693699 0.9913034068 0.0386425658 0.7476198379; ... - 4.3344588357 0.9925765901 0.0404028452 0.7275812069; ... - 4.3322427233 0.9936999422 0.0418211289 0.7119297395; ... - 4.3274956856 0.9956999819 0.0440891033 0.6876036122; ... - 4.3154119806 0.9997963248 0.0479397591 0.6474122281; ... - 4.2966706889 1.0050420842 0.0517347036 0.6074272113; ... - 4.2610464740 1.0135348381 0.0559827198 0.5575089518; ... - 4.2236953788 1.0214960279 0.0583763748 0.5199524066; ... - 4.2114971687 1.0239852380 0.0588603231 0.5093886860; ... - 4.2066140086 1.0249716121 0.0590196996 0.5053259810; ... - 4.2037569181 1.0255464522 0.0591042680 0.5029882761; ... - 4.2014068389 1.0260181267 0.0591691271 0.5010860112; ... - 4.1989300286 1.0265141820 0.0592329594 0.4991004650; ... - 4.1953346327 1.0272324868 0.0593174846 0.4962518855; ... - 4.1879862293 1.0286948426 0.0594609346 0.4905454113; ... - 4.1645925208 1.0333161410 0.0596666029 0.4732448565; ... - 4.0891753953 1.0483318888 0.0577756525 0.4227218693; ... - 3.9680003798 1.0758753256 0.0417283774 0.3430241217; ... - 3.9102015266 1.0926153718 0.0005492180 0.2990655449]; - case 'P4HO1' - T = [ ... - 8.2676291788 1.0619393673 -0.1999741894 -0.1749961447; ... - 8.3305139834 1.0681129545 -0.1976100673 -0.1868722998; ... - 8.3392497293 1.0738246899 -0.1946216607 -0.1984738200; ... - 8.3519917053 1.0786654098 -0.1919984087 -0.2083581121; ... - 8.3693678945 1.0833295002 -0.1893927253 -0.2179310067; ... - 8.3878867600 1.0871892065 -0.1871812140 -0.2258888455; ... - 8.4096026243 1.0909131908 -0.1850039744 -0.2335944979; ... - 8.4309304649 1.0940376079 -0.1831478491 -0.2400766927; ... - 8.4547978250 1.0971075676 -0.1813014490 -0.2464568317; ... - 8.4775801091 1.0997316965 -0.1797082140 -0.2519152440; ... - 8.5026193458 1.1023550439 -0.1781043165 -0.2573726538; ... - 8.5262347341 1.1046324969 -0.1767049403 -0.2621078568; ... - 8.5519768095 1.1069405119 -0.1752821493 -0.2669009993; ... - 8.5761129683 1.1089683638 -0.1740297442 -0.2711051374; ... - 8.6023105873 1.1110450669 -0.1727463692 -0.2754011052; ... - 8.6267952587 1.1128866230 -0.1716087730 -0.2792006165; ... - 8.6533063618 1.1147879601 -0.1704357869 -0.2831115888; ... - 8.6780368548 1.1164862743 -0.1693902247 -0.2865931191; ... - 8.7047739081 1.1182510772 -0.1683067141 -0.2901975105; ... - 8.7316043012 1.1199566179 -0.1672632043 -0.2936663490]; - case 'P4HO2' - T = [ ... - 9.0487552928 1.0749349573 -0.2021215571 -0.1920690379; ... - 9.0396341760 1.0578442875 -0.2047569352 -0.1678532831; ... - 9.0542791209 1.0524484578 -0.2047233749 -0.1606194256; ... - 9.0601028655 1.0505487899 -0.2045610040 -0.1581773718; ... - 9.0646190673 1.0491368344 -0.2043821958 -0.1564064537; ... - 9.0707205473 1.0472961590 -0.2040685317 -0.1541623227; ... - 9.0786237393 1.0450028327 -0.2035390634 -0.1514838937; ... - 9.0888613168 1.0421535217 -0.2026443039 -0.1483728351; ... - 9.1037491359 1.0382266494 -0.2009270031 -0.1445796433; ... - 9.1306372229 1.0319435819 -0.1968531223 -0.1401578788; ... - 9.2021506431 1.0215863525 -0.1864218292 -0.1399647167; ... - 9.3830140082 1.0126738542 -0.1723454500 -0.1560887631; ... - 9.5440425474 1.0113733655 -0.1648460585 -0.1739113299; ... - 9.6990195652 1.0129644895 -0.1593270671 -0.1918297866; ... - 10.2961547121 1.0286792032 -0.1457428883 -0.2522546927; ... - 11.7792576999 1.0624281188 -0.1343562748 -0.3275531781; ... - 13.7831848966 1.0850769369 -0.1379181777 -0.3653699405; ... - 16.0024825055 1.0955445136 -0.1591111395 -0.3767847569; ... - 18.1058513583 1.0954718438 -0.1988427540 -0.3675757435; ... - 20.0406128078 1.0829659022 -0.2591580427 -0.3388198392]; - otherwise - error('OTP:InvalidInput: %s is not a valid orbit family; try L2, P2HO1, P2HO2, P4HO1, or P4HO2', ... - orbittype) -end +st = load('haloorbits.mat'); + +T = st.orbits.(orbittype); end + +% below is a backup function to restart the Halo orbit table. +% function T = gethalotable(orbittype) +% % create internal table of halo orbit initial conditions +% % The first entry of each row is the orbital period, the second entry is +% % the $x$ initial condition, the third entry is the $z$ initial condition +% % and the fourth entry is the $y'$ initial condition. +% switch orbittype +% case 'L2' +% T = [ ... +% 1.3632096570 1.0110350588 -0.1731500000 -0.0780141199; ... +% 1.4748399512 1.0192741002 -0.1801324242 -0.0971927950; ... +% 1.5872714606 1.0277926091 -0.1858044184 -0.1154896637; ... +% 1.7008482705 1.0362652156 -0.1904417454 -0.1322667493; ... +% 1.8155211042 1.0445681848 -0.1942338538 -0.1473971442; ... +% 1.9311168544 1.0526805665 -0.1972878310 -0.1609628828; ... +% 2.0474562565 1.0606277874 -0.1996480091 -0.1731020372; ... +% 2.1741533495 1.0691059976 -0.2014140887 -0.1847950147; ... +% 2.2915829886 1.0768767277 -0.2022559057 -0.1943508955; ... +% 2.4093619266 1.0846726654 -0.2022295078 -0.2027817501; ... +% 2.5273849254 1.0925906981 -0.2011567058 -0.2101017213; ... +% 2.6455248145 1.1007585320 -0.1987609769 -0.2162644440; ... +% 2.7635889805 1.1093498794 -0.1946155759 -0.2211327592; ... +% 2.8909903824 1.1194130163 -0.1873686594 -0.2246002627; ... +% 3.0073088423 1.1297344316 -0.1769810336 -0.2254855800; ... +% 3.1205655022 1.1413664663 -0.1612996515 -0.2229158600; ... +% 3.2266000495 1.1542349115 -0.1379744940 -0.2147411949; ... +% 3.3173903769 1.1669663066 -0.1049833863 -0.1984458292; ... +% 3.3833013605 1.1766385512 -0.0621463948 -0.1748356762; ... +% 3.4154433338 1.1808881373 -0.0032736457 -0.1559184478]; +% case 'P2HO1' +% T = [ ... +% 2.7486723000 1.0124213729 -0.1739463747 -0.0801103670; ... +% 2.8235719903 0.9686100061 -0.1684646845 -0.0555868296; ... +% 2.9470915455 0.9428902766 -0.1621459176 -0.0367463465; ... +% 3.0718373938 0.9280339523 -0.1570215893 -0.0278535009; ... +% 3.2098634414 0.9177532010 -0.1525284780 -0.0253346758; ... +% 3.3390119532 0.9116052635 -0.1492702168 -0.0275809790; ... +% 3.4781621508 0.9074818076 -0.1466452811 -0.0334316688; ... +% 3.6071172693 0.9053614223 -0.1449312614 -0.0411597104; ... +% 3.7456075513 0.9044951611 -0.1437576475 -0.0512886367; ... +% 3.8738428281 0.9047472035 -0.1431979716 -0.0619366228; ... +% 4.0115944839 0.9059473646 -0.1430690015 -0.0743674595; ... +% 4.1392474276 0.9077780758 -0.1433117179 -0.0865211750; ... +% 4.2765261366 0.9103912789 -0.1438936424 -0.0999906823; ... +% 4.4039019522 0.9133170339 -0.1446839313 -0.1126257931; ... +% 4.5410586110 0.9169169198 -0.1457671148 -0.1261483024; ... +% 4.6684708483 0.9206100089 -0.1469654829 -0.1384542728; ... +% 4.8057987675 0.9249032213 -0.1484447194 -0.1512842987; ... +% 4.9334592177 0.9291366414 -0.1499785567 -0.1627035911; ... +% 5.0711130707 0.9339173119 -0.1517856322 -0.1743951579; ... +% 5.2089533970 0.9388868141 -0.1537336561 -0.1854168651]; +% case 'P2HO2' +% T = [ ... +% 4.5243776465 0.9902279661 0.0392479429 0.7431504374; ... +% 4.3364693699 0.9913034068 0.0386425658 0.7476198379; ... +% 4.3344588357 0.9925765901 0.0404028452 0.7275812069; ... +% 4.3322427233 0.9936999422 0.0418211289 0.7119297395; ... +% 4.3274956856 0.9956999819 0.0440891033 0.6876036122; ... +% 4.3154119806 0.9997963248 0.0479397591 0.6474122281; ... +% 4.2966706889 1.0050420842 0.0517347036 0.6074272113; ... +% 4.2610464740 1.0135348381 0.0559827198 0.5575089518; ... +% 4.2236953788 1.0214960279 0.0583763748 0.5199524066; ... +% 4.2114971687 1.0239852380 0.0588603231 0.5093886860; ... +% 4.2066140086 1.0249716121 0.0590196996 0.5053259810; ... +% 4.2037569181 1.0255464522 0.0591042680 0.5029882761; ... +% 4.2014068389 1.0260181267 0.0591691271 0.5010860112; ... +% 4.1989300286 1.0265141820 0.0592329594 0.4991004650; ... +% 4.1953346327 1.0272324868 0.0593174846 0.4962518855; ... +% 4.1879862293 1.0286948426 0.0594609346 0.4905454113; ... +% 4.1645925208 1.0333161410 0.0596666029 0.4732448565; ... +% 4.0891753953 1.0483318888 0.0577756525 0.4227218693; ... +% 3.9680003798 1.0758753256 0.0417283774 0.3430241217; ... +% 3.9102015266 1.0926153718 0.0005492180 0.2990655449]; +% case 'P4HO1' +% T = [ ... +% 8.2676291788 1.0619393673 -0.1999741894 -0.1749961447; ... +% 8.3305139834 1.0681129545 -0.1976100673 -0.1868722998; ... +% 8.3392497293 1.0738246899 -0.1946216607 -0.1984738200; ... +% 8.3519917053 1.0786654098 -0.1919984087 -0.2083581121; ... +% 8.3693678945 1.0833295002 -0.1893927253 -0.2179310067; ... +% 8.3878867600 1.0871892065 -0.1871812140 -0.2258888455; ... +% 8.4096026243 1.0909131908 -0.1850039744 -0.2335944979; ... +% 8.4309304649 1.0940376079 -0.1831478491 -0.2400766927; ... +% 8.4547978250 1.0971075676 -0.1813014490 -0.2464568317; ... +% 8.4775801091 1.0997316965 -0.1797082140 -0.2519152440; ... +% 8.5026193458 1.1023550439 -0.1781043165 -0.2573726538; ... +% 8.5262347341 1.1046324969 -0.1767049403 -0.2621078568; ... +% 8.5519768095 1.1069405119 -0.1752821493 -0.2669009993; ... +% 8.5761129683 1.1089683638 -0.1740297442 -0.2711051374; ... +% 8.6023105873 1.1110450669 -0.1727463692 -0.2754011052; ... +% 8.6267952587 1.1128866230 -0.1716087730 -0.2792006165; ... +% 8.6533063618 1.1147879601 -0.1704357869 -0.2831115888; ... +% 8.6780368548 1.1164862743 -0.1693902247 -0.2865931191; ... +% 8.7047739081 1.1182510772 -0.1683067141 -0.2901975105; ... +% 8.7316043012 1.1199566179 -0.1672632043 -0.2936663490]; +% case 'P4HO2' +% T = [ ... +% 9.0487552928 1.0749349573 -0.2021215571 -0.1920690379; ... +% 9.0396341760 1.0578442875 -0.2047569352 -0.1678532831; ... +% 9.0542791209 1.0524484578 -0.2047233749 -0.1606194256; ... +% 9.0601028655 1.0505487899 -0.2045610040 -0.1581773718; ... +% 9.0646190673 1.0491368344 -0.2043821958 -0.1564064537; ... +% 9.0707205473 1.0472961590 -0.2040685317 -0.1541623227; ... +% 9.0786237393 1.0450028327 -0.2035390634 -0.1514838937; ... +% 9.0888613168 1.0421535217 -0.2026443039 -0.1483728351; ... +% 9.1037491359 1.0382266494 -0.2009270031 -0.1445796433; ... +% 9.1306372229 1.0319435819 -0.1968531223 -0.1401578788; ... +% 9.2021506431 1.0215863525 -0.1864218292 -0.1399647167; ... +% 9.3830140082 1.0126738542 -0.1723454500 -0.1560887631; ... +% 9.5440425474 1.0113733655 -0.1648460585 -0.1739113299; ... +% 9.6990195652 1.0129644895 -0.1593270671 -0.1918297866; ... +% 10.2961547121 1.0286792032 -0.1457428883 -0.2522546927; ... +% 11.7792576999 1.0624281188 -0.1343562748 -0.3275531781; ... +% 13.7831848966 1.0850769369 -0.1379181777 -0.3653699405; ... +% 16.0024825055 1.0955445136 -0.1591111395 -0.3767847569; ... +% 18.1058513583 1.0954718438 -0.1988427540 -0.3675757435; ... +% 20.0406128078 1.0829659022 -0.2591580427 -0.3388198392]; +% otherwise +% error('OTP:InvalidInput: %s is not a valid orbit family; try L2, P2HO1, P2HO2, P4HO1, or P4HO2', ... +% orbittype) +% end +% +% end diff --git a/toolbox/+otp/+cr3bp/+presets/private/haloorbits.mat b/toolbox/+otp/+cr3bp/+presets/private/haloorbits.mat new file mode 100644 index 0000000000000000000000000000000000000000..9a1e5bd68f874040f51507fe976b067b1d90da97 GIT binary patch literal 3379 zcmV-34b1XQK~zjZLLfCRFd$7qR4ry{Y-KDUP;6mzW^ZzBIv`C!LqRq)EFeR2Wnpw> zWFT*DIv`hNWFSvtbRZxxATT*PFfckaFd#B8GBqF|ARr(hARr(hARr(hARr(hARr(h zARr(hARr(hARr(hARr(hARr(B00000000000ZB~{000213;+OloSjyAIMrMCHx-dO zq7pJirY0_pYHdoX3~`B4QZkj&RY|E(C{uOmmZB>ql?)XxU1D(@hjWfXrd$~_q#V9q z*KkUQSNHcm&+mQSKks^;&)#e8wf0{Be4fwRb3{Z$1hOI`GpBc{X^4sZpUwDBEc##g zQRC--$mi6_A9;zt;~(d$)37r7ach^6rLE!LtC7eq-Vd8Kjgaa3h>%`>ZuKbCR@brVgmF+iYR?St_zH##6*4ncPrx`QsO7+HAsCXPK-PbSp(94X>kRUz|bOo~vGZ?olVW|6sXw~mWPNt3RQVrvmO z86vjGZKTd#hNRi$tuR?mLEJMj{L6-dcgvj&ntq|6!1BPZsGAfF+caH^&Z9ut?$IQx ziGub7e|b9r1v6GIu9&4vgF~%H)ra*o=$D?;=I*6|Yp8f-O%M%QM?<;kNi?_`?p^vP z(f~)VhB5kSSaUd0#!-G4UaZ%9KV&`(x%*{!sb_}4{MRS)yl2A@ds_a4Pu(!ItF(33 zb`QfX0Z8@eve9-av3@LDfzBb)`DZICaMKB1v&Vr-)cHPm{eDX&rh@#u4B+60E8L~R z7!Gbb6tbg<&B44ELHa&j9E`nCY7*JWLCthiBe!A>YApWttj&*u8Pc}VQlBdEkZd^5XZT^buX@vN5(sVs~{kSTpy#Qms+PfRoW}~Id%P6JVBN(yj zsZ{xR4km>e?vC$xg>&k7@$uJ-uypW6SP;Ji*GCQt*0wV7)PvAjId@oCIq*Y)+TNF33y%LW%)=Ad?#X-Hmj>xo+Tr5At6uX|u!?fCy@8?)oVTCwj_2c3y zR6S$ubFZZeKh-F1`29*1`qWe_Wt`{Xjr1gwB2Nw`_=s`j4zkgV74SrVv<&|#Z;!aD z^=%W*Ae=7;`!jCMB*g>26sE>Y60T@vkXpPHS#)`4S5@XLaz9?JcW0h7X&%3rxwTG) z{Hy7*m0Qbf648Enq3F0Q`MgEKHbppx_&-VYZI@vXcekvLTh|yQr?@*^Pg#!Gj4nT$ za8Hgnv=nWT(Vs_RR3Q0s#yqlZTd(bsb@Id|`+4DqLU|(D>)u*xsX$hUD|L=HD3H+) zX|KDDC^)mb-*{kh5U$*?U9xG}AnYqlJa_#i9}YO?@tj@x;2RaEGHJkv{!B*jdTBna z8OuA_As7J07UNHq^Z=~bXlmR1*8t=tN84rz2cXby_sY!Ke7Jcc>A8X?A5!SZSef;F zuNSh{#b4{Zm>S>CN?Vv-PmpB2etFi*Lp(d>+L`A(vVnBZ{2IWe(l#`fmR9X zY*~=6_^K4OESzm`?kmFpL7wxacPzZRZ5V`GBkLydO@75|vGq)YGAnBqJ#$7bcR zXlFUHgBj*_ezzPkPbu_U;~+;~4EjpC*vb)8_P<;v%;m@#$-VGtkL!o0Zm4DySES8`k(gNW_4dHkq`K~flc2^-%sh^(sVu3R32NOW7hnNi6g z0j!+?!6gi0CZJp>;UJk#;j z5UfbEoYXHH0)rV|8(gY};8|ava>0io_-V8((Ry+S^c-~hQWPkd$hdpRdkqCUwzeCK zI#X~XgcV?Phyu<`%L}HS6i6xw+!;yHv6$BF(wPkk^A@DSq?x2LkNpr7xxJ%- znW-?l?MUEWuSa0|b;FrkH_|}F`J{nX{SV%D_RE6Z zjU(bAQCZ->z1?aMj9Qm>s{jyGf)M+I*y3Z{>Hx_L0;m$rB$zH~Qxm#imYJrn-9a*uxG8dmS+R z(W(Om6F+q(8gzh5tEl5?#SZWhOPXw#>HupsC#$t|JBTEI{`DHO9U8|1W#l~CVUf4v zy@c&;z*}X*sgrGiW4h;(=Z-YMfT#Pd_e<*FpYrzn8-*7$Bh^WP_JsPnr|QI4`s0fwg1~8Q+p>-4Wd4J=izZ#4bszRmz$)dL3mb2w+WVNkgl((vWJW{h*Wi#%@?9U zEIhjnF4}02$d{3ibX+uuLBTTh_ud*L<0`xHW`G7+eW8JuBcvd!Ky~5w1vIqRzmYw- zf`%yeu?b^y8Xh@*65r)W!)5KzC^Ziniu)G6aXUqW*O0XB@=zKARD~}~Vrhu?7UYH` z&~QM?y$~MKpf^{0{>=;;QeN4xrSfQ)*pwUXR7}Hs1MT0N*)$9-`=;mdmWCT0Gd3GH z(l9Ihs@7~}*ulZ-YUTwpJ=X=a)S2yZ zmB?{`jXkE7xI3vv>gY%XE-tNJ%gC(2QX$u6!-)zcOC+Z9f2u&&j&z66DK_47FbJl& zZ0t5IZtn&jnx+@5y~=Ea{ZjezW98`Fe|I3ix*W?J#6P=d zmZO60yy%aS<+xnR zcsCCZZ{B2Q(_e*Cfx$cLd*0x0>eVjh1#dBBqqL=SbTvBj>NzTEHCW6DREZ3%!Oug7 zglg;>ED`;s;i*C`X4^z}#q6oY#aDTvbaE}$O0(BlhbJD^bSWf z9h(oty~7NsZZ_t<`^S3gudh3;x0)o?fAioD1x+$z)4a1>Wld7!clk}ZswSxnSyVD# zLzBq;&I>Zp(j>a${jO%Znq=8L@eRStHA#fa&xxugnxsKk-PvWXNnVVK?g;VHB;U<0 z`Y5q9$=>5ni(?inA!!L`x_dm95EI#@w^MFw5yq*I)^FEz$OTiEB;gGMVli~(#RHL5 z@V!27Vt;#wEpgT;I_xk zm;zZlkLLcB6m-U$_&AwRVBF4^TG)w@=u5 zPxD(PC|1XmhVPosdq-P2NcOvnBpP$^#^>13^q;t>b77;vbp;oN<$s!7*5%@64;Mo- zRW2@xn2hDm=HjK5Z#^CiaqzL0y=Oo^2Pc+tU#Y8<1_rAN!OHso5OR(J9QVfZ?QdPL0g%<=S*2@z)7~SUg zMbV%J$I|}tw*0LUAD&!gi9zLP6`|sse~pF0?0s@CH?c6<-hWQI1Pi~VICp6jmf_C$ zxL*&rmtpJ2h_DkvCTi)+2bzR3krs=kw@{_n!S%ZLt-ln%8uGm!3d(SJcV}wKIyRm? zy*F7TrV?k$*6 Date: Wed, 1 Oct 2025 10:04:25 -1000 Subject: [PATCH 14/33] Corrected orbit types --- .../+cr3bp/+presets/+utils/halocorrector.m | 6 +- toolbox/+otp/+cr3bp/+presets/HaloOrbit.m | 262 +++++++++--------- toolbox/+otp/+cr3bp/f.m | 10 +- 3 files changed, 145 insertions(+), 133 deletions(-) diff --git a/toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m b/toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m index 7cdd52c0..eed0a26f 100644 --- a/toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m +++ b/toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m @@ -4,11 +4,11 @@ orbits = struct; for oti = 1:numel(orbittypes) - ot = orbittypes{oti}; + orbittype = orbittypes{oti}; otable = zeros(20, 4); for ind = 1:nIndices - problem = otp.cr3bp.presets.HaloOrbit('OrbitType', ot, 'Index', ind); + problem = otp.cr3bp.presets.HaloOrbit('OrbitType', orbittype, 'Index', ind); y0 = problem.Y0; y0 = [y0(1); y0(3); y0(5)]; @@ -25,7 +25,7 @@ otable(ind, :) = [period, ynew.']; end - orbits.(ot) = otable; + orbits.(orbittype) = otable; end function c = cost(yc, y0) diff --git a/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m b/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m index 9583e661..db56f585 100644 --- a/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m +++ b/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m @@ -22,16 +22,25 @@ % e. 'P4HO2' for the period quadrupling halo orbit of type 2 % Index : numeric(1, 1) % An integer index between 1 and 20. + % Corrected : boolean + % True (default) uses a modified version to better work + % with the physical constants and matlab time integrators. + % False uses the original table from :cite:p:`Spr21`. p = inputParser(); p.addParameter('OrbitType', 'L2', ... @(x) strcmp(x, 'L2') || strcmp(x, 'P2HO1') || strcmp(x, 'P2HO2') ... || strcmp(x, 'P4HO1') || strcmp(x, 'P4HO2')); p.addParameter('Index', 1, @(x) mod(x, 1) == 0 && x > 0 && x <= 20); + p.addParameter('Corrected', true, @islogical); p.parse(varargin{:}); results = p.Results; - T = gethalotable(results.OrbitType); + if results.Corrected + T = gethalotablecorrected(results.OrbitType); + else + T = gethalotable(results.OrbitType); + end pic = T(results.Index, :); @@ -56,133 +65,136 @@ end -function T = gethalotable(orbittype) +function T = gethalotablecorrected(orbittype) +% create internal table of halo orbit initial conditions +% The first entry of each row is the orbital period, the second entry is +% the $x$ initial condition, the third entry is the $z$ initial condition +% and the fourth entry is the $y'$ initial condition. st = load('haloorbits.mat'); T = st.orbits.(orbittype); end -% below is a backup function to restart the Halo orbit table. -% function T = gethalotable(orbittype) -% % create internal table of halo orbit initial conditions -% % The first entry of each row is the orbital period, the second entry is -% % the $x$ initial condition, the third entry is the $z$ initial condition -% % and the fourth entry is the $y'$ initial condition. -% switch orbittype -% case 'L2' -% T = [ ... -% 1.3632096570 1.0110350588 -0.1731500000 -0.0780141199; ... -% 1.4748399512 1.0192741002 -0.1801324242 -0.0971927950; ... -% 1.5872714606 1.0277926091 -0.1858044184 -0.1154896637; ... -% 1.7008482705 1.0362652156 -0.1904417454 -0.1322667493; ... -% 1.8155211042 1.0445681848 -0.1942338538 -0.1473971442; ... -% 1.9311168544 1.0526805665 -0.1972878310 -0.1609628828; ... -% 2.0474562565 1.0606277874 -0.1996480091 -0.1731020372; ... -% 2.1741533495 1.0691059976 -0.2014140887 -0.1847950147; ... -% 2.2915829886 1.0768767277 -0.2022559057 -0.1943508955; ... -% 2.4093619266 1.0846726654 -0.2022295078 -0.2027817501; ... -% 2.5273849254 1.0925906981 -0.2011567058 -0.2101017213; ... -% 2.6455248145 1.1007585320 -0.1987609769 -0.2162644440; ... -% 2.7635889805 1.1093498794 -0.1946155759 -0.2211327592; ... -% 2.8909903824 1.1194130163 -0.1873686594 -0.2246002627; ... -% 3.0073088423 1.1297344316 -0.1769810336 -0.2254855800; ... -% 3.1205655022 1.1413664663 -0.1612996515 -0.2229158600; ... -% 3.2266000495 1.1542349115 -0.1379744940 -0.2147411949; ... -% 3.3173903769 1.1669663066 -0.1049833863 -0.1984458292; ... -% 3.3833013605 1.1766385512 -0.0621463948 -0.1748356762; ... -% 3.4154433338 1.1808881373 -0.0032736457 -0.1559184478]; -% case 'P2HO1' -% T = [ ... -% 2.7486723000 1.0124213729 -0.1739463747 -0.0801103670; ... -% 2.8235719903 0.9686100061 -0.1684646845 -0.0555868296; ... -% 2.9470915455 0.9428902766 -0.1621459176 -0.0367463465; ... -% 3.0718373938 0.9280339523 -0.1570215893 -0.0278535009; ... -% 3.2098634414 0.9177532010 -0.1525284780 -0.0253346758; ... -% 3.3390119532 0.9116052635 -0.1492702168 -0.0275809790; ... -% 3.4781621508 0.9074818076 -0.1466452811 -0.0334316688; ... -% 3.6071172693 0.9053614223 -0.1449312614 -0.0411597104; ... -% 3.7456075513 0.9044951611 -0.1437576475 -0.0512886367; ... -% 3.8738428281 0.9047472035 -0.1431979716 -0.0619366228; ... -% 4.0115944839 0.9059473646 -0.1430690015 -0.0743674595; ... -% 4.1392474276 0.9077780758 -0.1433117179 -0.0865211750; ... -% 4.2765261366 0.9103912789 -0.1438936424 -0.0999906823; ... -% 4.4039019522 0.9133170339 -0.1446839313 -0.1126257931; ... -% 4.5410586110 0.9169169198 -0.1457671148 -0.1261483024; ... -% 4.6684708483 0.9206100089 -0.1469654829 -0.1384542728; ... -% 4.8057987675 0.9249032213 -0.1484447194 -0.1512842987; ... -% 4.9334592177 0.9291366414 -0.1499785567 -0.1627035911; ... -% 5.0711130707 0.9339173119 -0.1517856322 -0.1743951579; ... -% 5.2089533970 0.9388868141 -0.1537336561 -0.1854168651]; -% case 'P2HO2' -% T = [ ... -% 4.5243776465 0.9902279661 0.0392479429 0.7431504374; ... -% 4.3364693699 0.9913034068 0.0386425658 0.7476198379; ... -% 4.3344588357 0.9925765901 0.0404028452 0.7275812069; ... -% 4.3322427233 0.9936999422 0.0418211289 0.7119297395; ... -% 4.3274956856 0.9956999819 0.0440891033 0.6876036122; ... -% 4.3154119806 0.9997963248 0.0479397591 0.6474122281; ... -% 4.2966706889 1.0050420842 0.0517347036 0.6074272113; ... -% 4.2610464740 1.0135348381 0.0559827198 0.5575089518; ... -% 4.2236953788 1.0214960279 0.0583763748 0.5199524066; ... -% 4.2114971687 1.0239852380 0.0588603231 0.5093886860; ... -% 4.2066140086 1.0249716121 0.0590196996 0.5053259810; ... -% 4.2037569181 1.0255464522 0.0591042680 0.5029882761; ... -% 4.2014068389 1.0260181267 0.0591691271 0.5010860112; ... -% 4.1989300286 1.0265141820 0.0592329594 0.4991004650; ... -% 4.1953346327 1.0272324868 0.0593174846 0.4962518855; ... -% 4.1879862293 1.0286948426 0.0594609346 0.4905454113; ... -% 4.1645925208 1.0333161410 0.0596666029 0.4732448565; ... -% 4.0891753953 1.0483318888 0.0577756525 0.4227218693; ... -% 3.9680003798 1.0758753256 0.0417283774 0.3430241217; ... -% 3.9102015266 1.0926153718 0.0005492180 0.2990655449]; -% case 'P4HO1' -% T = [ ... -% 8.2676291788 1.0619393673 -0.1999741894 -0.1749961447; ... -% 8.3305139834 1.0681129545 -0.1976100673 -0.1868722998; ... -% 8.3392497293 1.0738246899 -0.1946216607 -0.1984738200; ... -% 8.3519917053 1.0786654098 -0.1919984087 -0.2083581121; ... -% 8.3693678945 1.0833295002 -0.1893927253 -0.2179310067; ... -% 8.3878867600 1.0871892065 -0.1871812140 -0.2258888455; ... -% 8.4096026243 1.0909131908 -0.1850039744 -0.2335944979; ... -% 8.4309304649 1.0940376079 -0.1831478491 -0.2400766927; ... -% 8.4547978250 1.0971075676 -0.1813014490 -0.2464568317; ... -% 8.4775801091 1.0997316965 -0.1797082140 -0.2519152440; ... -% 8.5026193458 1.1023550439 -0.1781043165 -0.2573726538; ... -% 8.5262347341 1.1046324969 -0.1767049403 -0.2621078568; ... -% 8.5519768095 1.1069405119 -0.1752821493 -0.2669009993; ... -% 8.5761129683 1.1089683638 -0.1740297442 -0.2711051374; ... -% 8.6023105873 1.1110450669 -0.1727463692 -0.2754011052; ... -% 8.6267952587 1.1128866230 -0.1716087730 -0.2792006165; ... -% 8.6533063618 1.1147879601 -0.1704357869 -0.2831115888; ... -% 8.6780368548 1.1164862743 -0.1693902247 -0.2865931191; ... -% 8.7047739081 1.1182510772 -0.1683067141 -0.2901975105; ... -% 8.7316043012 1.1199566179 -0.1672632043 -0.2936663490]; -% case 'P4HO2' -% T = [ ... -% 9.0487552928 1.0749349573 -0.2021215571 -0.1920690379; ... -% 9.0396341760 1.0578442875 -0.2047569352 -0.1678532831; ... -% 9.0542791209 1.0524484578 -0.2047233749 -0.1606194256; ... -% 9.0601028655 1.0505487899 -0.2045610040 -0.1581773718; ... -% 9.0646190673 1.0491368344 -0.2043821958 -0.1564064537; ... -% 9.0707205473 1.0472961590 -0.2040685317 -0.1541623227; ... -% 9.0786237393 1.0450028327 -0.2035390634 -0.1514838937; ... -% 9.0888613168 1.0421535217 -0.2026443039 -0.1483728351; ... -% 9.1037491359 1.0382266494 -0.2009270031 -0.1445796433; ... -% 9.1306372229 1.0319435819 -0.1968531223 -0.1401578788; ... -% 9.2021506431 1.0215863525 -0.1864218292 -0.1399647167; ... -% 9.3830140082 1.0126738542 -0.1723454500 -0.1560887631; ... -% 9.5440425474 1.0113733655 -0.1648460585 -0.1739113299; ... -% 9.6990195652 1.0129644895 -0.1593270671 -0.1918297866; ... -% 10.2961547121 1.0286792032 -0.1457428883 -0.2522546927; ... -% 11.7792576999 1.0624281188 -0.1343562748 -0.3275531781; ... -% 13.7831848966 1.0850769369 -0.1379181777 -0.3653699405; ... -% 16.0024825055 1.0955445136 -0.1591111395 -0.3767847569; ... -% 18.1058513583 1.0954718438 -0.1988427540 -0.3675757435; ... -% 20.0406128078 1.0829659022 -0.2591580427 -0.3388198392]; -% otherwise -% error('OTP:InvalidInput: %s is not a valid orbit family; try L2, P2HO1, P2HO2, P4HO1, or P4HO2', ... -% orbittype) -% end -% -% end +function T = gethalotable(orbittype) +% create internal table of halo orbit initial conditions +% The first entry of each row is the orbital period, the second entry is +% the $x$ initial condition, the third entry is the $z$ initial condition +% and the fourth entry is the $y'$ initial condition. +switch orbittype + case 'L2' + T = [ ... + 1.3632096570 1.0110350588 -0.1731500000 -0.0780141199; ... + 1.4748399512 1.0192741002 -0.1801324242 -0.0971927950; ... + 1.5872714606 1.0277926091 -0.1858044184 -0.1154896637; ... + 1.7008482705 1.0362652156 -0.1904417454 -0.1322667493; ... + 1.8155211042 1.0445681848 -0.1942338538 -0.1473971442; ... + 1.9311168544 1.0526805665 -0.1972878310 -0.1609628828; ... + 2.0474562565 1.0606277874 -0.1996480091 -0.1731020372; ... + 2.1741533495 1.0691059976 -0.2014140887 -0.1847950147; ... + 2.2915829886 1.0768767277 -0.2022559057 -0.1943508955; ... + 2.4093619266 1.0846726654 -0.2022295078 -0.2027817501; ... + 2.5273849254 1.0925906981 -0.2011567058 -0.2101017213; ... + 2.6455248145 1.1007585320 -0.1987609769 -0.2162644440; ... + 2.7635889805 1.1093498794 -0.1946155759 -0.2211327592; ... + 2.8909903824 1.1194130163 -0.1873686594 -0.2246002627; ... + 3.0073088423 1.1297344316 -0.1769810336 -0.2254855800; ... + 3.1205655022 1.1413664663 -0.1612996515 -0.2229158600; ... + 3.2266000495 1.1542349115 -0.1379744940 -0.2147411949; ... + 3.3173903769 1.1669663066 -0.1049833863 -0.1984458292; ... + 3.3833013605 1.1766385512 -0.0621463948 -0.1748356762; ... + 3.4154433338 1.1808881373 -0.0032736457 -0.1559184478]; + case 'P2HO1' + T = [ ... + 2.7486723000 1.0124213729 -0.1739463747 -0.0801103670; ... + 2.8235719903 0.9686100061 -0.1684646845 -0.0555868296; ... + 2.9470915455 0.9428902766 -0.1621459176 -0.0367463465; ... + 3.0718373938 0.9280339523 -0.1570215893 -0.0278535009; ... + 3.2098634414 0.9177532010 -0.1525284780 -0.0253346758; ... + 3.3390119532 0.9116052635 -0.1492702168 -0.0275809790; ... + 3.4781621508 0.9074818076 -0.1466452811 -0.0334316688; ... + 3.6071172693 0.9053614223 -0.1449312614 -0.0411597104; ... + 3.7456075513 0.9044951611 -0.1437576475 -0.0512886367; ... + 3.8738428281 0.9047472035 -0.1431979716 -0.0619366228; ... + 4.0115944839 0.9059473646 -0.1430690015 -0.0743674595; ... + 4.1392474276 0.9077780758 -0.1433117179 -0.0865211750; ... + 4.2765261366 0.9103912789 -0.1438936424 -0.0999906823; ... + 4.4039019522 0.9133170339 -0.1446839313 -0.1126257931; ... + 4.5410586110 0.9169169198 -0.1457671148 -0.1261483024; ... + 4.6684708483 0.9206100089 -0.1469654829 -0.1384542728; ... + 4.8057987675 0.9249032213 -0.1484447194 -0.1512842987; ... + 4.9334592177 0.9291366414 -0.1499785567 -0.1627035911; ... + 5.0711130707 0.9339173119 -0.1517856322 -0.1743951579; ... + 5.2089533970 0.9388868141 -0.1537336561 -0.1854168651]; + case 'P2HO2' + T = [ ... + 4.5243776465 0.9902279661 0.0392479429 0.7431504374; ... + 4.3364693699 0.9913034068 0.0386425658 0.7476198379; ... + 4.3344588357 0.9925765901 0.0404028452 0.7275812069; ... + 4.3322427233 0.9936999422 0.0418211289 0.7119297395; ... + 4.3274956856 0.9956999819 0.0440891033 0.6876036122; ... + 4.3154119806 0.9997963248 0.0479397591 0.6474122281; ... + 4.2966706889 1.0050420842 0.0517347036 0.6074272113; ... + 4.2610464740 1.0135348381 0.0559827198 0.5575089518; ... + 4.2236953788 1.0214960279 0.0583763748 0.5199524066; ... + 4.2114971687 1.0239852380 0.0588603231 0.5093886860; ... + 4.2066140086 1.0249716121 0.0590196996 0.5053259810; ... + 4.2037569181 1.0255464522 0.0591042680 0.5029882761; ... + 4.2014068389 1.0260181267 0.0591691271 0.5010860112; ... + 4.1989300286 1.0265141820 0.0592329594 0.4991004650; ... + 4.1953346327 1.0272324868 0.0593174846 0.4962518855; ... + 4.1879862293 1.0286948426 0.0594609346 0.4905454113; ... + 4.1645925208 1.0333161410 0.0596666029 0.4732448565; ... + 4.0891753953 1.0483318888 0.0577756525 0.4227218693; ... + 3.9680003798 1.0758753256 0.0417283774 0.3430241217; ... + 3.9102015266 1.0926153718 0.0005492180 0.2990655449]; + case 'P4HO1' + T = [ ... + 8.2676291788 1.0619393673 -0.1999741894 -0.1749961447; ... + 8.3305139834 1.0681129545 -0.1976100673 -0.1868722998; ... + 8.3392497293 1.0738246899 -0.1946216607 -0.1984738200; ... + 8.3519917053 1.0786654098 -0.1919984087 -0.2083581121; ... + 8.3693678945 1.0833295002 -0.1893927253 -0.2179310067; ... + 8.3878867600 1.0871892065 -0.1871812140 -0.2258888455; ... + 8.4096026243 1.0909131908 -0.1850039744 -0.2335944979; ... + 8.4309304649 1.0940376079 -0.1831478491 -0.2400766927; ... + 8.4547978250 1.0971075676 -0.1813014490 -0.2464568317; ... + 8.4775801091 1.0997316965 -0.1797082140 -0.2519152440; ... + 8.5026193458 1.1023550439 -0.1781043165 -0.2573726538; ... + 8.5262347341 1.1046324969 -0.1767049403 -0.2621078568; ... + 8.5519768095 1.1069405119 -0.1752821493 -0.2669009993; ... + 8.5761129683 1.1089683638 -0.1740297442 -0.2711051374; ... + 8.6023105873 1.1110450669 -0.1727463692 -0.2754011052; ... + 8.6267952587 1.1128866230 -0.1716087730 -0.2792006165; ... + 8.6533063618 1.1147879601 -0.1704357869 -0.2831115888; ... + 8.6780368548 1.1164862743 -0.1693902247 -0.2865931191; ... + 8.7047739081 1.1182510772 -0.1683067141 -0.2901975105; ... + 8.7316043012 1.1199566179 -0.1672632043 -0.2936663490]; + case 'P4HO2' + T = [ ... + 9.0487552928 1.0749349573 -0.2021215571 -0.1920690379; ... + 9.0396341760 1.0578442875 -0.2047569352 -0.1678532831; ... + 9.0542791209 1.0524484578 -0.2047233749 -0.1606194256; ... + 9.0601028655 1.0505487899 -0.2045610040 -0.1581773718; ... + 9.0646190673 1.0491368344 -0.2043821958 -0.1564064537; ... + 9.0707205473 1.0472961590 -0.2040685317 -0.1541623227; ... + 9.0786237393 1.0450028327 -0.2035390634 -0.1514838937; ... + 9.0888613168 1.0421535217 -0.2026443039 -0.1483728351; ... + 9.1037491359 1.0382266494 -0.2009270031 -0.1445796433; ... + 9.1306372229 1.0319435819 -0.1968531223 -0.1401578788; ... + 9.2021506431 1.0215863525 -0.1864218292 -0.1399647167; ... + 9.3830140082 1.0126738542 -0.1723454500 -0.1560887631; ... + 9.5440425474 1.0113733655 -0.1648460585 -0.1739113299; ... + 9.6990195652 1.0129644895 -0.1593270671 -0.1918297866; ... + 10.2961547121 1.0286792032 -0.1457428883 -0.2522546927; ... + 11.7792576999 1.0624281188 -0.1343562748 -0.3275531781; ... + 13.7831848966 1.0850769369 -0.1379181777 -0.3653699405; ... + 16.0024825055 1.0955445136 -0.1591111395 -0.3767847569; ... + 18.1058513583 1.0954718438 -0.1988427540 -0.3675757435; ... + 20.0406128078 1.0829659022 -0.2591580427 -0.3388198392]; + otherwise + error('OTP:InvalidInput: %s is not a valid orbit family; try L2, P2HO1, P2HO2, P4HO1, or P4HO2', ... + orbittype) +end + +end diff --git a/toolbox/+otp/+cr3bp/f.m b/toolbox/+otp/+cr3bp/f.m index e53b94c4..9f1630ac 100644 --- a/toolbox/+otp/+cr3bp/f.m +++ b/toolbox/+otp/+cr3bp/f.m @@ -7,15 +7,15 @@ r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); dddx = (x(1, :) + mu)./d; -dddy = x(2, :)./d; -dddz = x(3, :)./d; +dddy = x(2, :)./d; +dddz = x(3, :)./d; drdx = (x(1, :) - 1 + mu)./r; -drdy = x(2, :)./r; -drdz = x(3, :)./r; +drdy = x(2, :)./r; +drdz = x(3, :)./r; dUdx = x(1, :) - ((1 - mu).*dddx)./(d.^2) - (mu.*drdx)./(r.^2); dUdy = x(2, :) - ((1 - mu).*dddy)./(d.^2) - (mu.*drdy)./(r.^2); -dUdz = - ((1 - mu).*dddz)./(d.^2) - (mu.*drdz)./(r.^2); +dUdz = - ((1 - mu).*dddz)./(d.^2) - (mu.*drdz)./(r.^2); dx = v; dv = [2*v(2, :) + dUdx; -2*v(1, :) + dUdy; dUdz]; From 10cd9341c1c1732cf87743f4e1601a7156b47c29 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 1 Oct 2025 15:40:52 -1000 Subject: [PATCH 15/33] add jacobian --- toolbox/+otp/+cr3bp/CR3BPProblem.m | 4 +-- toolbox/+otp/+cr3bp/jacobian.m | 58 ++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+), 2 deletions(-) create mode 100644 toolbox/+otp/+cr3bp/jacobian.m diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index c4ba5da1..ca28ec18 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -88,10 +88,10 @@ function onSettingsChanged(obj) obj.JacobiConstant = @(y) otp.cr3bp.jacobiconstant(y, mu, soft); obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiconstantjacobian(y, mu, soft); - obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ... + 'Jacobian', @(t, y) otp.cr3bp.jacobian(t, y, mu, soft), ... 'Vectorized', 'on'); - + end function label = internalIndex2label(~, index) diff --git a/toolbox/+otp/+cr3bp/jacobian.m b/toolbox/+otp/+cr3bp/jacobian.m new file mode 100644 index 00000000..07cf9b22 --- /dev/null +++ b/toolbox/+otp/+cr3bp/jacobian.m @@ -0,0 +1,58 @@ +function J = jacobian(~, y, mu, soft) + +x = y(1:3, :); + +d = sqrt((x(1, :) + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); +r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); + +dddx = (x(1, :) + mu)./d; +dddy = x(2, :)./d; +dddz = x(3, :)./d; +drdx = (x(1, :) - 1 + mu)./r; +drdy = x(2, :)./r; +drdz = x(3, :)./r; + +dddxdx = (1 - dddx.^2)./d; +dddxdy = -(dddx.*x(2, :))./(d.^2); +dddxdz = -(dddx.*x(3, :))./(d.^2); +dddydy = (1 - dddy.^2)./d; +dddydz = -(dddy.*x(3, :))./(d.^2); +dddzdz = (1 - dddz.^2)./d; + +drdxdx = (1 - drdx.^2)./r; +drdxdy = -(drdx.*x(2, :))./(r.^2); +drdxdz = -(drdx.*x(3, :))./(r.^2); +drdydy = (1 - drdy.^2)./r; +drdydz = -(drdy.*x(3, :))./(r.^2); +drdzdz = (1 - drdz.^2)./r; + +dUdxdx = 1 + (2*(1 - mu)*dddx.^2)./(d.^3) ... + - ((1 - mu)*dddxdx)./(d.^2) ... + + (2*mu*drdx.^2)./(r.^3) ... + - (mu*drdxdx)./(r.^2); +dUdxdy = (2*(1 - mu)*dddx.*dddy)./(d.^3) ... + - ((1 - mu)*dddxdy)./(d.^2) ... + + (2*mu*drdx.*drdy)./(r.^3) ... + - (mu.*drdxdy)./(r.^2); +dUdxdz = (2*(1 - mu)*dddx.*dddz)./(d.^3) ... + - ((1 - mu)*dddxdz)./(d.^2) ... + + (2*mu*drdx.*drdz)./(r.^3) ... + - (mu.*drdxdz)./(r.^2); +dUdydy = 1 + (2*(1 - mu)*dddy.^2)./(d.^3) ... + - ((1 - mu)*dddydy)./(d.^2) ... + + (2*mu*drdy.^2)./(r.^3) ... + - (mu*drdydy)./(r.^2); +dUdydz = (2*(1 - mu)*dddy.*dddz)./(d.^3) ... + - ((1 - mu)*dddydz)./(d.^2) ... + + (2*mu*drdy.*drdz)./(r.^3) ... + - (mu.*drdydz)./(r.^2); +dUdzdz = (2*(1 - mu)*dddz.^2)./(d.^3) ... + - ((1 - mu)*dddzdz)./(d.^2) ... + + (2*mu*drdz.^2)./(r.^3) ... + - (mu*drdzdz)./(r.^2); + +dxdv = [dUdxdx, dUdxdy, dUdxdz; dUdxdy, dUdydy, dUdydz; dUdxdz, dUdydz, dUdzdz]; + +J = [zeros(3, 3), eye(3); dxdv, [0, 2, 0; -2, 0, 0; 0, 0, 0]]; + +end From 96830864a8e7e98b5c58832098db91bb01e4e3c0 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 1 Oct 2025 15:46:48 -1000 Subject: [PATCH 16/33] vectorized jacobian --- toolbox/+otp/+cr3bp/jacobian.m | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/toolbox/+otp/+cr3bp/jacobian.m b/toolbox/+otp/+cr3bp/jacobian.m index 07cf9b22..f6733adf 100644 --- a/toolbox/+otp/+cr3bp/jacobian.m +++ b/toolbox/+otp/+cr3bp/jacobian.m @@ -1,29 +1,31 @@ function J = jacobian(~, y, mu, soft) x = y(1:3, :); +x = reshape(x, 3, 1, []); +N = size(x, 3); -d = sqrt((x(1, :) + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); -r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); +d = sqrt((x(1, 1, :) + mu).^2 + x(2, 1, :).^2 + x(3, 1, :).^2 + soft^2); +r = sqrt((x(1, 1, :) - 1 + mu).^2 + x(2, 1, :).^2 + x(3, 1, :).^2 + soft^2); -dddx = (x(1, :) + mu)./d; -dddy = x(2, :)./d; -dddz = x(3, :)./d; -drdx = (x(1, :) - 1 + mu)./r; -drdy = x(2, :)./r; -drdz = x(3, :)./r; +dddx = (x(1, 1, :) + mu)./d; +dddy = x(2, 1, :)./d; +dddz = x(3, 1, :)./d; +drdx = (x(1, 1, :) - 1 + mu)./r; +drdy = x(2, 1, :)./r; +drdz = x(3, 1, :)./r; dddxdx = (1 - dddx.^2)./d; -dddxdy = -(dddx.*x(2, :))./(d.^2); -dddxdz = -(dddx.*x(3, :))./(d.^2); +dddxdy = -(dddx.*x(2, 1, :))./(d.^2); +dddxdz = -(dddx.*x(3, 1, :))./(d.^2); dddydy = (1 - dddy.^2)./d; -dddydz = -(dddy.*x(3, :))./(d.^2); +dddydz = -(dddy.*x(3, 1, :))./(d.^2); dddzdz = (1 - dddz.^2)./d; drdxdx = (1 - drdx.^2)./r; -drdxdy = -(drdx.*x(2, :))./(r.^2); -drdxdz = -(drdx.*x(3, :))./(r.^2); +drdxdy = -(drdx.*x(2, 1, :))./(r.^2); +drdxdz = -(drdx.*x(3, 1, :))./(r.^2); drdydy = (1 - drdy.^2)./r; -drdydz = -(drdy.*x(3, :))./(r.^2); +drdydz = -(drdy.*x(3, 1, :))./(r.^2); drdzdz = (1 - drdz.^2)./r; dUdxdx = 1 + (2*(1 - mu)*dddx.^2)./(d.^3) ... @@ -53,6 +55,7 @@ dxdv = [dUdxdx, dUdxdy, dUdxdz; dUdxdy, dUdydy, dUdydz; dUdxdz, dUdydz, dUdzdz]; -J = [zeros(3, 3), eye(3); dxdv, [0, 2, 0; -2, 0, 0; 0, 0, 0]]; +J = [zeros(3, 3, N), repmat(eye(3), 1, 1, N); ... + dxdv, repmat([0, 2, 0; -2, 0, 0; 0, 0, 0], 1, 1, N)]; end From bf483e4226a9348945d43c8f69ab0aa1941c569c Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Thu, 2 Oct 2025 07:51:34 -1000 Subject: [PATCH 17/33] add JVP --- toolbox/+otp/+cr3bp/CR3BPProblem.m | 3 +- toolbox/+otp/+cr3bp/jacobianVectorProduct.m | 59 +++++++++++++++++++++ 2 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 toolbox/+otp/+cr3bp/jacobianVectorProduct.m diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index ca28ec18..8e8cee41 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -90,8 +90,9 @@ function onSettingsChanged(obj) obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ... 'Jacobian', @(t, y) otp.cr3bp.jacobian(t, y, mu, soft), ... + 'JacobianVectorProduct', @(t, y, v) otp.cr3bp.jacobianVectorProduct(t, y, v, mu, soft), ... 'Vectorized', 'on'); - + end function label = internalIndex2label(~, index) diff --git a/toolbox/+otp/+cr3bp/jacobianVectorProduct.m b/toolbox/+otp/+cr3bp/jacobianVectorProduct.m new file mode 100644 index 00000000..2678509b --- /dev/null +++ b/toolbox/+otp/+cr3bp/jacobianVectorProduct.m @@ -0,0 +1,59 @@ +function Jvp = jacobianVectorProduct(~, y, v, mu, soft) + +x = y(1:3, :); + +d = sqrt((x(1, :) + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); +r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + x(3, :).^2 + soft^2); + +dddx = (x(1, :) + mu)./d; +dddy = x(2, :)./d; +dddz = x(3, :)./d; +drdx = (x(1, :) - 1 + mu)./r; +drdy = x(2, :)./r; +drdz = x(3, :)./r; + +dddxdx = (1 - dddx.^2)./d; +dddxdy = -(dddx.*x(2, :))./(d.^2); +dddxdz = -(dddx.*x(3, :))./(d.^2); +dddydy = (1 - dddy.^2)./d; +dddydz = -(dddy.*x(3, :))./(d.^2); +dddzdz = (1 - dddz.^2)./d; + +drdxdx = (1 - drdx.^2)./r; +drdxdy = -(drdx.*x(2, :))./(r.^2); +drdxdz = -(drdx.*x(3, :))./(r.^2); +drdydy = (1 - drdy.^2)./r; +drdydz = -(drdy.*x(3, :))./(r.^2); +drdzdz = (1 - drdz.^2)./r; + +dUdxdx = 1 + (2*(1 - mu)*dddx.^2)./(d.^3) ... + - ((1 - mu)*dddxdx)./(d.^2) ... + + (2*mu*drdx.^2)./(r.^3) ... + - (mu*drdxdx)./(r.^2); +dUdxdy = (2*(1 - mu)*dddx.*dddy)./(d.^3) ... + - ((1 - mu)*dddxdy)./(d.^2) ... + + (2*mu*drdx.*drdy)./(r.^3) ... + - (mu.*drdxdy)./(r.^2); +dUdxdz = (2*(1 - mu)*dddx.*dddz)./(d.^3) ... + - ((1 - mu)*dddxdz)./(d.^2) ... + + (2*mu*drdx.*drdz)./(r.^3) ... + - (mu.*drdxdz)./(r.^2); +dUdydy = 1 + (2*(1 - mu)*dddy.^2)./(d.^3) ... + - ((1 - mu)*dddydy)./(d.^2) ... + + (2*mu*drdy.^2)./(r.^3) ... + - (mu*drdydy)./(r.^2); +dUdydz = (2*(1 - mu)*dddy.*dddz)./(d.^3) ... + - ((1 - mu)*dddydz)./(d.^2) ... + + (2*mu*drdy.*drdz)./(r.^3) ... + - (mu.*drdydz)./(r.^2); +dUdzdz = (2*(1 - mu)*dddz.^2)./(d.^3) ... + - ((1 - mu)*dddzdz)./(d.^2) ... + + (2*mu*drdz.^2)./(r.^3) ... + - (mu*drdzdz)./(r.^2); + +Jvpv = [dUdxdx.*v(1, :) + dUdxdy.*v(2, :) + dUdxdz.*v(3, :) + 2*v(5, :); ... + dUdxdy.*v(1, :) + dUdydy.*v(2, :) + dUdydz.*v(3, :) - 2*v(4, :); ... + dUdxdz.*v(1, :) + dUdydz.*v(2, :) + dUdzdz.*v(3, :)]; +Jvp = [v(4:6, :); Jvpv]; + +end From fd2e1c59a3341552821e430eb88cbc2583d393a5 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Mon, 6 Oct 2025 04:34:24 -1000 Subject: [PATCH 18/33] added LEO preset --- .../+otp/+cr3bp/+presets/+utils/buildOrbit.m | 54 +++++++++++++++++++ toolbox/+otp/+cr3bp/+presets/LowEarthOrbit.m | 31 +++++++++++ toolbox/+otp/+cr3bp/CR3BPProblem.m | 3 +- 3 files changed, 87 insertions(+), 1 deletion(-) create mode 100644 toolbox/+otp/+cr3bp/+presets/+utils/buildOrbit.m create mode 100644 toolbox/+otp/+cr3bp/+presets/LowEarthOrbit.m diff --git a/toolbox/+otp/+cr3bp/+presets/+utils/buildOrbit.m b/toolbox/+otp/+cr3bp/+presets/+utils/buildOrbit.m new file mode 100644 index 00000000..05bf9521 --- /dev/null +++ b/toolbox/+otp/+cr3bp/+presets/+utils/buildOrbit.m @@ -0,0 +1,54 @@ +problem = otp.cr3bp.presets.LowEarthOrbit; + +y0 = problem.Y0; +tend = problem.TimeSpan(2); +f = problem.RHS.F; +mu = problem.Parameters.Mu; + +v0 = -7; + +J = @(v0) cost(v0, y0, f, tend, mu); +opts = optimoptions("fmincon", "Algorithm","interior-point", ... + "FunctionTolerance", 1e-14, "ConstraintTolerance", 1e-14, 'Display','off'); +[v0, fev] = fmincon(J, v0, [], [], [], [], [], [], [], opts); + +v0 + +% find period + +J = @(tend) cost2(v0, y0, f, tend, mu); +opts = optimoptions("fmincon", "Algorithm","interior-point", ... + "FunctionTolerance", 1e-14, "ConstraintTolerance", 1e-14, 'Display','off'); +[tend, fev] = fmincon(J, tend, [], [], [], [], 0.13, 0.15, [], opts); + +tend + + +function c = cost(v0, y0, f, tend, mu) + +y0(5) = v0; + +sol = ode45(f, [0, tend], y0, ... + odeset('AbsTol', 1e-14, 'RelTol', 100*eps)); + +yend = sol.y(:, :); + +earth = [-mu; 0; 0]; + +dy0 = norm(y0(1:3) - earth); +dyend = vecnorm(yend(1:3, :) - earth); + +c = sum((dy0 - dyend).^2); + +end + +function c = cost2(v0, y0, f, tend, ~) + +y0(5) = v0; + +sol = ode45(f, [0, tend], y0, ... + odeset('AbsTol', 1e-14, 'RelTol', 100*eps)); + +c = sum((y0 - sol.y(:, end)).^2); + +end diff --git a/toolbox/+otp/+cr3bp/+presets/LowEarthOrbit.m b/toolbox/+otp/+cr3bp/+presets/LowEarthOrbit.m new file mode 100644 index 00000000..26019da9 --- /dev/null +++ b/toolbox/+otp/+cr3bp/+presets/LowEarthOrbit.m @@ -0,0 +1,31 @@ +classdef LowEarthOrbit < otp.cr3bp.CR3BPProblem + % A trivial preset with a stable oscillating orbit around a lagrange point. + methods + function obj = LowEarthOrbit(varargin) + % Create the Canonical CR3BP problem object. + + + equatorialRadius = 6378; + earthMoonDist = 385000; + + delta = equatorialRadius/earthMoonDist; + + mE = otp.utils.PhysicalConstants.EarthMass; + mL = otp.utils.PhysicalConstants.MoonMass; + G = otp.utils.PhysicalConstants.GravitationalConstant; + + % derive mu + muE = G*mE; + muL = G*mL; + mu = muL/(muE + muL); + + x0 = -mu + delta; + y0 = -7.717617954315369e+00; + + y0 = [x0; 0; 0; 0; y0; 0]; + tspan = [0, 1.348715080895850e-01]; + params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 1e-3, varargin{:}); + obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); + end + end +end diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index 8e8cee41..e4d6073d 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -125,7 +125,8 @@ function onSettingsChanged(obj) fig = internalPlotPhaseSpace@otp.Problem(obj, t, y, ... 'Vars', 1:3, varargin{:}); hold on; - scatter3([mu, 1 - mu], [0, 0], [0, 0]); + scatter3([-mu, 1 - mu], [0, 0], [0, 0]); + axis equal end end From dfba98d0128f15d80cbdb7649f70f4f1b2d4c32c Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 7 Oct 2025 07:55:53 -1000 Subject: [PATCH 19/33] Replace LEO with CircularEarthOrbit --- docs/references.bib | 8 +++ .../+otp/+cr3bp/+presets/CircularEarthOrbit.m | 49 +++++++++++++++++++ toolbox/+otp/+cr3bp/+presets/LowEarthOrbit.m | 31 ------------ 3 files changed, 57 insertions(+), 31 deletions(-) create mode 100644 toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m delete mode 100644 toolbox/+otp/+cr3bp/+presets/LowEarthOrbit.m diff --git a/docs/references.bib b/docs/references.bib index f460e9a0..b60dac78 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -217,3 +217,11 @@ @phdthesis{Spr21 school={Purdue University} } +@book{She20, + title={A Preliminary Study of Leo to Geo Transfers for Inclination Changes Using Libration Point Orbits}, + author={Shepard, John Philip}, + year={2020}, + publisher={The University of North Dakota} +} + + diff --git a/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m b/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m new file mode 100644 index 00000000..985d2f5a --- /dev/null +++ b/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m @@ -0,0 +1,49 @@ +classdef CircularEarthOrbit < otp.cr3bp.CR3BPProblem + % This preset builds a circular earth orbit in CR3BP based on + % the equations derived in :cite:p:She20. + + methods + function obj = CircularEarthOrbit(varargin) + % Create the CircularEarthOrbit CR3BP problem object. + % + % Parameters + % ---------- + % OrbitalRadius : numeric(1, 1) + % The radius of the orbit above Earth's surface in km. + + p = inputParser(); + p.addParameter('OrbitalRadius', 6556, @isnumeric); + p.parse(varargin{:}); + results = p.Results; + + equatorialRadius = 6378; + earthMoonDist = 385000; + orbitalradius = results.OrbitalRadius; + + delta = (equatorialRadius + orbitalradius)/earthMoonDist; + + mE = otp.utils.PhysicalConstants.EarthMass; + mL = otp.utils.PhysicalConstants.MoonMass; + G = otp.utils.PhysicalConstants.GravitationalConstant; + + % derive mu + muE = G*mE; + muL = G*mL; + mu = muL/(muE + muL); + + % set initial distance to Earth + x0 = -mu + delta; + + % derive the initial velocity + y0 = -sqrt((1 - mu)/delta); + + % derive the orbital period + period = sqrt(((abs(delta))^3)/(1 - mu))*2*pi; + + y0 = [x0; 0; 0; 0; y0; 0]; + tspan = [0, period]; + params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 1e-3, varargin{:}); + obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); + end + end +end diff --git a/toolbox/+otp/+cr3bp/+presets/LowEarthOrbit.m b/toolbox/+otp/+cr3bp/+presets/LowEarthOrbit.m deleted file mode 100644 index 26019da9..00000000 --- a/toolbox/+otp/+cr3bp/+presets/LowEarthOrbit.m +++ /dev/null @@ -1,31 +0,0 @@ -classdef LowEarthOrbit < otp.cr3bp.CR3BPProblem - % A trivial preset with a stable oscillating orbit around a lagrange point. - methods - function obj = LowEarthOrbit(varargin) - % Create the Canonical CR3BP problem object. - - - equatorialRadius = 6378; - earthMoonDist = 385000; - - delta = equatorialRadius/earthMoonDist; - - mE = otp.utils.PhysicalConstants.EarthMass; - mL = otp.utils.PhysicalConstants.MoonMass; - G = otp.utils.PhysicalConstants.GravitationalConstant; - - % derive mu - muE = G*mE; - muL = G*mL; - mu = muL/(muE + muL); - - x0 = -mu + delta; - y0 = -7.717617954315369e+00; - - y0 = [x0; 0; 0; 0; y0; 0]; - tspan = [0, 1.348715080895850e-01]; - params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 1e-3, varargin{:}); - obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); - end - end -end From 5cd3e39995a55fc6451180880e6dc5597fcdf0aa Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 7 Oct 2025 08:29:47 -1000 Subject: [PATCH 20/33] deleted buildrobit script, but added a blacklist to codespell --- docs/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Makefile b/docs/Makefile index f108913d..0e340dea 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -19,7 +19,7 @@ serve: --ignore $(PROBLEMSDIR) lint: - pycodestyle .. && codespell .. + pycodestyle .. && codespell --skip="*.bib" .. clean: rm -rf $(BUILDDIR) $(PROBLEMSDIR) From 1790115e220ff3448f9c11b12294c2506b4d32fd Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 7 Oct 2025 11:56:26 -1000 Subject: [PATCH 21/33] fixed CircularEarthOrbit script --- toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m b/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m index 985d2f5a..a99dec0b 100644 --- a/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m +++ b/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m @@ -12,9 +12,11 @@ % The radius of the orbit above Earth's surface in km. p = inputParser(); - p.addParameter('OrbitalRadius', 6556, @isnumeric); + p.KeepUnmatched = true; + p.addParameter('OrbitalRadius', 340, @isnumeric); p.parse(varargin{:}); results = p.Results; + varargin = [fieldnames(p.Unmatched), struct2cell(p.Unmatched)].'; equatorialRadius = 6378; earthMoonDist = 385000; From 50e2974eba99d8bcbd6479d4ec0efedae1daa1f9 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 7 Oct 2025 15:23:27 -1000 Subject: [PATCH 22/33] add radar measurement --- .../+otp/+cr3bp/+presets/CircularEarthOrbit.m | 5 ++++- toolbox/+otp/+cr3bp/CR3BPProblem.m | 4 +++- toolbox/+otp/+cr3bp/radarMeasurement.m | 16 ++++++++++++++++ 3 files changed, 23 insertions(+), 2 deletions(-) create mode 100644 toolbox/+otp/+cr3bp/radarMeasurement.m diff --git a/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m b/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m index a99dec0b..8a98e15a 100644 --- a/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m +++ b/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m @@ -10,10 +10,12 @@ % ---------- % OrbitalRadius : numeric(1, 1) % The radius of the orbit above Earth's surface in km. + % The default radius is 340 km, which is approximately where + % the ISS is/was located. p = inputParser(); p.KeepUnmatched = true; - p.addParameter('OrbitalRadius', 340, @isnumeric); + p.addParameter('OrbitalRadius', 340, @(x) isnumeric(x) && isscalar(x)); p.parse(varargin{:}); results = p.Results; varargin = [fieldnames(p.Unmatched), struct2cell(p.Unmatched)].'; @@ -22,6 +24,7 @@ earthMoonDist = 385000; orbitalradius = results.OrbitalRadius; + % set the relative distance from earth to the moon delta = (equatorialRadius + orbitalradius)/earthMoonDist; mE = otp.utils.PhysicalConstants.EarthMass; diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index e4d6073d..6352d1a6 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -77,7 +77,7 @@ properties (SetAccess = private) JacobiConstant JacobiConstantJacobian - RHSStab + RadarMeasurement end methods (Access = protected) @@ -88,6 +88,8 @@ function onSettingsChanged(obj) obj.JacobiConstant = @(y) otp.cr3bp.jacobiconstant(y, mu, soft); obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiconstantjacobian(y, mu, soft); + obj.RadarMeasurement = @(y, radary) otp.cr3bp.radarMeasurement(t, y, mu, soft, radary); + obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ... 'Jacobian', @(t, y) otp.cr3bp.jacobian(t, y, mu, soft), ... 'JacobianVectorProduct', @(t, y, v) otp.cr3bp.jacobianVectorProduct(t, y, v, mu, soft), ... diff --git a/toolbox/+otp/+cr3bp/radarMeasurement.m b/toolbox/+otp/+cr3bp/radarMeasurement.m new file mode 100644 index 00000000..4cbbdd72 --- /dev/null +++ b/toolbox/+otp/+cr3bp/radarMeasurement.m @@ -0,0 +1,16 @@ +function rmeas = radarMeasurement(~, y, ~, ~, radary) + +x = y(1:3, :); +v = y(4:6, :); + +relpos = x - radary; + +r = sqrt(relpos(1, :).^2 + relpos(2, :).^2 + relpos(3, :).^2); +drdt = sum(relpos.*v, 1)/r; +theta = atan2(relpos(2, :), relpos(1, :)); +phi = atan2(relpos(3, :), r); + + +rmeas = [r; drdt; theta; phi]; + +end From e609e5f2e3602351663161e9347b05ef66f986fb Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 7 Oct 2025 15:26:28 -1000 Subject: [PATCH 23/33] renamed some files --- toolbox/+otp/+cr3bp/CR3BPProblem.m | 4 ++-- toolbox/+otp/+cr3bp/{jacobiconstant.m => jacobiConstant.m} | 0 .../{jacobiconstantjacobian.m => jacobiConstantJacobian.m} | 0 3 files changed, 2 insertions(+), 2 deletions(-) rename toolbox/+otp/+cr3bp/{jacobiconstant.m => jacobiConstant.m} (100%) rename toolbox/+otp/+cr3bp/{jacobiconstantjacobian.m => jacobiConstantJacobian.m} (100%) diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index 6352d1a6..a919b323 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -85,8 +85,8 @@ function onSettingsChanged(obj) mu = obj.Parameters.Mu; soft = obj.Parameters.SoftFactor; - obj.JacobiConstant = @(y) otp.cr3bp.jacobiconstant(y, mu, soft); - obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiconstantjacobian(y, mu, soft); + obj.JacobiConstant = @(y) otp.cr3bp.jacobiConstant(y, mu, soft); + obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiConstantJacobian(y, mu, soft); obj.RadarMeasurement = @(y, radary) otp.cr3bp.radarMeasurement(t, y, mu, soft, radary); diff --git a/toolbox/+otp/+cr3bp/jacobiconstant.m b/toolbox/+otp/+cr3bp/jacobiConstant.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobiconstant.m rename to toolbox/+otp/+cr3bp/jacobiConstant.m diff --git a/toolbox/+otp/+cr3bp/jacobiconstantjacobian.m b/toolbox/+otp/+cr3bp/jacobiConstantJacobian.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobiconstantjacobian.m rename to toolbox/+otp/+cr3bp/jacobiConstantJacobian.m From 99cc582c2a97d02451342ca22306f75ea4975885 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 7 Oct 2025 16:11:20 -1000 Subject: [PATCH 24/33] temp --- .../{jacobiConstantJacobian.m => jacobiConstantJacobian2.m} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename toolbox/+otp/+cr3bp/{jacobiConstantJacobian.m => jacobiConstantJacobian2.m} (100%) diff --git a/toolbox/+otp/+cr3bp/jacobiConstantJacobian.m b/toolbox/+otp/+cr3bp/jacobiConstantJacobian2.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobiConstantJacobian.m rename to toolbox/+otp/+cr3bp/jacobiConstantJacobian2.m From adbb3d8d21af757bc1ed5d484666ac66df530589 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 7 Oct 2025 16:11:37 -1000 Subject: [PATCH 25/33] temp2 --- .../{jacobiConstantJacobian2.m => jacobiConstantJacobian.m} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename toolbox/+otp/+cr3bp/{jacobiConstantJacobian2.m => jacobiConstantJacobian.m} (100%) diff --git a/toolbox/+otp/+cr3bp/jacobiConstantJacobian2.m b/toolbox/+otp/+cr3bp/jacobiConstantJacobian.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobiConstantJacobian2.m rename to toolbox/+otp/+cr3bp/jacobiConstantJacobian.m From 9a8d016f8232c7723bdac07d0c0bdeb3e6d357b1 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Tue, 7 Oct 2025 16:14:47 -1000 Subject: [PATCH 26/33] potential fix? --- toolbox/+otp/+cr3bp/jacobiConstant.m | 2 +- toolbox/+otp/+cr3bp/jacobiConstantJacobian.m | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/toolbox/+otp/+cr3bp/jacobiConstant.m b/toolbox/+otp/+cr3bp/jacobiConstant.m index 18cc2bbd..054a3245 100644 --- a/toolbox/+otp/+cr3bp/jacobiConstant.m +++ b/toolbox/+otp/+cr3bp/jacobiConstant.m @@ -1,4 +1,4 @@ -function J = jacobiconstant(y, mu, soft) +function J = jacobiConstant(y, mu, soft) x = y(1:3, :); v = y(4:6, :); diff --git a/toolbox/+otp/+cr3bp/jacobiConstantJacobian.m b/toolbox/+otp/+cr3bp/jacobiConstantJacobian.m index eb19fe9b..8a114554 100644 --- a/toolbox/+otp/+cr3bp/jacobiConstantJacobian.m +++ b/toolbox/+otp/+cr3bp/jacobiConstantJacobian.m @@ -1,4 +1,4 @@ -function dJdy = jacobiconstantjacobian(y, mu, soft) +function dJdy = jacobiConstantJacobian(y, mu, soft) x = y(1:3, :); v = y(4:6, :); From 64424beedc9e2d3947c0a7c248e91a6dcc3ece01 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 8 Oct 2025 08:32:33 -1000 Subject: [PATCH 27/33] Initial commit of Arenstorf problem into cr3bp. --- toolbox/+otp/+cr3bp/+presets/Arenstorf.m | 24 ++++ toolbox/+otp/+cr3bp/+presets/Canonical.m | 5 +- toolbox/+otp/+cr3bp/CR3BPProblem.m | 106 ++++++++++++------ toolbox/+otp/+cr3bp/fPlanar.m | 22 ++++ toolbox/+otp/+cr3bp/jacobianPlanar.m | 41 +++++++ .../+otp/+cr3bp/jacobianVectorProductPlanar.m | 38 +++++++ 6 files changed, 203 insertions(+), 33 deletions(-) create mode 100644 toolbox/+otp/+cr3bp/+presets/Arenstorf.m create mode 100644 toolbox/+otp/+cr3bp/fPlanar.m create mode 100644 toolbox/+otp/+cr3bp/jacobianPlanar.m create mode 100644 toolbox/+otp/+cr3bp/jacobianVectorProductPlanar.m diff --git a/toolbox/+otp/+cr3bp/+presets/Arenstorf.m b/toolbox/+otp/+cr3bp/+presets/Arenstorf.m new file mode 100644 index 00000000..4f12c362 --- /dev/null +++ b/toolbox/+otp/+cr3bp/+presets/Arenstorf.m @@ -0,0 +1,24 @@ +classdef Arenstorf < otp.cr3bp.CR3BPProblem + % One period of a satellite moving in an Earth-Moon system on a planar + % orbit. See pages 129--130 of :cite:p:`HNW93` for more details. + + methods + function obj = Arenstorf(varargin) + % Create the Arenstorf CR3BP problem object. + + params = otp.cr3bp.CR3BPParameters('Mu', 0.012277471, 'SoftFactor', 0, varargin{:}); + cls = class(params.Mu); + + % Decimals converted to rational to support multiple data types + y0 = [cast(497, cls) / cast(500, cls); ... + cast(0, cls); ... + cast(0, cls); ... + cast(-823970832321143, cls) / cast(411659154384760, cls)]; + tspan = [cast(0, cls); ... + cast(4541277234950502, cls) / cast(266113073862361, cls)]; + + obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); + end + + end +end diff --git a/toolbox/+otp/+cr3bp/+presets/Canonical.m b/toolbox/+otp/+cr3bp/+presets/Canonical.m index b374efc7..a84d72b3 100644 --- a/toolbox/+otp/+cr3bp/+presets/Canonical.m +++ b/toolbox/+otp/+cr3bp/+presets/Canonical.m @@ -1,5 +1,6 @@ classdef Canonical < otp.cr3bp.CR3BPProblem % A trivial preset with a stable oscillating orbit around a lagrange point. + methods function obj = Canonical(varargin) % Create the Canonical CR3BP problem object. @@ -8,7 +9,9 @@ y0 = [0; 0; 0; 0; 0; 1]; tspan = [0, 10]; - params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 1e-3, varargin{:}); + params = otp.cr3bp.CR3BPParameters('Mu', mu, ... + 'SoftFactor', 1e-3, ... + varargin{:}); obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); end end diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index a919b323..f2ef452a 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -16,17 +16,26 @@ % $$ % x'' &= 2y' + \frac{\partial U}{\partial x},\\ % y'' &= -2x' + \frac{\partial U}{\partial y},\\ - % z'' &= \frac{\partial U}{\partial z},\\ + % z'' &= \frac{\partial U}{\partial z}, + % $$ + % + % where + % + % $$ % U &= \frac{1}{2} (x^2 + y^2) + \frac{1 - \mu}{d} + \frac{mu}{r},\\ % d &= \sqrt{(x + \mu)^2 + y^2 + z^2},\\ % r &= \sqrt{(x - 1 + \mu)^2 + y^2 + z^2}, % $$ % - % where the system is converted to a differential equation in six + % and where the system is converted to a differential equation in six % variables in the standard fashion. The distances $d$ and $r$ can % cause numerical instability as they approach zero, thus a softening % factor of $s^2$ is typically added under both of the square-roots. % + % When the object under consideration is on an orbit that is co-planar + % to the orbit of the two other objects, then the system of equations + % can reduce by two dimensions, removing the $z''$ and $z'$ terms. + % % The system of equations is energy-preserving, meaning that the Jacobi % constant of the system, % @@ -42,7 +51,7 @@ % +---------------------+-----------------------------------------------+ % | Type | ODE | % +---------------------+-----------------------------------------------+ - % | Number of Variables | 6 | + % | Number of Variables | 4 for planar or 6 for non-planar | % +---------------------+-----------------------------------------------+ % | Stiff | no | % +---------------------+-----------------------------------------------+ @@ -65,12 +74,12 @@ % ---------- % timeSpan : numeric(1, 2) % The start and final time. - % y0 : numeric(6, 1) - % The initial conditions. + % y0 : numeric(n, 1) + % The initial conditions. Either n = 4 or n = 6. % parameters : otp.cr3bp.CR3BPParameters % The parameters. - obj@otp.Problem('Circular Restricted 3 Body Problem', 6, timeSpan, y0, parameters); + obj@otp.Problem('Circular Restricted 3 Body Problem', [], timeSpan, y0, parameters); end end @@ -82,35 +91,59 @@ methods (Access = protected) function onSettingsChanged(obj) - mu = obj.Parameters.Mu; - soft = obj.Parameters.SoftFactor; + mu = obj.Parameters.Mu; + soft = obj.Parameters.SoftFactor; obj.JacobiConstant = @(y) otp.cr3bp.jacobiConstant(y, mu, soft); obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiConstantJacobian(y, mu, soft); obj.RadarMeasurement = @(y, radary) otp.cr3bp.radarMeasurement(t, y, mu, soft, radary); - obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ... - 'Jacobian', @(t, y) otp.cr3bp.jacobian(t, y, mu, soft), ... - 'JacobianVectorProduct', @(t, y, v) otp.cr3bp.jacobianVectorProduct(t, y, v, mu, soft), ... - 'Vectorized', 'on'); + spatialdim = numel(obj.Y0)/2; + if ~(spatialdim == 2 || spatialdim == 3) + error('TODO') + end + if spatialdim == 3 + obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ... + 'Jacobian', @(t, y) otp.cr3bp.jacobian(t, y, mu, soft), ... + 'JacobianVectorProduct', @(t, y, v) otp.cr3bp.jacobianVectorProduct(t, y, v, mu, soft), ... + 'Vectorized', 'on'); + else + obj.RHS = otp.RHS(@(t, y) otp.cr3bp.fPlanar(t, y, mu, soft), ... + 'Jacobian', @(t, y) otp.cr3bp.jacobianPlanar(t, y, mu, soft), ... + 'JacobianVectorProduct', @(t, y, v) otp.cr3bp.jacobianVectorProductPlanar(t, y, v, mu, soft), ... + 'Vectorized', 'on'); + end end - - function label = internalIndex2label(~, index) - switch index - case 1 - label = 'x'; - case 2 - label = 'y'; - case 3 - label = 'z'; - case 4 - label = 'dx'; - case 5 - label = 'dy'; - case 6 - label = 'dz'; + + function label = internalIndex2label(obj, index) + if numel(obj.Y0) == 6 + switch index + case 1 + label = 'x'; + case 2 + label = 'y'; + case 3 + label = 'z'; + case 4 + label = 'dx'; + case 5 + label = 'dy'; + case 6 + label = 'dz'; + end + elseif numel(obj.Y0) == 4 + switch index + case 1 + label = 'x'; + case 2 + label = 'y'; + case 3 + label = 'dx'; + case 4 + label = 'dy'; + end end end @@ -124,11 +157,20 @@ function onSettingsChanged(obj) function fig = internalPlotPhaseSpace(obj, t, y, varargin) mu = obj.Parameters.Mu; - fig = internalPlotPhaseSpace@otp.Problem(obj, t, y, ... - 'Vars', 1:3, varargin{:}); - hold on; - scatter3([-mu, 1 - mu], [0, 0], [0, 0]); - axis equal + + if numel(obj.Y0) == 6 + fig = internalPlotPhaseSpace@otp.Problem(obj, t, y, ... + 'Vars', 1:3, varargin{:}); + hold on; + scatter3([-mu, 1 - mu], [0, 0], [0, 0]); + axis equal + elseif numel(obj.Y0) == 4 + fig = internalPlotPhaseSpace@otp.Problem(obj, t, y, ... + 'Vars', 1:2, varargin{:}); + hold on; + scatter([-mu, 1 - mu], [0, 0]); + axis equal + end end end diff --git a/toolbox/+otp/+cr3bp/fPlanar.m b/toolbox/+otp/+cr3bp/fPlanar.m new file mode 100644 index 00000000..c4efdfe3 --- /dev/null +++ b/toolbox/+otp/+cr3bp/fPlanar.m @@ -0,0 +1,22 @@ +function dy = fPlanar(~, y, mu, soft) + +x = y(1:2, :); +v = y(3:4, :); + +d = sqrt((x(1, :) + mu).^2 + x(2, :).^2 + soft^2); +r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + soft^2); + +dddx = (x(1, :) + mu)./d; +dddy = x(2, :)./d; +drdx = (x(1, :) - 1 + mu)./r; +drdy = x(2, :)./r; + +dUdx = x(1, :) - ((1 - mu).*dddx)./(d.^2) - (mu.*drdx)./(r.^2); +dUdy = x(2, :) - ((1 - mu).*dddy)./(d.^2) - (mu.*drdy)./(r.^2); + +dx = v; +dv = [2*v(2, :) + dUdx; -2*v(1, :) + dUdy]; + +dy = [dx; dv]; + +end diff --git a/toolbox/+otp/+cr3bp/jacobianPlanar.m b/toolbox/+otp/+cr3bp/jacobianPlanar.m new file mode 100644 index 00000000..945363e2 --- /dev/null +++ b/toolbox/+otp/+cr3bp/jacobianPlanar.m @@ -0,0 +1,41 @@ +function J = jacobianPlanar(~, y, mu, soft) + +x = y(1:2, :); +x = reshape(x, 2, 1, []); +N = size(x, 3); + +d = sqrt((x(1, 1, :) + mu).^2 + x(2, 1, :).^2 + soft^2); +r = sqrt((x(1, 1, :) - 1 + mu).^2 + x(2, 1, :).^2 + soft^2); + +dddx = (x(1, 1, :) + mu)./d; +dddy = x(2, 1, :)./d; +drdx = (x(1, 1, :) - 1 + mu)./r; +drdy = x(2, 1, :)./r; + +dddxdx = (1 - dddx.^2)./d; +dddxdy = -(dddx.*x(2, 1, :))./(d.^2); +dddydy = (1 - dddy.^2)./d; + +drdxdx = (1 - drdx.^2)./r; +drdxdy = -(drdx.*x(2, 1, :))./(r.^2); +drdydy = (1 - drdy.^2)./r; + +dUdxdx = 1 + (2*(1 - mu)*dddx.^2)./(d.^3) ... + - ((1 - mu)*dddxdx)./(d.^2) ... + + (2*mu*drdx.^2)./(r.^3) ... + - (mu*drdxdx)./(r.^2); +dUdxdy = (2*(1 - mu)*dddx.*dddy)./(d.^3) ... + - ((1 - mu)*dddxdy)./(d.^2) ... + + (2*mu*drdx.*drdy)./(r.^3) ... + - (mu.*drdxdy)./(r.^2); +dUdydy = 1 + (2*(1 - mu)*dddy.^2)./(d.^3) ... + - ((1 - mu)*dddydy)./(d.^2) ... + + (2*mu*drdy.^2)./(r.^3) ... + - (mu*drdydy)./(r.^2); + +dxdv = [dUdxdx, dUdxdy; dUdxdy, dUdydy]; + +J = [zeros(2, 2, N), repmat(eye(2), 1, 1, N); ... + dxdv, repmat([0, 2; -2, 0], 1, 1, N)]; + +end diff --git a/toolbox/+otp/+cr3bp/jacobianVectorProductPlanar.m b/toolbox/+otp/+cr3bp/jacobianVectorProductPlanar.m new file mode 100644 index 00000000..ecbf95be --- /dev/null +++ b/toolbox/+otp/+cr3bp/jacobianVectorProductPlanar.m @@ -0,0 +1,38 @@ +function Jvp = jacobianVectorProductPlanar(~, y, v, mu, soft) + +x = y(1:2, :); + +d = sqrt((x(1, :) + mu).^2 + x(2, :).^2 + soft^2); +r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + soft^2); + +dddx = (x(1, :) + mu)./d; +dddy = x(2, :)./d; +drdx = (x(1, :) - 1 + mu)./r; +drdy = x(2, :)./r; + +dddxdx = (1 - dddx.^2)./d; +dddxdy = -(dddx.*x(2, :))./(d.^2); +dddydy = (1 - dddy.^2)./d; + +drdxdx = (1 - drdx.^2)./r; +drdxdy = -(drdx.*x(2, :))./(r.^2); +drdydy = (1 - drdy.^2)./r; + +dUdxdx = 1 + (2*(1 - mu)*dddx.^2)./(d.^3) ... + - ((1 - mu)*dddxdx)./(d.^2) ... + + (2*mu*drdx.^2)./(r.^3) ... + - (mu*drdxdx)./(r.^2); +dUdxdy = (2*(1 - mu)*dddx.*dddy)./(d.^3) ... + - ((1 - mu)*dddxdy)./(d.^2) ... + + (2*mu*drdx.*drdy)./(r.^3) ... + - (mu.*drdxdy)./(r.^2); +dUdydy = 1 + (2*(1 - mu)*dddy.^2)./(d.^3) ... + - ((1 - mu)*dddydy)./(d.^2) ... + + (2*mu*drdy.^2)./(r.^3) ... + - (mu*drdydy)./(r.^2); + +Jvpv = [dUdxdx.*v(1, :) + dUdxdy.*v(2, :) + 2*v(4, :); ... + dUdxdy.*v(1, :) + dUdydy.*v(2, :) - 2*v(3, :)]; +Jvp = [v(3:4, :); Jvpv]; + +end From 7a1f4cc21fdc9bed906ae7fd5624b85427141ca2 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 8 Oct 2025 09:38:33 -1000 Subject: [PATCH 28/33] more planar additions --- toolbox/+otp/+cr3bp/CR3BPProblem.m | 35 ++++++++++++------- .../+cr3bp/jacobiConstantJacobianPlanar.m | 22 ++++++++++++ toolbox/+otp/+cr3bp/jacobiConstantPlanar.m | 13 +++++++ toolbox/+otp/+cr3bp/radarMeasurementPlanar.m | 14 ++++++++ 4 files changed, 72 insertions(+), 12 deletions(-) create mode 100644 toolbox/+otp/+cr3bp/jacobiConstantJacobianPlanar.m create mode 100644 toolbox/+otp/+cr3bp/jacobiConstantPlanar.m create mode 100644 toolbox/+otp/+cr3bp/radarMeasurementPlanar.m diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+cr3bp/CR3BPProblem.m index f2ef452a..048b9e8d 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+cr3bp/CR3BPProblem.m @@ -94,26 +94,37 @@ function onSettingsChanged(obj) mu = obj.Parameters.Mu; soft = obj.Parameters.SoftFactor; - obj.JacobiConstant = @(y) otp.cr3bp.jacobiConstant(y, mu, soft); - obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiConstantJacobian(y, mu, soft); - - obj.RadarMeasurement = @(y, radary) otp.cr3bp.radarMeasurement(t, y, mu, soft, radary); - spatialdim = numel(obj.Y0)/2; if ~(spatialdim == 2 || spatialdim == 3) - error('TODO') + error('OTP:invalidY0', 'Expected Y0 to have 4 or 6 components but has %d', numel(obj.Y0)); end if spatialdim == 3 + obj.JacobiConstant = @(y) otp.cr3bp.jacobiConstant(y, mu, soft); + obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiConstantJacobian(y, mu, soft); + + obj.RadarMeasurement = @(y, radary) otp.cr3bp.radarMeasurement(t, y, mu, soft, radary); + obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ... - 'Jacobian', @(t, y) otp.cr3bp.jacobian(t, y, mu, soft), ... - 'JacobianVectorProduct', @(t, y, v) otp.cr3bp.jacobianVectorProduct(t, y, v, mu, soft), ... - 'Vectorized', 'on'); + 'Jacobian', ... + @(t, y) otp.cr3bp.jacobian(t, y, mu, soft), ... + 'JacobianVectorProduct', ... + @(t, y, v) otp.cr3bp.jacobianVectorProduct(t, y, v, mu, soft), ... + 'Vectorized', ... + 'on'); else + obj.JacobiConstant = @(y) otp.cr3bp.jacobiConstantPlanar(y, mu, soft); + obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiConstantJacobianPlanar(y, mu, soft); + + obj.RadarMeasurement = @(y, radary) otp.cr3bp.radarMeasurementPlanar(t, y, mu, soft, radary); + obj.RHS = otp.RHS(@(t, y) otp.cr3bp.fPlanar(t, y, mu, soft), ... - 'Jacobian', @(t, y) otp.cr3bp.jacobianPlanar(t, y, mu, soft), ... - 'JacobianVectorProduct', @(t, y, v) otp.cr3bp.jacobianVectorProductPlanar(t, y, v, mu, soft), ... - 'Vectorized', 'on'); + 'Jacobian', ... + @(t, y) otp.cr3bp.jacobianPlanar(t, y, mu, soft), ... + 'JacobianVectorProduct', ... + @(t, y, v) otp.cr3bp.jacobianVectorProductPlanar(t, y, v, mu, soft), ... + 'Vectorized', ... + 'on'); end end diff --git a/toolbox/+otp/+cr3bp/jacobiConstantJacobianPlanar.m b/toolbox/+otp/+cr3bp/jacobiConstantJacobianPlanar.m new file mode 100644 index 00000000..f5210966 --- /dev/null +++ b/toolbox/+otp/+cr3bp/jacobiConstantJacobianPlanar.m @@ -0,0 +1,22 @@ +function dJdy = jacobiConstantJacobianPlanar(y, mu, soft) + +x = y(1:2, :); +v = y(3:4, :); + +d = sqrt((x(1, :) + mu).^2 + x(2, :).^2 + soft^2); +r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + soft^2); + +dddx = (x(1, :) + mu)./d; +dddy = x(2, :)./d; +drdx = (x(1, :) - 1 + mu)./r; +drdy = x(2, :)./r; + +dUdx = x(1, :) - ((1 - mu).*dddx)./(d.^2) - (mu.*drdx)./(r.^2); +dUdy = x(2, :) - ((1 - mu).*dddy)./(d.^2) - (mu.*drdy)./(r.^2); + +dJdy = [reshape(2*dUdx, 1, 1, []), ... + reshape(2*dUdy, 1, 1, []), ... + reshape(-2*v(1, :), 1, 1, []), ... + reshape(-2*v(2, :), 1, 1, [])]; + +end diff --git a/toolbox/+otp/+cr3bp/jacobiConstantPlanar.m b/toolbox/+otp/+cr3bp/jacobiConstantPlanar.m new file mode 100644 index 00000000..8068ae32 --- /dev/null +++ b/toolbox/+otp/+cr3bp/jacobiConstantPlanar.m @@ -0,0 +1,13 @@ +function J = jacobiConstantPlanar(y, mu, soft) + +x = y(1:2, :); +v = y(3:4, :); + +d = sqrt((x(1, :) + mu).^2 + x(2, :).^2 + soft^2); +r = sqrt((x(1, :) - 1 + mu).^2 + x(2, :).^2 + soft^2); + +U = 0.5*(x(1, :).^2 + x(2, :).^2) + (1 - mu)./d + mu./r; + +J = 2*U - v(1, :).^2 - v(2, :).^2; + +end diff --git a/toolbox/+otp/+cr3bp/radarMeasurementPlanar.m b/toolbox/+otp/+cr3bp/radarMeasurementPlanar.m new file mode 100644 index 00000000..02ac25fc --- /dev/null +++ b/toolbox/+otp/+cr3bp/radarMeasurementPlanar.m @@ -0,0 +1,14 @@ +function rmeas = radarMeasurementPlanar(~, y, ~, ~, radary) + +x = y(1:2, :); +v = y(3:4, :); + +relpos = x - radary; + +r = sqrt(relpos(1, :).^2 + relpos(2, :).^2); +drdt = sum(relpos.*v, 1)/r; +theta = atan2(relpos(2, :), relpos(1, :)); + +rmeas = [r; drdt; theta]; + +end From 35a3a34c6c50af43621771fd749db983ff30a10e Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 8 Oct 2025 13:05:41 -1000 Subject: [PATCH 29/33] renamed some files --- .../+presets/+utils/halocorrector.m | 2 +- .../+presets/Arenstorf.m | 9 ++++-- .../+presets/Canonical.m | 7 +++-- .../+presets/CircularEarthOrbit.m | 9 ++++-- .../+presets/HaloOrbit.m | 8 +++-- .../+presets/private/haloorbits.mat | Bin .../CR3BPParameters.m | 0 .../CR3BPProblem.m | 28 +++++++++--------- .../{+cr3bp => +circularrestricted3body}/f.m | 0 .../fPlanar.m | 0 .../jacobiConstant.m | 0 .../jacobiConstantJacobian.m | 0 .../jacobiConstantJacobianPlanar.m | 0 .../jacobiConstantPlanar.m | 0 .../jacobian.m | 0 .../jacobianPlanar.m | 0 .../jacobianVectorProduct.m | 0 .../jacobianVectorProductPlanar.m | 0 .../radarMeasurement.m | 0 .../radarMeasurementPlanar.m | 0 20 files changed, 36 insertions(+), 27 deletions(-) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/+presets/+utils/halocorrector.m (90%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/+presets/Arenstorf.m (69%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/+presets/Canonical.m (60%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/+presets/CircularEarthOrbit.m (83%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/+presets/HaloOrbit.m (97%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/+presets/private/haloorbits.mat (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/CR3BPParameters.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/CR3BPProblem.m (80%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/f.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/fPlanar.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/jacobiConstant.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/jacobiConstantJacobian.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/jacobiConstantJacobianPlanar.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/jacobiConstantPlanar.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/jacobian.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/jacobianPlanar.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/jacobianVectorProduct.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/jacobianVectorProductPlanar.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/radarMeasurement.m (100%) rename toolbox/+otp/{+cr3bp => +circularrestricted3body}/radarMeasurementPlanar.m (100%) diff --git a/toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m b/toolbox/+otp/+circularrestricted3body/+presets/+utils/halocorrector.m similarity index 90% rename from toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m rename to toolbox/+otp/+circularrestricted3body/+presets/+utils/halocorrector.m index eed0a26f..8594f77c 100644 --- a/toolbox/+otp/+cr3bp/+presets/+utils/halocorrector.m +++ b/toolbox/+otp/+circularrestricted3body/+presets/+utils/halocorrector.m @@ -8,7 +8,7 @@ otable = zeros(20, 4); for ind = 1:nIndices - problem = otp.cr3bp.presets.HaloOrbit('OrbitType', orbittype, 'Index', ind); + problem = otp.circularrestricted3body.presets.HaloOrbit('OrbitType', orbittype, 'Index', ind); y0 = problem.Y0; y0 = [y0(1); y0(3); y0(5)]; diff --git a/toolbox/+otp/+cr3bp/+presets/Arenstorf.m b/toolbox/+otp/+circularrestricted3body/+presets/Arenstorf.m similarity index 69% rename from toolbox/+otp/+cr3bp/+presets/Arenstorf.m rename to toolbox/+otp/+circularrestricted3body/+presets/Arenstorf.m index 4f12c362..401ff12b 100644 --- a/toolbox/+otp/+cr3bp/+presets/Arenstorf.m +++ b/toolbox/+otp/+circularrestricted3body/+presets/Arenstorf.m @@ -1,4 +1,4 @@ -classdef Arenstorf < otp.cr3bp.CR3BPProblem +classdef Arenstorf < otp.circularrestricted3body.CR3BPProblem % One period of a satellite moving in an Earth-Moon system on a planar % orbit. See pages 129--130 of :cite:p:`HNW93` for more details. @@ -6,7 +6,10 @@ function obj = Arenstorf(varargin) % Create the Arenstorf CR3BP problem object. - params = otp.cr3bp.CR3BPParameters('Mu', 0.012277471, 'SoftFactor', 0, varargin{:}); + params = otp.circularrestricted3body.CR3BPParameters(... + 'Mu', 0.012277471, ... + 'SoftFactor', 0, ... + varargin{:}); cls = class(params.Mu); % Decimals converted to rational to support multiple data types @@ -17,7 +20,7 @@ tspan = [cast(0, cls); ... cast(4541277234950502, cls) / cast(266113073862361, cls)]; - obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); + obj = obj@otp.circularrestricted3body.CR3BPProblem(tspan, y0, params); end end diff --git a/toolbox/+otp/+cr3bp/+presets/Canonical.m b/toolbox/+otp/+circularrestricted3body/+presets/Canonical.m similarity index 60% rename from toolbox/+otp/+cr3bp/+presets/Canonical.m rename to toolbox/+otp/+circularrestricted3body/+presets/Canonical.m index a84d72b3..ebcbec28 100644 --- a/toolbox/+otp/+cr3bp/+presets/Canonical.m +++ b/toolbox/+otp/+circularrestricted3body/+presets/Canonical.m @@ -1,4 +1,4 @@ -classdef Canonical < otp.cr3bp.CR3BPProblem +classdef Canonical < otp.circularrestricted3body.CR3BPProblem % A trivial preset with a stable oscillating orbit around a lagrange point. methods @@ -9,10 +9,11 @@ y0 = [0; 0; 0; 0; 0; 1]; tspan = [0, 10]; - params = otp.cr3bp.CR3BPParameters('Mu', mu, ... + params = otp.circularrestricted3body.CR3BPParameters(... + 'Mu', mu, ... 'SoftFactor', 1e-3, ... varargin{:}); - obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); + obj = obj@otp.circularrestricted3body.CR3BPProblem(tspan, y0, params); end end end diff --git a/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m b/toolbox/+otp/+circularrestricted3body/+presets/CircularEarthOrbit.m similarity index 83% rename from toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m rename to toolbox/+otp/+circularrestricted3body/+presets/CircularEarthOrbit.m index 8a98e15a..dd6c1430 100644 --- a/toolbox/+otp/+cr3bp/+presets/CircularEarthOrbit.m +++ b/toolbox/+otp/+circularrestricted3body/+presets/CircularEarthOrbit.m @@ -1,4 +1,4 @@ -classdef CircularEarthOrbit < otp.cr3bp.CR3BPProblem +classdef CircularEarthOrbit < otp.circularrestricted3body.CR3BPProblem % This preset builds a circular earth orbit in CR3BP based on % the equations derived in :cite:p:She20. @@ -47,8 +47,11 @@ y0 = [x0; 0; 0; 0; y0; 0]; tspan = [0, period]; - params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 1e-3, varargin{:}); - obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); + params = otp.circularrestricted3body.CR3BPParameters(... + 'Mu', mu, ... + 'SoftFactor', 1e-3, ... + varargin{:}); + obj = obj@otp.circularrestricted3body.CR3BPProblem(tspan, y0, params); end end end diff --git a/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m b/toolbox/+otp/+circularrestricted3body/+presets/HaloOrbit.m similarity index 97% rename from toolbox/+otp/+cr3bp/+presets/HaloOrbit.m rename to toolbox/+otp/+circularrestricted3body/+presets/HaloOrbit.m index db56f585..4ef66d8c 100644 --- a/toolbox/+otp/+cr3bp/+presets/HaloOrbit.m +++ b/toolbox/+otp/+circularrestricted3body/+presets/HaloOrbit.m @@ -1,4 +1,4 @@ -classdef HaloOrbit < otp.cr3bp.CR3BPProblem +classdef HaloOrbit < otp.circularrestricted3body.CR3BPProblem % This preset builds a halo orbit preset for the CR3BP based on % a table of reference initial conditions and periods. There are four % possible different orbits given by this preset, $L_2$, $P2HO_1$, @@ -58,8 +58,10 @@ y0 = [ic(1); 0; ic(2); 0; ic(3); 0]; tspan = [0, orbitalperiod]; - params = otp.cr3bp.CR3BPParameters('Mu', mu, 'SoftFactor', 0); - obj = obj@otp.cr3bp.CR3BPProblem(tspan, y0, params); + params = otp.circularrestricted3body.CR3BPParameters(... + 'Mu', mu, ... + 'SoftFactor', 0); + obj = obj@otp.circularrestricted3body.CR3BPProblem(tspan, y0, params); end end end diff --git a/toolbox/+otp/+cr3bp/+presets/private/haloorbits.mat b/toolbox/+otp/+circularrestricted3body/+presets/private/haloorbits.mat similarity index 100% rename from toolbox/+otp/+cr3bp/+presets/private/haloorbits.mat rename to toolbox/+otp/+circularrestricted3body/+presets/private/haloorbits.mat diff --git a/toolbox/+otp/+cr3bp/CR3BPParameters.m b/toolbox/+otp/+circularrestricted3body/CR3BPParameters.m similarity index 100% rename from toolbox/+otp/+cr3bp/CR3BPParameters.m rename to toolbox/+otp/+circularrestricted3body/CR3BPParameters.m diff --git a/toolbox/+otp/+cr3bp/CR3BPProblem.m b/toolbox/+otp/+circularrestricted3body/CR3BPProblem.m similarity index 80% rename from toolbox/+otp/+cr3bp/CR3BPProblem.m rename to toolbox/+otp/+circularrestricted3body/CR3BPProblem.m index 048b9e8d..0f1910e4 100644 --- a/toolbox/+otp/+cr3bp/CR3BPProblem.m +++ b/toolbox/+otp/+circularrestricted3body/CR3BPProblem.m @@ -58,7 +58,7 @@ % % Example % ------- - % >>> problem = otp.cr3bp.presets.HaloOrbit('OrbitType', 'L2', 'Index', 10); + % >>> problem = otp.circularrestricted3body.presets.HaloOrbit('OrbitType', 'L2', 'Index', 10); % >>> sol = model.solve(); % >>> problem.plotPhaseSpace(sol); % @@ -76,7 +76,7 @@ % The start and final time. % y0 : numeric(n, 1) % The initial conditions. Either n = 4 or n = 6. - % parameters : otp.cr3bp.CR3BPParameters + % parameters : otp.circularrestricted3body.CR3BPParameters % The parameters. obj@otp.Problem('Circular Restricted 3 Body Problem', [], timeSpan, y0, parameters); @@ -100,29 +100,29 @@ function onSettingsChanged(obj) end if spatialdim == 3 - obj.JacobiConstant = @(y) otp.cr3bp.jacobiConstant(y, mu, soft); - obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiConstantJacobian(y, mu, soft); + obj.JacobiConstant = @(y) otp.circularrestricted3body.jacobiConstant(y, mu, soft); + obj.JacobiConstantJacobian = @(y) otp.circularrestricted3body.jacobiConstantJacobian(y, mu, soft); - obj.RadarMeasurement = @(y, radary) otp.cr3bp.radarMeasurement(t, y, mu, soft, radary); + obj.RadarMeasurement = @(y, radary) otp.circularrestricted3body.radarMeasurement(t, y, mu, soft, radary); - obj.RHS = otp.RHS(@(t, y) otp.cr3bp.f(t, y, mu, soft), ... + obj.RHS = otp.RHS(@(t, y) otp.circularrestricted3body.f(t, y, mu, soft), ... 'Jacobian', ... - @(t, y) otp.cr3bp.jacobian(t, y, mu, soft), ... + @(t, y) otp.circularrestricted3body.jacobian(t, y, mu, soft), ... 'JacobianVectorProduct', ... - @(t, y, v) otp.cr3bp.jacobianVectorProduct(t, y, v, mu, soft), ... + @(t, y, v) otp.circularrestricted3body.jacobianVectorProduct(t, y, v, mu, soft), ... 'Vectorized', ... 'on'); else - obj.JacobiConstant = @(y) otp.cr3bp.jacobiConstantPlanar(y, mu, soft); - obj.JacobiConstantJacobian = @(y) otp.cr3bp.jacobiConstantJacobianPlanar(y, mu, soft); + obj.JacobiConstant = @(y) otp.circularrestricted3body.jacobiConstantPlanar(y, mu, soft); + obj.JacobiConstantJacobian = @(y) otp.circularrestricted3body.jacobiConstantJacobianPlanar(y, mu, soft); - obj.RadarMeasurement = @(y, radary) otp.cr3bp.radarMeasurementPlanar(t, y, mu, soft, radary); + obj.RadarMeasurement = @(y, radary) otp.circularrestricted3body.radarMeasurementPlanar(t, y, mu, soft, radary); - obj.RHS = otp.RHS(@(t, y) otp.cr3bp.fPlanar(t, y, mu, soft), ... + obj.RHS = otp.RHS(@(t, y) otp.circularrestricted3body.fPlanar(t, y, mu, soft), ... 'Jacobian', ... - @(t, y) otp.cr3bp.jacobianPlanar(t, y, mu, soft), ... + @(t, y) otp.circularrestricted3body.jacobianPlanar(t, y, mu, soft), ... 'JacobianVectorProduct', ... - @(t, y, v) otp.cr3bp.jacobianVectorProductPlanar(t, y, v, mu, soft), ... + @(t, y, v) otp.circularrestricted3body.jacobianVectorProductPlanar(t, y, v, mu, soft), ... 'Vectorized', ... 'on'); end diff --git a/toolbox/+otp/+cr3bp/f.m b/toolbox/+otp/+circularrestricted3body/f.m similarity index 100% rename from toolbox/+otp/+cr3bp/f.m rename to toolbox/+otp/+circularrestricted3body/f.m diff --git a/toolbox/+otp/+cr3bp/fPlanar.m b/toolbox/+otp/+circularrestricted3body/fPlanar.m similarity index 100% rename from toolbox/+otp/+cr3bp/fPlanar.m rename to toolbox/+otp/+circularrestricted3body/fPlanar.m diff --git a/toolbox/+otp/+cr3bp/jacobiConstant.m b/toolbox/+otp/+circularrestricted3body/jacobiConstant.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobiConstant.m rename to toolbox/+otp/+circularrestricted3body/jacobiConstant.m diff --git a/toolbox/+otp/+cr3bp/jacobiConstantJacobian.m b/toolbox/+otp/+circularrestricted3body/jacobiConstantJacobian.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobiConstantJacobian.m rename to toolbox/+otp/+circularrestricted3body/jacobiConstantJacobian.m diff --git a/toolbox/+otp/+cr3bp/jacobiConstantJacobianPlanar.m b/toolbox/+otp/+circularrestricted3body/jacobiConstantJacobianPlanar.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobiConstantJacobianPlanar.m rename to toolbox/+otp/+circularrestricted3body/jacobiConstantJacobianPlanar.m diff --git a/toolbox/+otp/+cr3bp/jacobiConstantPlanar.m b/toolbox/+otp/+circularrestricted3body/jacobiConstantPlanar.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobiConstantPlanar.m rename to toolbox/+otp/+circularrestricted3body/jacobiConstantPlanar.m diff --git a/toolbox/+otp/+cr3bp/jacobian.m b/toolbox/+otp/+circularrestricted3body/jacobian.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobian.m rename to toolbox/+otp/+circularrestricted3body/jacobian.m diff --git a/toolbox/+otp/+cr3bp/jacobianPlanar.m b/toolbox/+otp/+circularrestricted3body/jacobianPlanar.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobianPlanar.m rename to toolbox/+otp/+circularrestricted3body/jacobianPlanar.m diff --git a/toolbox/+otp/+cr3bp/jacobianVectorProduct.m b/toolbox/+otp/+circularrestricted3body/jacobianVectorProduct.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobianVectorProduct.m rename to toolbox/+otp/+circularrestricted3body/jacobianVectorProduct.m diff --git a/toolbox/+otp/+cr3bp/jacobianVectorProductPlanar.m b/toolbox/+otp/+circularrestricted3body/jacobianVectorProductPlanar.m similarity index 100% rename from toolbox/+otp/+cr3bp/jacobianVectorProductPlanar.m rename to toolbox/+otp/+circularrestricted3body/jacobianVectorProductPlanar.m diff --git a/toolbox/+otp/+cr3bp/radarMeasurement.m b/toolbox/+otp/+circularrestricted3body/radarMeasurement.m similarity index 100% rename from toolbox/+otp/+cr3bp/radarMeasurement.m rename to toolbox/+otp/+circularrestricted3body/radarMeasurement.m diff --git a/toolbox/+otp/+cr3bp/radarMeasurementPlanar.m b/toolbox/+otp/+circularrestricted3body/radarMeasurementPlanar.m similarity index 100% rename from toolbox/+otp/+cr3bp/radarMeasurementPlanar.m rename to toolbox/+otp/+circularrestricted3body/radarMeasurementPlanar.m From 63459344eb77bf8cbfb0c9b603bd41118a092a6f Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 8 Oct 2025 13:06:11 -1000 Subject: [PATCH 30/33] delete this file finally. --- .../+otp/+cr3bp/+presets/+utils/buildOrbit.m | 54 ------------------- 1 file changed, 54 deletions(-) delete mode 100644 toolbox/+otp/+cr3bp/+presets/+utils/buildOrbit.m diff --git a/toolbox/+otp/+cr3bp/+presets/+utils/buildOrbit.m b/toolbox/+otp/+cr3bp/+presets/+utils/buildOrbit.m deleted file mode 100644 index 05bf9521..00000000 --- a/toolbox/+otp/+cr3bp/+presets/+utils/buildOrbit.m +++ /dev/null @@ -1,54 +0,0 @@ -problem = otp.cr3bp.presets.LowEarthOrbit; - -y0 = problem.Y0; -tend = problem.TimeSpan(2); -f = problem.RHS.F; -mu = problem.Parameters.Mu; - -v0 = -7; - -J = @(v0) cost(v0, y0, f, tend, mu); -opts = optimoptions("fmincon", "Algorithm","interior-point", ... - "FunctionTolerance", 1e-14, "ConstraintTolerance", 1e-14, 'Display','off'); -[v0, fev] = fmincon(J, v0, [], [], [], [], [], [], [], opts); - -v0 - -% find period - -J = @(tend) cost2(v0, y0, f, tend, mu); -opts = optimoptions("fmincon", "Algorithm","interior-point", ... - "FunctionTolerance", 1e-14, "ConstraintTolerance", 1e-14, 'Display','off'); -[tend, fev] = fmincon(J, tend, [], [], [], [], 0.13, 0.15, [], opts); - -tend - - -function c = cost(v0, y0, f, tend, mu) - -y0(5) = v0; - -sol = ode45(f, [0, tend], y0, ... - odeset('AbsTol', 1e-14, 'RelTol', 100*eps)); - -yend = sol.y(:, :); - -earth = [-mu; 0; 0]; - -dy0 = norm(y0(1:3) - earth); -dyend = vecnorm(yend(1:3, :) - earth); - -c = sum((dy0 - dyend).^2); - -end - -function c = cost2(v0, y0, f, tend, ~) - -y0(5) = v0; - -sol = ode45(f, [0, tend], y0, ... - odeset('AbsTol', 1e-14, 'RelTol', 100*eps)); - -c = sum((y0 - sol.y(:, end)).^2); - -end From afe2d10389d3a339625a9407316c4d4b84a07378 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 8 Oct 2025 14:15:47 -1000 Subject: [PATCH 31/33] modified documentation --- docs/references.bib | 7 +++ .../+presets/Arenstorf.m | 5 +- .../+circularrestricted3body/CR3BPProblem.m | 53 +++++++++++++++---- 3 files changed, 54 insertions(+), 11 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index b60dac78..e7ceda0d 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -224,4 +224,11 @@ @book{She20 publisher={The University of North Dakota} } +@book{Are63, + title={Periodic solutions of the restricted three body problem representing analytic continuations of Keplerian elliptic motions}, + author={Arenstorf, Richard F}, + year={1963}, + publisher={National Aeronautics and Space Administration} +} + diff --git a/toolbox/+otp/+circularrestricted3body/+presets/Arenstorf.m b/toolbox/+otp/+circularrestricted3body/+presets/Arenstorf.m index 401ff12b..c4e12914 100644 --- a/toolbox/+otp/+circularrestricted3body/+presets/Arenstorf.m +++ b/toolbox/+otp/+circularrestricted3body/+presets/Arenstorf.m @@ -1,6 +1,9 @@ classdef Arenstorf < otp.circularrestricted3body.CR3BPProblem % One period of a satellite moving in an Earth-Moon system on a planar - % orbit. See pages 129--130 of :cite:p:`HNW93` for more details. + % orbit. The original orbit was derived in :cite:p:`Are63` for the + % planar circular restricted three body problem, and the setting from + % which this preset is derived is taken form pages 129--130 of + % :cite:p:`HNW93`. methods function obj = Arenstorf(varargin) diff --git a/toolbox/+otp/+circularrestricted3body/CR3BPProblem.m b/toolbox/+otp/+circularrestricted3body/CR3BPProblem.m index 0f1910e4..620e5597 100644 --- a/toolbox/+otp/+circularrestricted3body/CR3BPProblem.m +++ b/toolbox/+otp/+circularrestricted3body/CR3BPProblem.m @@ -8,8 +8,8 @@ % total mass of the system is represented by the non-dimensional % constant $\mu$. The reference frame of the system is fixed to the % rotationg frame of the two objects, meaning that the objects have - % fixed constant positions of $(\mu,0,0)^T$ for the first object, and - % $(1 - \mu,0,0)^T$ for the second object. The evolution of the third + % fixed constant positions of $[-\mu, 0, 0]^T$ for the first object, and + % $[1 - \mu, 0, 0]^T$ for the second object. The evolution of the third % object of negligent mass is given by the following second-order % non-dimensionalized differential equation: % @@ -19,18 +19,32 @@ % z'' &= \frac{\partial U}{\partial z}, % $$ % - % where + % where the energy and distances are defined as, % % $$ - % U &= \frac{1}{2} (x^2 + y^2) + \frac{1 - \mu}{d} + \frac{mu}{r},\\ + % U &= \frac{1}{2} (x^2 + y^2) + \frac{1 - \mu}{d} + \frac{\mu}{r},\\ % d &= \sqrt{(x + \mu)^2 + y^2 + z^2},\\ % r &= \sqrt{(x - 1 + \mu)^2 + y^2 + z^2}, % $$ % - % and where the system is converted to a differential equation in six - % variables in the standard fashion. The distances $d$ and $r$ can - % cause numerical instability as they approach zero, thus a softening - % factor of $s^2$ is typically added under both of the square-roots. + % and where the system is converted to a first order differential + % equation in six variables in the standard fashion as, + % + % $$ + % x' = v_x,\\ + % y' = v_y,\\ + % z' = v_z,\\ + % v_x' = x'',\\ + % v_y' = y'',\\ + % v_z' = z'', + % $$ + % + % where the new variables $v_x$, $v_y$, and $v_z$ represent the + % velocity vector. + % + % The distances $d$ and $r$ can cause numerical instability as they + % approach zero, thus a softening factor of $s^2$ is typically added + % under both of the square-roots of the distances $d$ and $r$. % % When the object under consideration is on an orbit that is co-planar % to the orbit of the two other objects, then the system of equations @@ -40,11 +54,30 @@ % constant of the system, % % $$ - % J = 2U - x'^2 - y'^2 - z'^2, + % J = 2U - v_x^2 - v_y^2 - v_z^2, % $$ % % is preserved throughout the evolution of the equations, though this - % is typically not true numerically. + % is typically not true numerically. The Jacobi constant is provided as + % a function. + % + % For state estimation purposes a radar measurement is also provided, + % + % $$ + % h(x, y, z, v_x, v_y, v_z) = \begin{bmatrix} + % \sqrt{(x - s_x)^2 + (y - s_y)^2 + (z - s_z)^2}\\ + % \frac{(x - s_x) \cdot v_x + (y - s_y) \cdot v_y + (z - s_z) \cdot v_z} + % {\sqrt{(x - s_x)^2 + (y - s_y)^2 + (z - s_z)^2}}\\ + % \tan^{-1} \frac{y}{x + \mu}\\ + % \tan^{-1} \frac{z}{\sqrt{(x - s_x)^2 + (y - s_y)^2 + (z - s_z)^2}} + % \end{bmatrix} + % $$ + % + % where $s_x$, $s_y$, and $s_z$ are the locations of the radar sensor, + % the first term is the range to the object from the sensor, the second + % term is the rate of change of the range (range-rate), and the third + % and fourth terms are the two angles the radar must be pointing. In + % the planar case, the fourth term is absent. % % Notes % ----- From 0e4fd7e82ef6874e8ad23b908c12058b6964231f Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 8 Oct 2025 15:02:25 -1000 Subject: [PATCH 32/33] minor fixes for docs --- .../+presets/CircularEarthOrbit.m | 2 +- .../+circularrestricted3body/+presets/HaloOrbit.m | 15 ++++++++------- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/toolbox/+otp/+circularrestricted3body/+presets/CircularEarthOrbit.m b/toolbox/+otp/+circularrestricted3body/+presets/CircularEarthOrbit.m index dd6c1430..6cf75420 100644 --- a/toolbox/+otp/+circularrestricted3body/+presets/CircularEarthOrbit.m +++ b/toolbox/+otp/+circularrestricted3body/+presets/CircularEarthOrbit.m @@ -1,6 +1,6 @@ classdef CircularEarthOrbit < otp.circularrestricted3body.CR3BPProblem % This preset builds a circular earth orbit in CR3BP based on - % the equations derived in :cite:p:She20. + % the equations derived in :cite:p:`She20`. methods function obj = CircularEarthOrbit(varargin) diff --git a/toolbox/+otp/+circularrestricted3body/+presets/HaloOrbit.m b/toolbox/+otp/+circularrestricted3body/+presets/HaloOrbit.m index 4ef66d8c..e0d66a11 100644 --- a/toolbox/+otp/+circularrestricted3body/+presets/HaloOrbit.m +++ b/toolbox/+otp/+circularrestricted3body/+presets/HaloOrbit.m @@ -1,11 +1,12 @@ classdef HaloOrbit < otp.circularrestricted3body.CR3BPProblem - % This preset builds a halo orbit preset for the CR3BP based on - % a table of reference initial conditions and periods. There are four - % possible different orbits given by this preset, $L_2$, $P2HO_1$, - % $P2HO_2$, $P4HO_1$, and $P4HO_2$. The tables contain $20$ entries, and - % thus there are $100$ different possible orbits given by this preset. - % The tables of $L_2$ halo orbits are given as Table A.1--A.5 on pages - % 218--222 of :cite:p:`Spr21`. + % This preset builds a halo orbit around the Moon in the Earth-Moon + % system preset for the CR3BP based on a table of reference initial + % conditions and periods. There are four possible different orbits + % given by this preset, $L_2$, $P2HO_1$, $P2HO_2$, $P4HO_1$, and + % $P4HO_2$. The tables contain $20$ entries, and thus there are $100$ + % different possible orbits given by this preset. The tables of $L_2$ + % halo orbits are given as Table A.1--A.5 on pages 218--222 of + % :cite:p:`Spr21`. methods function obj = HaloOrbit(varargin) From b5e4d94335ac5cbf1115d36ed7e3d1c31180e35d Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Wed, 8 Oct 2025 15:02:44 -1000 Subject: [PATCH 33/33] Delete old Arenstorf problem --- toolbox/+otp/+arenstorf/+presets/Canonical.m | 34 ---------- toolbox/+otp/+arenstorf/ArenstorfParameters.m | 22 ------ toolbox/+otp/+arenstorf/ArenstorfProblem.m | 67 ------------------- toolbox/+otp/+arenstorf/f.m | 18 ----- toolbox/+otp/+arenstorf/jacobian.m | 27 -------- 5 files changed, 168 deletions(-) delete mode 100644 toolbox/+otp/+arenstorf/+presets/Canonical.m delete mode 100644 toolbox/+otp/+arenstorf/ArenstorfParameters.m delete mode 100644 toolbox/+otp/+arenstorf/ArenstorfProblem.m delete mode 100644 toolbox/+otp/+arenstorf/f.m delete mode 100644 toolbox/+otp/+arenstorf/jacobian.m diff --git a/toolbox/+otp/+arenstorf/+presets/Canonical.m b/toolbox/+otp/+arenstorf/+presets/Canonical.m deleted file mode 100644 index 6ce4917c..00000000 --- a/toolbox/+otp/+arenstorf/+presets/Canonical.m +++ /dev/null @@ -1,34 +0,0 @@ -classdef Canonical < otp.arenstorf.ArenstorfProblem - %CANONICAL One period of a satellite moving in an Earth-Moon system - % Source: - % - % Hairer, Ernst, et al. Solving Ordinary Differential Equations I: - % Nonstiff Problems. 2nd ed., vol. 8, Springer-Verlag, 1993, pp. 129-130. - % - % See also otp.arenstorf.ArenstorfProblem - - methods - function obj = Canonical(varargin) - %CANONICAL Construct a canonical Arenstorf problem - % OBJ = CANONICAL(MU) Uses a Moon of mass MU. The time span and - % initial conditions use the same data type as MU. Single - % precision, vpa, and other types are supported. - % - % OBJ = CANONICAL() Uses a MU corresponding to the Moon - - params = otp.arenstorf.ArenstorfParameters('Mu', 0.012277471, varargin{:}); - cls = class(params.Mu); - - % Decimals converted to rational to support multiple data types - y0 = [cast(497, cls) / cast(500, cls); ... - cast(0, cls); ... - cast(0, cls); ... - cast(-823970832321143, cls) / cast(411659154384760, cls)]; - tspan = [cast(0, cls); ... - cast(4541277234950502, cls) / cast(266113073862361, cls)]; - - obj = obj@otp.arenstorf.ArenstorfProblem(tspan, y0, params); - end - - end -end diff --git a/toolbox/+otp/+arenstorf/ArenstorfParameters.m b/toolbox/+otp/+arenstorf/ArenstorfParameters.m deleted file mode 100644 index 449a6845..00000000 --- a/toolbox/+otp/+arenstorf/ArenstorfParameters.m +++ /dev/null @@ -1,22 +0,0 @@ -classdef ArenstorfParameters < otp.Parameters - % Parameters for the Arenstorf problem. - - properties - % The mass of one body, while the other body has mass $\mu' = 1 - \mu$. - Mu %MATLAB ONLY: (1, 1) {otp.utils.validation.mustBeNumerical} - end - - methods - function obj = ArenstorfParameters(varargin) - % Create an Arenstorf parameters object. - % - % Parameters - % ---------- - % varargin - % A variable number of name-value pairs. A name can be any property of this class, and the subsequent - % value initializes that property. - - obj = obj@otp.Parameters(varargin{:}); - end - end -end diff --git a/toolbox/+otp/+arenstorf/ArenstorfProblem.m b/toolbox/+otp/+arenstorf/ArenstorfProblem.m deleted file mode 100644 index f53dd718..00000000 --- a/toolbox/+otp/+arenstorf/ArenstorfProblem.m +++ /dev/null @@ -1,67 +0,0 @@ -classdef ArenstorfProblem < otp.Problem - %ARENSTORFPROBLEM A simplified 3-body problem - % This problem models the trajectory of a satellite of negligible mass - % under the gravitational pull of two orbiting bodies. These can be - % considered the Moon and Earth with masses MU and MU'=1-MU, respectively. - % The reference frame is chosen such that the Earth is at position - % (-MU, 0) on the plane, and the moon is at position (MU', 0). The - % satellite's motion is described by the second-order ODEs - % - % y1'' = y1 + 2 y2' - MU' (y1 + MU) / D1 - MU (y1 - MU') / D2, - % y2'' = y2 - 2 y1' - MU' y2 / D1 - MU y2 / D2, - % - % where - % - % D1 = ((y1 + MU)^2 + y2^2)^(3/2), D2 = ((y1 - MU')^2 + y2^2)^(3/2). - % - % This is converted into a system of first order ODEs in which the state - % contains two position variables and two velocity variables: - % - % y = [y1; y2; y1'; y2']. - % - % Sources: - % - % Arenstorf, Richard F. "Periodic Solutions of the Restricted Three Body - % Problem Representing Analytic Continuations of Keplerian Elliptic - % Motions." American Journal of Mathematics, vol. 85, no. 1, 1963, p. 27. - % - % Hairer, Ernst, et al. Solving Ordinary Differential Equations I: - % Nonstiff Problems. 2nd ed., vol. 8, Springer-Verlag, 1993, pp. 129-130. - % - % See also otp.arenstorf.ArenstorfParameters - - methods - %ARENSTORFPROBLEM Construct an Arenstorf problem - function obj = ArenstorfProblem(timeSpan, y0, parameters) - obj@otp.Problem('Arenstorf Orbit', 4, timeSpan, y0, parameters); - end - end - - methods (Access=protected) - function onSettingsChanged(obj) - mu = obj.Parameters.Mu; - - obj.RHS = otp.RHS(@(t, y) otp.arenstorf.f(t, y, mu), ... - 'Jacobian', @(t, y) otp.arenstorf.jacobian(t, y, mu)); - end - - function label = internalIndex2label(~, index) - if index <= 2 - label = sprintf('y_{%d}', index); - else - label = sprintf('y''_{%d}', index - 2); - end - end - - function fig = internalPlotPhaseSpace(obj, t, y, varargin) - fig = internalPlotPhaseSpace@otp.Problem(obj, t, y, ... - 'Vars', 1:2, varargin{:}); - end - - function mov = internalMovie(obj, t, y, varargin) - mov = otp.utils.movie.PhaseSpaceMovie('title', obj.Name, ... - 'Vars', 1:2, varargin{:}); - mov.record(t, y); - end - end -end diff --git a/toolbox/+otp/+arenstorf/f.m b/toolbox/+otp/+arenstorf/f.m deleted file mode 100644 index 18a7d6bf..00000000 --- a/toolbox/+otp/+arenstorf/f.m +++ /dev/null @@ -1,18 +0,0 @@ -function dy = f(~, y, mu) -%F The right-hand side function for the Arenstorf problem -% -% See also otp.arenstorf.ArenstorfProblem - -muPrime = 1 - mu; - -p1 = mu + y(1); -p2 = y(1) - muPrime; - -d1 = (p1^2 + y(2)^2)^-1.5; -d2 = (p2^2 + y(2)^2)^-1.5; - -dy = [y(3:4); ... - y(1) + 2 * y(4) - muPrime * d1 * p1 - mu * d2 * p2; ... - y(2) - 2 * y(3) - (muPrime * d1 + mu * d2) * y(2)]; - -end diff --git a/toolbox/+otp/+arenstorf/jacobian.m b/toolbox/+otp/+arenstorf/jacobian.m deleted file mode 100644 index 29b8803a..00000000 --- a/toolbox/+otp/+arenstorf/jacobian.m +++ /dev/null @@ -1,27 +0,0 @@ -function j = jacobian(~, y, mu) -%JACOBIAN The Jacobian for the Arenstorf problem -% -% See also otp.arenstorf.f, otp.arenstorf.ArenstorfProblem - -muPrime = 1 - mu; - -p1 = mu + y(1); -p2 = y(1) - muPrime; - -distSq1 = p1^2 + y(2)^2; -distSq2 = p2^2 + y(2)^2; - -d1 = muPrime * distSq1^-1.5; -d2 = mu * distSq2^-1.5; - -l1 = 3 * muPrime * distSq1^-2.5; -l2 = 3 * mu * distSq2^-2.5; - -offDiag = y(2) * (l1 * p1 + l2 * p2); - -j = [0, 0, 1, 0; ... - 0, 0, 0, 1; ... - 1 + l1 * p1^2 - d1 + l2 * p2^2 - d2, offDiag, 0, 2; ... - offDiag, 1 + y(2)^2 * (l1 + l2) - d1 - d2, -2, 0]; - -end