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/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) diff --git a/docs/references.bib b/docs/references.bib index 15a4c373..e7ceda0d 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -209,3 +209,26 @@ @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} +} + +@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} +} + +@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/+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/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 diff --git a/toolbox/+otp/+circularrestricted3body/+presets/+utils/halocorrector.m b/toolbox/+otp/+circularrestricted3body/+presets/+utils/halocorrector.m new file mode 100644 index 00000000..8594f77c --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/+presets/+utils/halocorrector.m @@ -0,0 +1,50 @@ +orbittypes = {'L2', 'P2HO1', 'P2HO2', 'P4HO1', 'P4HO2'}; +nIndices = 20; + +orbits = struct; + +for oti = 1:numel(orbittypes) + orbittype = orbittypes{oti}; + otable = zeros(20, 4); + for ind = 1:nIndices + + problem = otp.circularrestricted3body.presets.HaloOrbit('OrbitType', orbittype, '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.(orbittype) = 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/+circularrestricted3body/+presets/Arenstorf.m b/toolbox/+otp/+circularrestricted3body/+presets/Arenstorf.m new file mode 100644 index 00000000..c4e12914 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/+presets/Arenstorf.m @@ -0,0 +1,30 @@ +classdef Arenstorf < otp.circularrestricted3body.CR3BPProblem + % One period of a satellite moving in an Earth-Moon system on a planar + % 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) + % Create the Arenstorf CR3BP problem object. + + params = otp.circularrestricted3body.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.circularrestricted3body.CR3BPProblem(tspan, y0, params); + end + + end +end diff --git a/toolbox/+otp/+circularrestricted3body/+presets/Canonical.m b/toolbox/+otp/+circularrestricted3body/+presets/Canonical.m new file mode 100644 index 00000000..ebcbec28 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/+presets/Canonical.m @@ -0,0 +1,19 @@ +classdef Canonical < otp.circularrestricted3body.CR3BPProblem + % A trivial preset with a stable oscillating orbit around a lagrange point. + + methods + function obj = Canonical(varargin) + % Create the Canonical CR3BP problem object. + + mu = 0.5; + + y0 = [0; 0; 0; 0; 0; 1]; + tspan = [0, 10]; + 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/+circularrestricted3body/+presets/CircularEarthOrbit.m b/toolbox/+otp/+circularrestricted3body/+presets/CircularEarthOrbit.m new file mode 100644 index 00000000..6cf75420 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/+presets/CircularEarthOrbit.m @@ -0,0 +1,57 @@ +classdef CircularEarthOrbit < otp.circularrestricted3body.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. + % The default radius is 340 km, which is approximately where + % the ISS is/was located. + + p = inputParser(); + p.KeepUnmatched = true; + p.addParameter('OrbitalRadius', 340, @(x) isnumeric(x) && isscalar(x)); + p.parse(varargin{:}); + results = p.Results; + varargin = [fieldnames(p.Unmatched), struct2cell(p.Unmatched)].'; + + equatorialRadius = 6378; + earthMoonDist = 385000; + orbitalradius = results.OrbitalRadius; + + % set the relative distance from earth to the moon + 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.circularrestricted3body.CR3BPParameters(... + 'Mu', mu, ... + 'SoftFactor', 1e-3, ... + varargin{:}); + obj = obj@otp.circularrestricted3body.CR3BPProblem(tspan, y0, params); + end + end +end diff --git a/toolbox/+otp/+circularrestricted3body/+presets/HaloOrbit.m b/toolbox/+otp/+circularrestricted3body/+presets/HaloOrbit.m new file mode 100644 index 00000000..e0d66a11 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/+presets/HaloOrbit.m @@ -0,0 +1,203 @@ +classdef HaloOrbit < otp.circularrestricted3body.CR3BPProblem + % 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) + % 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. + % 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; + + if results.Corrected + T = gethalotablecorrected(results.OrbitType); + else + T = gethalotable(results.OrbitType); + end + + 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.circularrestricted3body.CR3BPParameters(... + 'Mu', mu, ... + 'SoftFactor', 0); + obj = obj@otp.circularrestricted3body.CR3BPProblem(tspan, y0, params); + end + end +end + + +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 + +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/+circularrestricted3body/+presets/private/haloorbits.mat b/toolbox/+otp/+circularrestricted3body/+presets/private/haloorbits.mat new file mode 100644 index 00000000..9a1e5bd6 Binary files /dev/null and b/toolbox/+otp/+circularrestricted3body/+presets/private/haloorbits.mat differ diff --git a/toolbox/+otp/+arenstorf/ArenstorfParameters.m b/toolbox/+otp/+circularrestricted3body/CR3BPParameters.m similarity index 53% rename from toolbox/+otp/+arenstorf/ArenstorfParameters.m rename to toolbox/+otp/+circularrestricted3body/CR3BPParameters.m index 449a6845..77232e81 100644 --- a/toolbox/+otp/+arenstorf/ArenstorfParameters.m +++ b/toolbox/+otp/+circularrestricted3body/CR3BPParameters.m @@ -1,14 +1,17 @@ -classdef ArenstorfParameters < otp.Parameters - % Parameters for the Arenstorf problem. - +classdef CR3BPParameters < otp.Parameters + % Parameters for the CR3BP problem. + properties - % The mass of one body, while the other body has mass $\mu' = 1 - \mu$. + % Relative mass of the system. Mu %MATLAB ONLY: (1, 1) {otp.utils.validation.mustBeNumerical} + + % A factor to avoid singularity when computing distances. + SoftFactor %MATLAB ONLY: (1, 1) {otp.utils.validation.mustBeNumerical} end methods - function obj = ArenstorfParameters(varargin) - % Create an Arenstorf parameters object. + function obj = CR3BPParameters(varargin) + % Create a CR3BP parameters object. % % Parameters % ---------- diff --git a/toolbox/+otp/+circularrestricted3body/CR3BPProblem.m b/toolbox/+otp/+circularrestricted3body/CR3BPProblem.m new file mode 100644 index 00000000..620e5597 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/CR3BPProblem.m @@ -0,0 +1,222 @@ +classdef CR3BPProblem < otp.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 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: + % + % $$ + % x'' &= 2y' + \frac{\partial U}{\partial x},\\ + % y'' &= -2x' + \frac{\partial U}{\partial y},\\ + % z'' &= \frac{\partial U}{\partial z}, + % $$ + % + % where the energy and distances are defined as, + % + % $$ + % 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 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 + % 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, + % + % $$ + % 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. 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 + % ----- + % +---------------------+-----------------------------------------------+ + % | Type | ODE | + % +---------------------+-----------------------------------------------+ + % | Number of Variables | 4 for planar or 6 for non-planar | + % +---------------------+-----------------------------------------------+ + % | Stiff | no | + % +---------------------+-----------------------------------------------+ + % + % Example + % ------- + % >>> problem = otp.circularrestricted3body.presets.HaloOrbit('OrbitType', 'L2', 'Index', 10); + % >>> sol = model.solve(); + % >>> 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(n, 1) + % The initial conditions. Either n = 4 or n = 6. + % parameters : otp.circularrestricted3body.CR3BPParameters + % The parameters. + + obj@otp.Problem('Circular Restricted 3 Body Problem', [], timeSpan, y0, parameters); + end + end + + properties (SetAccess = private) + JacobiConstant + JacobiConstantJacobian + RadarMeasurement + end + + methods (Access = protected) + function onSettingsChanged(obj) + mu = obj.Parameters.Mu; + soft = obj.Parameters.SoftFactor; + + spatialdim = numel(obj.Y0)/2; + if ~(spatialdim == 2 || spatialdim == 3) + 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.circularrestricted3body.jacobiConstant(y, mu, soft); + obj.JacobiConstantJacobian = @(y) otp.circularrestricted3body.jacobiConstantJacobian(y, mu, soft); + + obj.RadarMeasurement = @(y, radary) otp.circularrestricted3body.radarMeasurement(t, y, mu, soft, radary); + + obj.RHS = otp.RHS(@(t, y) otp.circularrestricted3body.f(t, y, mu, soft), ... + 'Jacobian', ... + @(t, y) otp.circularrestricted3body.jacobian(t, y, mu, soft), ... + 'JacobianVectorProduct', ... + @(t, y, v) otp.circularrestricted3body.jacobianVectorProduct(t, y, v, mu, soft), ... + 'Vectorized', ... + 'on'); + else + obj.JacobiConstant = @(y) otp.circularrestricted3body.jacobiConstantPlanar(y, mu, soft); + obj.JacobiConstantJacobian = @(y) otp.circularrestricted3body.jacobiConstantJacobianPlanar(y, mu, soft); + + obj.RadarMeasurement = @(y, radary) otp.circularrestricted3body.radarMeasurementPlanar(t, y, mu, soft, radary); + + obj.RHS = otp.RHS(@(t, y) otp.circularrestricted3body.fPlanar(t, y, mu, soft), ... + 'Jacobian', ... + @(t, y) otp.circularrestricted3body.jacobianPlanar(t, y, mu, soft), ... + 'JacobianVectorProduct', ... + @(t, y, v) otp.circularrestricted3body.jacobianVectorProductPlanar(t, y, v, mu, soft), ... + 'Vectorized', ... + 'on'); + end + end + + 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 + + function sol = internalSolve(obj, 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; + + 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 + +end diff --git a/toolbox/+otp/+circularrestricted3body/f.m b/toolbox/+otp/+circularrestricted3body/f.m new file mode 100644 index 00000000..9f1630ac --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/f.m @@ -0,0 +1,25 @@ +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); + +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); + +dx = v; +dv = [2*v(2, :) + dUdx; -2*v(1, :) + dUdy; dUdz]; + +dy = [dx; dv]; + +end diff --git a/toolbox/+otp/+circularrestricted3body/fPlanar.m b/toolbox/+otp/+circularrestricted3body/fPlanar.m new file mode 100644 index 00000000..c4efdfe3 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/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/+circularrestricted3body/jacobiConstant.m b/toolbox/+otp/+circularrestricted3body/jacobiConstant.m new file mode 100644 index 00000000..054a3245 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/jacobiConstant.m @@ -0,0 +1,13 @@ +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 diff --git a/toolbox/+otp/+circularrestricted3body/jacobiConstantJacobian.m b/toolbox/+otp/+circularrestricted3body/jacobiConstantJacobian.m new file mode 100644 index 00000000..8a114554 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/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 diff --git a/toolbox/+otp/+circularrestricted3body/jacobiConstantJacobianPlanar.m b/toolbox/+otp/+circularrestricted3body/jacobiConstantJacobianPlanar.m new file mode 100644 index 00000000..f5210966 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/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/+circularrestricted3body/jacobiConstantPlanar.m b/toolbox/+otp/+circularrestricted3body/jacobiConstantPlanar.m new file mode 100644 index 00000000..8068ae32 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/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/+circularrestricted3body/jacobian.m b/toolbox/+otp/+circularrestricted3body/jacobian.m new file mode 100644 index 00000000..f6733adf --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/jacobian.m @@ -0,0 +1,61 @@ +function J = jacobian(~, y, mu, soft) + +x = y(1:3, :); +x = reshape(x, 3, 1, []); +N = size(x, 3); + +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, 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, 1, :))./(d.^2); +dddxdz = -(dddx.*x(3, 1, :))./(d.^2); +dddydy = (1 - dddy.^2)./d; +dddydz = -(dddy.*x(3, 1, :))./(d.^2); +dddzdz = (1 - dddz.^2)./d; + +drdxdx = (1 - drdx.^2)./r; +drdxdy = -(drdx.*x(2, 1, :))./(r.^2); +drdxdz = -(drdx.*x(3, 1, :))./(r.^2); +drdydy = (1 - drdy.^2)./r; +drdydz = -(drdy.*x(3, 1, :))./(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, N), repmat(eye(3), 1, 1, N); ... + dxdv, repmat([0, 2, 0; -2, 0, 0; 0, 0, 0], 1, 1, N)]; + +end diff --git a/toolbox/+otp/+circularrestricted3body/jacobianPlanar.m b/toolbox/+otp/+circularrestricted3body/jacobianPlanar.m new file mode 100644 index 00000000..945363e2 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/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/+circularrestricted3body/jacobianVectorProduct.m b/toolbox/+otp/+circularrestricted3body/jacobianVectorProduct.m new file mode 100644 index 00000000..2678509b --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/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 diff --git a/toolbox/+otp/+circularrestricted3body/jacobianVectorProductPlanar.m b/toolbox/+otp/+circularrestricted3body/jacobianVectorProductPlanar.m new file mode 100644 index 00000000..ecbf95be --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/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 diff --git a/toolbox/+otp/+circularrestricted3body/radarMeasurement.m b/toolbox/+otp/+circularrestricted3body/radarMeasurement.m new file mode 100644 index 00000000..4cbbdd72 --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/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 diff --git a/toolbox/+otp/+circularrestricted3body/radarMeasurementPlanar.m b/toolbox/+otp/+circularrestricted3body/radarMeasurementPlanar.m new file mode 100644 index 00000000..02ac25fc --- /dev/null +++ b/toolbox/+otp/+circularrestricted3body/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 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