clear
clf

rxs=0;     %x coordinate of the Sun's position
rys=0;     %y coordinate of the Sun's position
ryp=0;          %y coordinate of the planet's position
rxp = 20;       %x coordinate of the planet's position
global mSun mPlanet G mComet
mSun    = 1e+7;
mPlanet = 1e+4;
G=1;
mComet = 1;
f=0;

R=20;    %distance from rs_initial to rp_initial
rs_initial=[rxs 0];
rp_initial=[rxp 0];
rcm=(rp_initial.*mPlanet+rs_initial.*mSun)/(mPlanet+mSun); %center of mass

rs=rs_initial-rcm;
rp=rp_initial-rcm;

rc = [(R-norm(rcm))*cosd(45) (R-norm(rcm))*sind(45)];  %position of comet

rxc=rc(1);    %x coordinate of the comet's position
ryc=rc(2);    %y coordinate of the comet's position



m=mSun*mPlanet/(mSun+mPlanet);
omega=sqrt(G*(mSun+mPlanet)/R^3); %angular velocity

vs=[0 -norm(rs)].*omega; %velocity of sun

vp=[0 norm(rp)].*omega;  %velocity of planet

vc=[-R*sind(45)*0.3 R*cosd(45)*0.3].*omega; %velocity of comet


rc2 = [(R-norm(rcm))*cosd(45)+0.3 (R-norm(rcm))*sind(45)+0.3];  %position of second comet

vc2 = [-R*sind(45)*0.35 R*cosd(45)*0.35].*omega; %velocity of second comet

y=[rs(1); rs(2); vs(1); vs(2);
    rp(1);rp(2);vp(1);vp(2);
    rc(1);rc(2);vc(1);vc(2);
    rc2(1);rc2(2);vc2(1);vc2(2)];



T=sqrt(4*(pi)^2*(R)^3/(G*(mSun+mPlanet)));
dt=T/1000;
t=0;
tmax=T/5.2;
n=tmax/dt;
n = int64(n);
tt = zeros(n+1,1);
yt = zeros(n+1,size(y,1));
tt(1)   = t;
yt(1,:) = y';



rc3 = [(R-norm(rcm))*cosd(120)+2 (R-norm(rcm))*sind(120)+2];  %position of third comet

rxc3=rc3(1);    %x coordinate of the comet's position
ryc3=rc3(2);    %y coordinate of the comet's position

vc3=[-R*sind(120)*0.5 R*cosd(120)*0.5].*omega;  %velocity of third comet

rc4 = [(R-norm(rcm))*cosd(180)+3 (R-norm(rcm))*sind(180)+3]; %position of fourth comet

vc4=[-R*sind(180)*0.55 R*cosd(180)*0.55].*omega;  %velocity of fourth comet

y_new=[rs(1); rs(2); vs(1); vs(2);
    rp(1);rp(2);vp(1);vp(2);
    rc3(1);rc3(2);vc3(1);vc3(2);
    rc4(1);rc4(2);vc4(1);vc4(2)];

t=0;
t_newmax=T;
m=t_newmax/dt;
m = int64(m);
tt = zeros(m+1,1);
y_newt = zeros(m+1,size(y_new,1));
tt(1)   = t;
y_newt(1,:) = y_new';


for i=1:m
    if i <=n
        f1 = kepler_ode_khan(t     , y                );
        f2 = kepler_ode_khan(t+dt/2, y+f1 .*  (dt/2)  );
        f3 = kepler_ode_khan(t+dt/2, y+f2 .*  (dt/2)  );
        f4 = kepler_ode_khan(t+dt  , y+f3 .*   dt     );
        
        y = y + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
        t = t+dt;
        tt(i+1)   = t;
        yt(i+1,:) = y';
     
        figure(1)
        plot(0,0,'X',yt(:,1),yt(:,2),'-or',yt(:,9),yt(:,10),'k.',yt(:,13),yt(:,14),'k.');
        xlim([-20 15])
        ylim([-10 25])
        hold on
        pause(0.01)
        
        F1 = kepler_ode_khan2(t     , y_new                );
        F2 = kepler_ode_khan2(t+dt/2, y_new+F1 .*  (dt/2)  );
        F3 = kepler_ode_khan2(t+dt/2, y_new+F2 .*  (dt/2)  );
        F4 = kepler_ode_khan2(t+dt  , y_new+F3 .*   dt     );
        
        y_new = y_new + (F1 + F2.*2 + F3.*2 + F4) .* (dt/6);
        t = t+dt;
        tt(i+1)   = t;
        y_newt(i+1,:) = y_new';
        
        plot(0,0,'X',y_newt(:,1),y_newt(:,2),'-or',y_newt(:,9),y_newt(:,10),'m.',y_newt(:,13),y_newt(:,14),'g.');
        hold on
        pause(0.01)
        if mod(i+9,10) == 0 || i==n || i==m
            f=f+1;
            saveas(gcf, ['~/Desktop/images_comet/comets_',sprintf('%03d',f),'.png'])
            pause(0.01)
        end

    end
    if i > n
        F1 = kepler_ode_khan2(t     , y_new                );
        F2 = kepler_ode_khan2(t+dt/2, y_new+F1 .*  (dt/2)  );
        F3 = kepler_ode_khan2(t+dt/2, y_new+F2 .*  (dt/2)  );
        F4 = kepler_ode_khan2(t+dt  , y_new+F3 .*   dt     );
        
        y_new = y_new + (F1 + F2.*2 + F3.*2 + F4) .* (dt/6);
        t = t+dt;
        tt(i+1)   = t;
        y_newt(i+1,:) = y_new';
        
        plot(0,0,'X',y_newt(:,1),y_newt(:,2),'-or',y_newt(:,9),y_newt(:,10),'m.',y_newt(:,13),y_newt(:,14),'g.');
        hold on
        pause(0.01)
        
        if mod(i+9,10) == 0 
            f=f+1;
            fname = ['~/Desktop/images_comet/comets_',sprintf('%03d',f),'.png'];
            disp(['saving ' fname])
            saveas(gcf, fname)
            %pause(0.01)
        end
        
    end
    
    
end


