added LEO preset

Author: AndreyAPopov

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

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

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