%% Problem Set D: Problem 6
%% (a)
% Let w(1), w(2) represent the coordinates x1 and y1, w(3), w(4) represent the coordinates x2 and y2, w(5), w(6) represent the components of the velocity of the planet, w(7), w(8) represent the components of the velocity of the comet. First we assume no interactions.

opts = odeset('RelTol',1e-5,'AbsTol',1e-7);
[tsol, wsol] = ode45(@(t, w) [w(5);w(6);w(7);w(8); ...
     -w(1)./(w(1).^2+w(2).^2).^(3/2); ...
     -w(2)./(w(1).^2+w(2).^2).^(3/2); ...
     -w(3)./(w(3).^2+w(4).^2).^(3/2); ...
     -w(4)./(w(3).^2+w(4).^2).^(3/2)], ...
     [0, 10], [1; 0; 2; 0; 0; 1; 0; 0.1], opts);
plot(wsol(:,1), wsol(:,2), 'r'), hold on;
plot(wsol(:,3), wsol(:,4), 'b'), hold off; axis equal
title('Motion of a comet (blue) and a planet (red)')

%
% Let's blow up the picture to make sure the comet doesn't crash into the sun. Here the sun is shown as a red circle.

axis([-.1, .1, -.1, .1]), hold on
plot(0, 0, 'or', 'MarkerSize', 10), hold off
title('Motion of a comet (blue) close to the sun')

%
%% (b)
% Now we put in the interaction term. This affects the equation for the comet by addition of a term involving the vector difference between the positions of the planet and the comet. The equation of motion of the planet is unchanged.

[tsol, wsol] = ode45(@(t, w) [w(5);w(6);w(7);w(8); ...
     -w(1)./(w(1).^2+w(2).^2).^(3/2); ...
     -w(2)./(w(1).^2+w(2).^2).^(3/2); ...
     -w(3)./(w(3).^2+w(4).^2).^(3/2) - ...
  0.1*(w(3)-w(1))./((w(3)-w(1)).^2+(w(4)-w(2)).^2).^(3/2); ...
     -w(4)./(w(3).^2+w(4).^2).^(3/2) - ...
  0.1*(w(4)-w(2))./((w(3)-w(1)).^2+(w(4)-w(2)).^2).^(3/2)], ...
     [0, 50], [1; 0; 2; 0; 0; 1; 0; 0.1], opts);
plot(wsol(:,1), wsol(:,2), 'r'), hold on;
plot(wsol(:,3), wsol(:,4), 'b'), hold off; axis equal
title(['Motion of a comet (blue) and a planet (red) ', ...
  'with interactions'])

% Note that the comet ``precesses" around instead of staying in an elliptical orbit.