clc
clear all
global M Re J2 %for state space
Rvektor0=1000*[-2.3689 4.2406 6.502]; %km
Vvektor0=[-6.8283,-1.5083,-0.428];%km/sec
M=398603;%km^3/sec^2
Re=6378;%km
J2=1.082*10^(-3);
R0=norm(Rvektor0);%|R|
V0=norm(Vvektor0);%|V|
Tamount=864e3; % Simulation duration [sec]
t=0:8.64:864000;

%Q1a-

opts = odeset('RelTol',1e-8,'AbsTol',1e-10); % Acceptable error
[t1,Rvektor]=ode45(@tr3statespace,t,[Rvektor0 Vvektor0]',opts);
rvektor=[Rvektor(:,1),Rvektor(:,2),Rvektor(:,3)]; %the position vector in time [km]
vvektor=[Rvektor(:,4),Rvektor(:,5),Rvektor(:,6)]; %the velocity vector in time [km/sec]
for i=1:(length(t)/10)
r(i)=norm(rvektor(i,:)); %radius value in time km
v(i)=norm(vvektor(i,:)); %velocity value in time km/sec
epsilon(i)=((v(i))^2)/2-M/r(i); %specific energy km^2/sec^2
hv(i,:)=cross(rvektor(i,:),vvektor(i,:)); %Specific Angular Momentum Vector km^2/sec
h(i)=norm(hv(i,:)); %Specific Angular Momentum km^2/sec
%Orbit Angle [rad]:
phi(i)=(acos(h(i)/(r(i)*v(i))));
if((rvektor(i,:)*(vvektor(i,:))')<0)
   phi(i)=-phi(i);
end
t_sug(i)=t(i);
end
figure(1)
p=plot3(rvektor(:,1),rvektor(:,2),rvektor(:,3))
p.LineWidth = 1;
pbaspect([1 1 1]);
title('The Satellite movement around earth')
xlabel('x [km]');
ylabel('y [km]');
zlabel('z [km]');

figure(2)
plot(t_sug,epsilon)
title('Satellite Specific Energy - \epsilon(t)[km^2/sec^2]')
ylabel('epsilon [km^2/sec^2]');
xlabel('time [sec]');
grid on
figure(3)
plot(t_sug,h)
title('Satellite Specific Angular Momentum  h(t)[km^2/sec]')
ylabel('h [km^2/sec]');
xlabel('time [sec]');
grid on
figure(4)
plot(t_sug,r)
title('Satellite Motion Radius  r(t) [km]')
ylabel('r [km]');
xlabel('time [sec]');
grid on
figure(5)
plot(t_sug,v)
title('Satellite Motion Velocity  v(t) [km/sec]')
ylabel('v [km/sec]');
xlabel('time [sec]');
grid on
figure(6)
plot(t_sug,phi)
title('Satellite Orbit Angle - \phi (t) [rad]')
ylabel('phi [rad]');
xlabel('time [sec]');
grid on

%Q1b-

opts = odeset('RelTol',1e-8,'AbsTol',1e-10); % Acceptable error
[t,Rvektor2]=ode45(@tr3statespacenonsphere,t,[Rvektor0 Vvektor0]',opts);
rvektor2=[Rvektor2(:,1),Rvektor2(:,2),Rvektor2(:,3)]; %the position vector in time [km]
vvektor2=[Rvektor2(:,4),Rvektor2(:,5),Rvektor2(:,6)]; %the velocity value in time [km/sec]
for i=1:(length(t)/10)
r2(i)=norm(rvektor2(i,:)); %radius value in time km
v2(i)=norm(vvektor2(i,:)); %velocity value in time km/sec
epsilon2(i)=((v2(i))^2)/2-M/r2(i); %specific energy km^2/sec^2
hv2(i,:)=cross(rvektor2(i,:),vvektor2(i,:)); %Specific Angular Momentum Vector km^2/sec
h2(i)=norm(hv2(i,:)); %Specific Angular Momentum km^2/sec
%Orbit Angle [rad]:
phi2(i)=(acos(h2(i)/(r2(i)*v2(i))));
if((rvektor2(i,:)*(vvektor2(i,:))')<0)
   phi2(i)=-phi2(i);
end
t_sug2(i)=t(i);
end

figure(7)
p=plot3(rvektor2(:,1),rvektor2(:,2),rvektor2(:,3))
p.LineWidth = 1;
pbaspect([1 1 1]);
title('The Satellite movement around earth')
xlabel('x [km]');
ylabel('y [km]');
zlabel('z [km]')

figure(8)
plot(t_sug2,epsilon2)
title('Satellite Specific Energy - \epsilon(t)[km^2/sec^2]')
ylabel('epsilon [km^2/sec^2]');
xlabel('time [sec]');
grid on

figure(9)
plot(t_sug2,h2)
title('Satellite Specific Angular Momentum  h(t)[km^2/sec]')
ylabel('h [km^2/sec]');
xlabel('time [sec]');
grid on
figure(10)
plot(t_sug2,r2)
title('Satellite Motion Radius  r(t) [km]')
ylabel('r [km]');
xlabel('time [sec]');
grid on
figure(11)
plot(t_sug2,v2)
title('Satellite Motion Velocity  v(t) [km/sec]')
ylabel('v [km/sec]');
xlabel('time [sec]');
grid on
figure(12)
plot(t_sug2,phi2)
title('Satellite Orbit Angle  \phi (t) [rad]')
ylabel('\phi [rad]');
xlabel('time [sec]');
grid on