Merge pull request #139 from ComputationalScienceLaboratory/cr3bp

Circular Restricted 3 Body Problem

Author: AndreyAPopov

.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/+arenstorf/+presets/Canonical.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/+arenstorf/ArenstorfProblem.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

toolbox/+otp/+arenstorf/f.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