Circular Restricted 3 Body Problem

.gitignore

@@ -10,3 +10,4 @@ docs/build/
 docs/problems/
 
 .DS_Store
+tests/octave-workspace

docs/Makefile

@@ -19,7 +19,7 @@ serve:
 	--ignore $(PROBLEMSDIR)
 
 lint:
-	pycodestyle .. && codespell ..
+	pycodestyle .. && codespell --skip="*.bib" ..
 
 clean:
 	rm -rf $(BUILDDIR) $(PROBLEMSDIR)

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}
+}
+
+

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
+
+

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 = [email protected](tspan, y0, params);
+        end
+
+    end
+end

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 = [email protected](tspan, y0, params);            
+        end
+    end
+end

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 = [email protected](tspan, y0, params);            
+        end
+    end
+end