
Re = 6371e3; % Earth radius (m)
g = 9.81; % Earth surface gravity (m/s^2)

% Initial conditions (LEO)
r0 = [1.05*Re; 0];
v0 = [0; sqrt(Re*g)];
q0 = [r0; v0];
t0 = 0.0;

% First integration (almost one orbit)
t1 = 90 * 60;
Npoints = 1000;
t01 = linspace(t0, t1, Npoints);
q01 = lsode(@earth, q0, t01);

% Plot
figure(1)
x01 = q01(:,1);
y01 = q01(:,2);
plot(x01, y01)
axis equal

% Second integration (boost to elliptical orbit)
t2 = t1 + 2.4*t1;
t12 = linspace(t1, t2, Npoints);
r1 = q01(end,1:2)';
v1 = q01(end,3:4)';
delta_vee = 2000;
v1b = v1 * (1 + delta_vee/norm(v1)); % post-burn
q1b = [r1; v1b];
q12 = lsode(@earth, q1b, t12);

% Plot
figure(2)
hold off
plot(x01, y01, 'b')
hold on
x12 = q12(:,1);
y12 = q12(:,2);
plot(x12, y12, 'r')
axis equal

% Third integration (circularize orbit)
t3 = t2 + 9*t1;
t23 = linspace(t2, t3, Npoints);
r2 = q12(end,1:2)';
v2 = q12(end,3:4)';
delta_vee = 1400;
v2b = v2 * (1 + delta_vee/norm(v2)); % post-burn
q2b = [r2; v2b];
q23 = lsode(@earth, q2b, t23);

% Plot
figure(3)
hold off
plot(x01, y01, 'b')
hold on
plot(x12, y12, 'r')
x23 = q23(:,1);
y23 = q23(:,2);
plot(x23, y23, 'k')
axis equal
