clf
clear
global ms G mp ma
N=1000;

ms=10^7;
G=1;
mp=10^4;
ma=10^7;
x=zeros(1,N+1);
y=zeros(1,N+1);
vy=zeros(1,N+1);
vx=zeros(1,N+1);
x(1)=1;
y(1)=0;
vy(1)=2;
vx(1)=0;
r=[x,y];
v=[vx,vy];
% ybar=[x,y,vx,vy];
%Initial Conditions
xs=0;
xp=20;

rs=[xs,0];
rp=[xp,0];
R=norm(rs-rp);

rs=[rs(1)+6,0];

astroid=[norm(rs)*cosd(90), norm(rs)*sind(90)];
nue=(ms*ma)/(ms+ma);
W= sqrt((G*ms*ma)./((R.^3)*nue));
vs=[0,W*-norm(rs)];
va=[W*-norm(astroid)*cosd(90), W*norm(astroid)*sind(90)];
%ybar for the function
ybar=[rs(1),rs(2),vs(1),vs(2),astroid(1), astroid(2),va(1),va(2)];
%T and dt for the time step
T=sqrt((4*((pi)^2)*(R^3))/(G*(ms+ma)));
dt=T/3500;
n=ceil((T/dt));

f=0;
startFrame = 900;

hold on
for j=1:   n/2.33
    f1 = run_homework11(j*dt,ybar(j,:));
    f2 = run_homework11((j*dt)+ (dt/2), ybar(j,:)+f1 .* (dt/2));
    f3 = run_homework11((j*dt)+ (dt/2), ybar(j,:)+f2 .*  (dt/2));
    f4 = run_homework11((j*dt)+ dt, ybar(j,:)+f3 .* dt );
    ybar(j+1,:)= ybar(j,:)+ (dt/6).*(f1+2.*f2 + 2.*f3+f4);
    if (j >= startFrame)
        plot3(ybar(j,1),ybar(j,2),ybar(j,5)*0,'k-*')
        % plot(ybar(:,5),ybar(:,6),'k*')
        plot3(ybar(j,5),ybar(j,6), ybar(j,5)*0,'r*-')
        drawnow
        f=f+1
        saveas(gcf, ['final_', num2str(f, '%04d'), '.png'])
    end
    
end




