function elliptical()
clear, close all,
% +-----------------------------------------------+ %
% |------------------ Constants ------------------| %
% +-----------------------------------------------+ %
global dt;
global points;
global G;
global m_Saturn;
global r_Saturn;
global m_Prometheus;
global s_Prometheus;
global v_Prometheus;
global m_Pan;
global s_Pan;
global v_Pan;
%Constants(Saturn)
dt	= 1000;
points = 900;
G = 6.674*(10^-11);	        % (m^3)/(kg*(s^2))
m_Saturn = 5.683*(10^26);   % kg
r_Saturn = 58232000; 	    % m
%Constants(Shepard Moons)
m_Prometheus = 1.595*10^17;%7.347*10.^22;
s_Prometheus = 139380000;
v_Prometheus = 16540;
m_Pan = 4.95*10^15;
s_Pan = 133584000;
v_Pan = 16900;
%Constants(Test Particle(s))
distSat = 81248000;
mass = 14000;
velocity = 15000;
% +-----------------------------------------------+ %
% |--------------- Initial Values  ---------------| %
% +-----------------------------------------------+ %
%Initial values of the first test particle
pos_test1_x	= s_Prometheus+1000;% m
pos_test1_y	= -10000; %-distSat;		    % m
vel_test1_x = 0;				% m/s
vel_test1_y	= velocity;			% m/s
m_test1 	= mass;				% kg
%Initial values of the second test particle
pos_test2_x	= s_Pan+1000;   % m
pos_test2_y	= -10000; %-distSat;		    % m
vel_test2_x = 0;				% m/s
vel_test2_y	= velocity;			% m/s
m_test2 	= mass;				% kg
%Initial values of the Satellite Prometheus
distPro = s_Prometheus;
pos_Pro_x = s_Prometheus;	    % m
pos_Pro_y = 0;					% m
vel_Pro_x = 0;					% m/s
vel_Pro_y = v_Prometheus;		% m/s
m_Pro = m_Prometheus;			% kg
%Initial values of the Satellite Pandora
distPan = s_Pan;
pos_Pan_x = -s_Pan;
pos_Pan_y = 0;
vel_Pan_x = 0;
vel_Pan_y = -v_Pan;
m_Pan = m_Pan;
% +-----------------------------------------------+ %
% |--------------- Start the graph ---------------| %
% +-----------------------------------------------+ %
myvideo = VideoWriter('myvid.avi');
%opens figure to record video
open(myvideo);
figure(900);
figure('Units','normalized','Position',[0 0 1 1]);

%figure
xlabel('pos_sat_x [m]');
ylabel('pos_sat_y [m]');
title(['Orbit Analysis']);
% Draw the planet Saturn %
for point=0:361
	hold on
	x = r_Saturn*cosd(point);
	y = r_Saturn*sind(point);
	plot(x,y,'.')
end
text(-4000,0,'Saturn')
plot(250000000,250000000,'.');
plot(250000000,-250000000,'.');
plot(-250000000,250000000,'.');
plot(-250000000,-250000000,'.');
% +-----------------------------------------------+ %
% |---------------  Draw the orbit ---------------| %
% +-----------------------------------------------+ %
for point=0:points
%Changing position and velocity of the satellite Prometheus
    %hold on
	%plot(pos_Pro_x,pos_Pro_y,'r.')
    t= text(pos_Pro_x,pos_Pro_y,'.','Color','red','FontSize',30);
    set(t,'Visible','on')
    pause(0.05)
    
    [vel_Pro_x, vel_Pro_y, pos_Pro_x, pos_Pro_y] = MoonAccel(vel_Pro_x, vel_Pro_y, pos_Pro_x, pos_Pro_y);
%Changing position and velocity of the satellite Pandora
    %hold on
	%plot(pos_Pro_x,pos_Pro_y,'r.')
    s= text(pos_Pan_x,pos_Pan_y,'.','Color','black','FontSize',30);
    set(s,'Visible','on')

    [vel_Pan_x, vel_Pan_y, pos_Pan_x, pos_Pan_y] = MoonAccel(vel_Pan_x, vel_Pan_y, pos_Pan_x, pos_Pan_y);
%Changing position and velocity of the first test particle	
    plot(pos_test1_x,pos_test1_y,'b.')
    [vel_test1_x, vel_test1_y, pos_test1_x, pos_test1_y] = acceleration(vel_test1_x, vel_test1_y, pos_test1_x, pos_test1_y, pos_Pro_x, pos_Pro_y, pos_Pan_x, pos_Pan_y);
	%plot(pos_test_x,pos_test_y,'b.')
%Changing position and velocity of the second test particle	  
    plot(pos_test2_x,pos_test2_y,'m.')
    [vel_test2_x, vel_test2_y, pos_test2_x, pos_test2_y] = acceleration(vel_test2_x, vel_test2_y, pos_test2_x, pos_test2_y, pos_Pro_x, pos_Pro_y, pos_Pan_x, pos_Pan_y);
	%plot(pos_test_x,pos_test_y,'b.')    

	pause(0.000000001)
    writeVideo(myvideo, getframe(gcf));
    set(t,'Visible','off')
    set(s,'Visible','off')
end
set(t,'Visible','on')
set(s,'Visible','on')

close(900);
close(myvideo);
end

function [vel_test_xf, vel_test_yf, pos_test_xf, pos_test_yf] = acceleration(vel_test_xi, vel_test_yi, pos_test_xi, pos_test_yi, pos_Pro_x, pos_Pro_y, pos_Pan_x, pos_Pan_y)
global dt;
global points;
global G;
global m_Saturn;
global r_Saturn;
global m_Prometheus;
global s_Prometheus;
global v_Prometheus;
global m_Pan;
global s_Pan;
global v_Pan;

%function
pos_test_xf = pos_test_xi + vel_test_xi*dt;
pos_test_yf = pos_test_yi + vel_test_yi*dt;
distSat = sqrt(pos_test_yf^2+pos_test_xf^2);
if (distSat<r_Saturn)
    text(pos_test_xf,pos_test_yf,'X');
    return;
    %break;
end
%acceleration from Saturn
gravSat  = -((G*(m_Saturn))/(distSat^2));
angle_Sat = atan(pos_test_yf/pos_test_xf);
accSaty = gravSat*abs(sin(angle_Sat));
accSatx = gravSat*abs(cos(angle_Sat));
if(pos_test_xf<0)
    accSatx=-accSatx;
end
if(pos_test_yf<0)
    accSaty=-accSaty;
end
%acceleration from Prometheus
pos_testPro_x = pos_test_xf - pos_Pro_x;
pos_testPro_y = pos_test_yf - pos_Pro_y;
distance_Pro = sqrt(pos_testPro_y^2+pos_testPro_x^2);
gravPro = -((G*(m_Prometheus))/(distance_Pro^2));
angle_Pro = atan(pos_testPro_y/pos_testPro_x);
accProy = gravPro*abs(sin(angle_Pro));
accProx = gravPro*abs(cos(angle_Pro));
if(pos_testPro_x<0)
    accProx=-accProx;
end
if((pos_testPro_y<0))
    accProy=-accProy;
end
%acceleration from Pandora
pos_testPan_x = pos_test_xf - pos_Pan_x;
pos_testPan_y = pos_test_yf - pos_Pan_y;
distance_Pan = sqrt(pos_testPan_y^2+pos_testPan_x^2);
gravPan = -((G*(m_Prometheus))/(distance_Pan^2));
angle_Pan = atan(pos_testPan_y/pos_testPan_x);
accPany = gravPan*abs(sin(angle_Pan));
accPanx = gravPan*abs(cos(angle_Pan));
if(pos_testPan_x<0)
    accPanx=-accPanx;
end
if((pos_testPan_y<0))
    accPany=-accPany;
end
%total acceleration
acctotx = accSatx + accProx + accPanx;
acctoty = accSaty + accProy + accPany;
vel_test_xf  = vel_test_xi + acctotx*dt;
vel_test_yf  = vel_test_yi + acctoty*dt;
end

function [vel_moon_xf, vel_moon_yf, pos_moon_xf, pos_moon_yf] = MoonAccel(vel_moon_xi, vel_moon_yi, pos_moon_xi, pos_moon_yi)
global dt;
global points;
global G;
global m_Saturn;
global r_Saturn;
global m_Prometheus;
global s_Prometheus;
global v_Prometheus;
global m_Pan;
global s_Pan;
global v_Pan;

pos_moon_xf = pos_moon_xi + vel_moon_xi*dt;
pos_moon_yf = pos_moon_yi + vel_moon_yi*dt;
distPro = sqrt(pos_moon_yf^2+pos_moon_xf^2);
if (distPro<r_Saturn)
    text(pos_moon_xf,pos_moon_yf,'X');
    return;
end
gravity  = -((G*(m_Saturn))/(distPro^2));
angle = atan(pos_moon_yf/pos_moon_xf);
accy = gravity*abs(sin(angle));
accx = gravity*abs(cos(angle));
if(pos_moon_xf<0)
    accx=-accx;
end
if(pos_moon_yf<0)
    accy=-accy;
end
vel_moon_xf  = vel_moon_xi + accx*dt;
vel_moon_yf  = vel_moon_yi + accy*dt;
end