added Arentsorf Orbit Problem, needs custom Phaseplot with masses as point plots

Author: elswit

src/+otp/+arenstorf/+presets/Canonical.m

@@ -0,0 +1,14 @@
+classdef Canonical < otp.arenstorf.ArenstorfProblem
+    methods
+        function obj = Canonical
+            params.m1 = 0.012277471;
+            params.m2 = 1-0.012277471;
+
+            y0 = [0.994; 0; 0; -2.00158510637908252240537862224];
+            tspan = [0,20];
+            
+            obj = obj@otp.arenstorf.ArenstorfProblem(tspan, y0, params);
+        end
+
+    end
+end

src/+otp/+arenstorf/ArenstorfProblem.m

@@ -0,0 +1,38 @@
+classdef ArenstorfProblem < otp.Problem
+  % The classic  3-body  problem with parameters from
+  % Reference:  Hairer, Ernst, Syvert Paul Nørsett, and Gerhard Wanner.
+  % Solving Ordinary Differential Equations I: Nonstiff Problems Springer-Verlag, 1987.
+  % CH II, p. 129.
+  methods
+    function obj = ArenstorfProblem(timeSpan, y0, parameters)
+      obj@otp.Problem('Arenstorf Orbit Problem', 4, timeSpan, y0, parameters);
+    end
+  end
+  
+  methods (Access=protected)
+    
+    function onSettingsChanged(obj)
+      m1 = obj.Parameters.m1;
+      m2 = obj.Parameters.m2;
+      
+      obj.Rhs = otp.Rhs( @(t, y) otp.arenstorf.f(t, y, m1, m2), ...
+        otp.Rhs.FieldNames.Jacobian, @(t, y) otp.arenstorf.jac(t, y, m1, m2));
+    end
+    
+    
+    function validateNewState(obj, newTimeSpan, newY0, newParameters)
+      validateNewState@otp.Problem(obj, newTimeSpan, newY0, newParameters);
+      
+      otp.utils.StructParser(newParameters) ...
+        .checkField('m1', 'scalar', 'real', 'finite', 'positive') ...
+        .checkField('m2', 'scalar', 'real', 'finite', 'positive');
+      
+    end
+    
+    
+%     function fig = internalPlotPhaseSpace(obj, sol, varargin)
+%       fig = internalPlotPhaseSpace@otp.Problem(obj, sol, 1:2, varargin{:});
+%     end
+    
+  end
+end

src/+otp/+arenstorf/f.m

@@ -0,0 +1,12 @@
+function dy = f(~, y, m1 ,m2)
+
+    
+%     a=0.012277471; b=1.0-a; 
+    
+    D1= ((y(1)+m1)^2+y(2)^2)^(3/2);
+    D2= ((y(1)-m2)^2+y(2)^2)^(3/2);
+    
+    dy = [y(3);y(4);y(1)+2*y(4)-m2/D1*(y(1)+m1)-m1/D2*(y(1)-m2);y(2)-2*y(3)-m2/D1*y(2)-m1/D2*y(2)];
+
+
+end

src/+otp/+arenstorf/jac.m

@@ -0,0 +1,9 @@
+function dfy = jac(~, y, m1 ,m2)
+
+dfy =[
+                                                                                                                                                                                                 0,                                                                                                                                                   0,  1, 0;
+                                                                                                                                                                                                 0,                                                                                                                                                   0,  0, 1;
+ (3*m1*(2*m2 - 2*y(1))*(m2 - y(1)))/(2*((m2 - y(1))^2 + y(2)^2)^(5/2)) - m1/((m2 - y(1))^2 + y(2)^2)^(3/2) - m2/((m1 + y(1))^2 + y(2)^2)^(3/2) + (3*m2*(2*m1 + 2*y(1))*(m1 + y(1)))/(2*((m1 + y(1))^2 + y(2)^2)^(5/2)) + 1,                                                     (3*m2*y(2)*(m1 + y(1)))/((m1 + y(1))^2 + y(2)^2)^(5/2) - (3*m1*y(2)*(m2 - y(1)))/((m2 - y(1))^2 + y(2)^2)^(5/2),  0, 2;
+                                                                                           (3*m2*y(2)*(m1 + y(1)))/((m1 + y(1))^2 + y(2)^2)^(5/2) - (3*m1*y(2)*(2*m2 - 2*y(1)))/(2*((m2 - y(1))^2 + y(2)^2)^(5/2)), (3*m1*y(2)^2)/((m2 - y(1))^2 + y(2)^2)^(5/2) - m1/((m2 - y(1))^2 + y(2)^2)^(3/2) - m2/((m1 + y(1))^2 + y(2)^2)^(3/2) + (3*m2*y(2)^2)/((m1 + y(1))^2 + y(2)^2)^(5/2) + 1, -2, 0];
+
+end