% orbit_sim
% See: https://codepen.io/BobsPen/pen/VQgaYL


close all
clear
% Ellipse
a = 200;
b = 120;
e = sqrt(1-b^2/a^2);
c = e*a; % Linear eccentricity

u = 1; % standard grav. param?
T = 2*pi*a^(1.5)/sqrt(u);

dt = T/1000;
timestep = 1/60; % fps

ba = a-0.001; %Boundary for a.
xi = -ba;
vx = 0; % x-component of velocity.
s = 1; % sign of y-coord.
px = 0; % x and y coords on ellipse.
py = 0;


%% Plot ellipse
% Parametric definition
t = linspace(0, 2*pi, 1000);
x = a*cos(t);
y = b*sin(t);

% Plot the ellipse.
plot(x,y,'k');
hold on
grid on
axis equal

%% Plot satellite
h = plot(px, py, 'r*')
h.XDataSource = 'px';
h.YDataSource = 'py';

while(1)
    px = xi + vx*dt;
    if(s*px > ba)
        px = s*ba;
        s = -s;
    end

    py = s*b * sqrt(1 - px^2/a^2);
    vx = py/b * sqrt(u*a/((px - c)^2 +py^2));
    xi = px;

    refreshdata(h)
    
    pause(timestep)
end
