


clear
clf;
global ms mp G
ms=1.65087*(10^30);
mp=1.009268*(10^25);

G=6.67384*(10^(-11));

whitebg('black')
xs=0;  % in the future, xs could equal 0-rs
ys=0;
%initial position of planet
xp=1495978707;
yp=0;

rs=[xs,ys];
rp=[xp,yp];
R=1495978707; % distance between Sun and Planet


%center of mass
rcm=((rp.*mp+rs.*ms)/(mp+ms));



wsp=sqrt((G*(ms+mp))/(norm(rs-rp)^3));
w=wsp;

rs=rs-rcm;
vsx=0;
vsy=w*rs(1);

rp=rp-rcm;
vpx=0;
vpy=w*rp(1);

y=[rs vsx vsy rp vpx vpy]';




tt = zeros(1000+1,1);
yt = zeros(1000+1,size(y,1));
for i=1:1000
    rs=[y(1) y(2)];
    
    
    rp=[y(5) y(6)];
    
    
    T=sqrt((4*pi^2)*(norm(rs-rp)^3)/(G*(ms+mp )));
    t    =  0.0;
    tmax = 200*T;
    dt=T/10;
    
    tt(1)   = t;
    yt(1,:) = y';
    
    
    f1 = kepler_ode(t,y);
    f2 = kepler_ode(t+dt/2, y+f1 .*  (dt/2)  );
    f3 = kepler_ode(t+dt/2, y+f2 .*  (dt/2)  );
    f4 = kepler_ode(t+dt  , y+f3 .*   dt     );
    
    y = y + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    
    t = t+dt;
    tt(i+1)   = t;  % this records the time at each step
    yt(i+1,:) = y'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
    
    
    plot(yt(i,5),yt(i,6),'-w',yt(i,1),yt(i,2),'or','MarkerSize',200,'MarkerEdgeColor','k','MarkerFaceColor','r')
    %plot(yt(i,5),yt(i,6),'ow','MarkerSize',6,'MarkerEdgeColor','w','MarkerFaceColor','w');
    axis([-0.0006*10^9 0.0006*10^9 -0.0006*10^9 0.0006*10^9 ])
    f=1;
        if mod(i,5)==0
            f=f+1;
             saveas(gcf, ['kepler78b',sprintf('%03d',i),'.png'])
    
         end
    
    
    hold on
    
    drawnow;
end


title('Kepler 78b trajectory ','fontsize',20);
%axis([-1*10^9 1*10^9 -1*10^9 1*10^9 ])
axis equal
 
