Merge pull request #7 from ComputationalScienceLaboratory/ArenstorfOrbit

Added Arenstorf orbit problem

Author: elswit

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

@@ -0,0 +1,20 @@
+classdef Canonical < otp.arenstorf.ArenstorfProblem
+    methods
+        function obj = Canonical(varargin)
+                 
+            p = inputParser;
+            p.addParameter('m1', 0.012277471);
+            p.parse(varargin{:});
+            s = p.Results;
+            
+            params.m1 = s.m1;
+            params.m2 = 1-params.m1;
+                
+            y0 = [0.994; 0; 0; -2.00158510637908252240537862224];
+            tspan = [0,17.0652165601579625588917206249];
+            
+            obj = obj@otp.arenstorf.ArenstorfProblem(tspan, y0, params);
+        end
+
+    end
+end

src/+otp/+arenstorf/ArenstorfProblem.m

@@ -0,0 +1,37 @@
+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', 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 mov = internalMovie(obj, t, y, varargin)
+      mov = otp.utils.movie.PhasePlaneMovie(obj.Name, @obj.index2label, varargin{:});
+      mov.record(t, y(:,1:2));
+    end
+    
+  end
+end

src/+otp/+arenstorf/f.m

@@ -0,0 +1,10 @@
+function dy = f(~, y, m1 ,m2)
+
+    D1 = ((y(1) + m1)^2 + y(2)^2)^(3/2);
+    D2 = ((y(1) - m2)^2 + y(2)^2)^(3/2);
+    
+    dy = [y(3: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,14 @@
+function dfy = jac(~, y, m1 ,m2)
+
+D1 = m1/((m2 - y(1))^2 + y(2)^2)^(3/2);
+D2 = m2/((m1 + y(1))^2 + y(2)^2)^(3/2);
+L1 = (2*((m2 - y(1))^2 + y(2)^2)^(5/2));
+L2 = ((m1 + y(1))^2 + y(2)^2)^(5/2);
+
+
+dfy =[ 0,  0,  1, 0;
+       0,  0,  0, 1;
+       (3*m1*(2*m2 - 2*y(1))*(m2 - y(1)))/L1 - D1 - D2 + (3*m2*(2*m1 + 2*y(1))*(m1 + y(1)))/(2*L2) + 1,  (3*m2*y(2)*(m1 + y(1)))/L2 - (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)))/L2 - (3*m1*y(2)*(2*m2 - 2*y(1)))/L1, (3*m1*y(2)^2)/((m2 - y(1))^2 + y(2)^2)^(5/2) - D1 - D2 + (3*m2*y(2)^2)/L2 + 1, -2, 0];
+
+end