Question

6. [20 points] The following nonlinear differential equations describe the motion of a body in orbit around two much heavier bodies. An example would be an Apollo capsule in an Earth-moon orbit. The three bodies determine a two-dimensional Cartesian plane in space The origin is at the center of mass of the two heavy bodies, the r-axis is the line through these two bodies, and the distance between their centers is taken as the unit. Thus, if μ is the ratio of the mass of the moon to that of Earth, then the centers of the moon and the earth are located at coordinates (1-1, 0) and (-1,0), respectively, and the coordinate system moves as the moon rotates about Earth The third body, the Apollo spacecraft, is assumed to have a mass that is negligible compared to the other two bodies, and its position as a function of time is (x(t), y(t)). The equations are derived from Newton's law of motion and the inverse square law of gravitation; the first-order derivatives in the equations come from the rotating coordinate system. The governing equations for motion of the spacecraft are 3 where 82.45 and the initial conditions are r(0)1.2, r00, y(00, and y'(0)1.04935751 The solution will be periodic with a period T = 6.19216933. This means that Apollo starts on the far side of the moon with an altitude about 0.2 times the earth-moon distance and a given initial velocity. The resulting orbit brings Apollo in close proximity to earth, out in a big loop in the opposite side the the earth from the moon, back in close to the earth again, and finally back to its original position and velocity on the far side of the moon (a) Solve the equations ofmotion for the given initial conditions in the domain 0 < t < 2T using the MATLAB ode45 function. Use the following ode45 call in your program: options odeset ('RelTol' 1e-6); [TIME, z] - ode45 ('filename', [0 2*T], [initial conditions], options); and plot Apollo's orbit y(x) for the above conditions. In your plot, clearly mark the position of the moon and the earth, and set the axis properties using axis([-1.5 1.5 -1.5 1.5]); axis equal, grid on Use the same axis properties as in part (a). What happens? can the velocity change and the craft still make it back to (near) Earth? (b) Try half and double the initial velocity and plot the resulting orbits for 0 t 2T. (c) How sensitive is the orbit of Apollo to the initial velocity? In other words, how much

Answer 1

ANSWER:

Here I'm writing the code

Code:

odefunc.m

function dydt = odefunc(t,a)
mu = 1/82.45;
mu_ = 1-mu;
r1 = sqrt((a(1)+mu)^2 + a(2)^2);
r2 = sqrt((a(1)-mu_)^2 + a(2)^2);
dydt = [a(3);
a(4);
2*a(4) + a(1) - mu_*(a(1)+mu)/(r1^3) - mu*(a(1)-mu_)/(r2^3);
-2*a(3) + a(2) - mu_*a(2)/(r1^3) - mu*a(2)/(r2^3)];
end

odefunc_half.m

function dydt = odefunc_half(t,a)
mu = 1/82.45;
mu_ = 1-mu;
r1 = sqrt((a(1)+mu)^2 + a(2)^2);
r2 = sqrt((a(1)-mu_)^2 + a(2)^2);
dydt = [a(3);
a(4);
a(4) + a(1) - mu_*(a(1)+mu)/(r1^3) - mu*(a(1)-mu_)/(r2^3);
-a(3) + a(2) - mu_*a(2)/(r1^3) - mu*a(2)/(r2^3)];
end

odefunc_double.m

function dydt = odefunc_double(t,a)
mu = 1/82.45;
mu_ = 1-mu;
r1 = sqrt((a(1)+mu)^2 + a(2)^2);
r2 = sqrt((a(1)-mu_)^2 + a(2)^2);
dydt = [a(3);
a(4);
4*a(4) + a(1) - mu_*(a(1)+mu)/(r1^3) - mu*(a(1)-mu_)/(r2^3);
-4*a(3) + a(2) - mu_*a(2)/(r1^3) - mu*a(2)/(r2^3)];
end

code.m

close all
clear
clc

T = 6.19216933;
ic = [0;
-1.04935751;
-1.84;
0];
mu = 1/82.45;
mu_ = 1-mu;

options = odeset('RelTol', 1e-6);
[TIME,Z] = ode45(@odefunc, [0,2*T], ic, options);
figure, plot(Z(:,1),Z(:,2)); hold on
scatter(mu_,0), scatter(-mu,0);
legend('show');
legend('Orbit', 'Moon', 'Earth')
axis([-1.5 1.5 -1.5 1.5]); axis equal, grid on;
title('y(x)');

[TIME,Z] = ode45(@odefunc_half, [0,2*T], ic, options);
figure, plot(Z(:,1),Z(:,2)); hold on
scatter(mu_,0), scatter(-mu,0);
legend('show');
legend('Orbit', 'Moon', 'Earth')
grid on;
title('y(x) (Half the initial velocity)');
% Note that the axis is not same as other two plots.

[TIME,Z] = ode45(@odefunc_double, [0,2*T], ic, options);
figure, plot(Z(:,1),Z(:,2)); hold on
scatter(mu_,0), scatter(-mu,0);
legend('show');
legend('Orbit', 'Moon', 'Earth')
axis([-1.5 1.5 -1.5 1.5]); axis equal, grid on;
title('y(x) (Double the initial velocity)');

Resultant Graph Plots:

y(x) 1.5 Orbit O Moon O Earth 0.5 -0.5 -1 21.5-1 0.5 0 0.5 1.5

y(x) (Half the initial velocity) 2000 Orbit O Moon O Earth 0 -2000 4000 6000 -8000 -10000 -12000 14000 16000 -3 -2.5 -2 1.5 -1 -0.5 0.5 x 10

y(x) (Double the initial velocity) 1.5 Orbit Moon O Earth 0.5 -0.5 -1 1.5-1 -0.5 0 0.5 1.5
Note: Still if you have any doubts let me know in the comment box. Thank You. Good Luck>