% ’†�S—Í‚ª�ì—p‚·‚éŽ¿“_

global a;
a = 5.00;
m = 1.00;

timeinterval=[0,100]; % ‰Â•ÏƒXƒeƒbƒv
q0=[2; 0; 1; 0.20]; % shifted ellipsoid
q0=[2; 0; 1; 0.40]; % T0+U0 < 0 elllipsoid
q0=[2; 0; -1; 0.30]; % double elllipsoid

%q0=[2; 0; 1; 1.00]; % T0+U0 = 0
%q0=[2; 0; 0.5; 1.09]; % T0+U0
%q0=[2; 0; 1; 1.20]; % T0+U0 > 0
%q0=[2; 0; 0; 1.00]; % T0+U0 < 0

r0=q0(1); theta0=q0(2); v0=q0(3); omega0=q0(4);
T0 = (1/2)*m*(v0*v0+r0*r0*omega0*omega0);
U0 = -a/r0;
E0= T0 + U0;
E0

[time,q]=ode45(@central_force,timeinterval,q0);
%[time,q]=ode23t(@central_force,timeinterval,q0);
%plot(time,q(:,1),time,q(:,2),'--');
%xlabel('time');
%pause;

set(0,'defaultAxesFontSize',16);
set(0,'defaultTextFontSize',16);

hold on;
x = q(:,1).*cos(q(:,2));
y = q(:,1).*sin(q(:,2));
x0 = r0*cos(theta0);
y0 = r0*sin(theta0);
pbaspect([1 1 1]);

set(0,'defaultAxesFontSize',16);
set(0,'defaultTextFontSize',16);

plot(x,y,'-k');
plot(x(1),y(1),'ok');
quiver(x0, y0, v0, r0*omega0, 'r');
xlabel('x');
ylabel('y');
epsilon = 1e-2;
maxx = ceil(max(x)+epsilon); minx = floor(min(x)-epsilon);
maxy = ceil(max(y)+epsilon); miny = floor(min(y)-epsilon);
if maxx - minx < maxy - miny
    h = (maxy - miny) - (maxx - minx);
    maxx = maxx + h/2;
    minx = minx - h/2;
else
    h = (maxx - minx) - (maxy - miny);
    maxy = maxy + h/2;
    miny = miny - h/2;
end
xlim([minx,maxx]);
ylim([miny,maxy]);
pbaspect([1 1 1]);
filename = strcat('central-force-', num2str(v0), '-', num2str(omega0), '.eps');
saveas(gcf, filename, 'eps');
hold off;

%pause;
%
%angular_momentum = m*(q(:,1).^2).*q(:,4);
%plot(time,angular_momentum,'-k');
