initial commit of CR3BP

Author: AndreyAPopov

toolbox/+otp/+cr3bp/+presets/Canonical.m

@@ -0,0 +1,3 @@
+classdef Canonical < otp.cr3bp.CR3BPProblem
+    
+end

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

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

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

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