% integrates equation of relative motion between two masses m1 and m2. m1
% is the Earth, m2 a spacecraft
clear all
clc
close all
Re = 6378; %km Radius of Earth

r0 = [8000; 0; 6000]; %km 
v0 = [0; 7; 0]; %km/s 
Y0 = [r0; v0]; %initial state vector (km & km/s)
tf = 1*4*3600; %s days
TSPAN = [0 tf]; %duration of integration in seconds
OPTIONS = odeset('Maxstep', 10);

%unperturbed orbit
[TOUT,YOUT] = ode45(@relacc,TSPAN,Y0,OPTIONS);
figure(1)
plot3(YOUT(:,1),YOUT(:,2),YOUT(:,3),'r','linewidth',2)
grid on
hold on
axis equal
plot3(YOUT(1,1),YOUT(1,2),YOUT(1,3),'r*')
xlabel('x (km)'); ylabel('y (km)'); zlabel('z (km)');
[X,Y,Z] = sphere;
surf(6378*X, 6378*Y, 6378*Z)

[min_altitude, min_index] = min(sqrt(YOUT(:,1).^2+YOUT(:,2).^2+YOUT(:,3).^2));
min_altitude = min_altitude - Re
speed_at_min = norm(YOUT(min_index,4:6))

[max_altitude, max_index] = max(sqrt(YOUT(:,1).^2+YOUT(:,2).^2+YOUT(:,3).^2));
max_altitude = max_altitude - Re
speed_at_max = norm(YOUT(max_index,4:6))

% %orbit with J2
% [TOUT_J2,YOUT_J2] = ode45(@relacc_with_j2,TSPAN,Y0,OPTIONS); % 
% plot3(YOUT_J2(:,1),YOUT_J2(:,2),YOUT_J2(:,3),'b')

