% Program 6.4 Plotting program for three-body problem
% Inputs: time interval inter, initial conditions
% ic = [x1 vx1 y1 vy1 x2 vx2 y2 vy2 x3 vx3 y3 vy3], x1-x3 pos, x1-x3 vel, y1-y3 pos, y1-y3 vel
% m = [m1 m2 m3] (masses of the three bodies)
% number of steps n, steps per point plotted p
% Outputs animated gif to fname (unless fname is ‘’)
% Case 1: orbit3b([0 100],[2 .2 2 -.2 0 0 0 0 -2 -.2 -2 .2], [.03 .30 .03], 50000,500,'case1.gif')
% Case 2: orbit3b([0 200],[2 .20001 2 -.2 0 0 0 0 -2 -.2 -2 .2], [.03 .30 .03], 100000,500,'case2.gif')
% Figure 8: orbit3b([0 50],[-.97 -.466 .243 -.433 .97 -.466 -.243 -.433 0 .932 0 .866],[1 1 1],50000,250, 'fig8_orig.gif')
% Chaotic but controlled: orbit3b([0 200],[2 .2 2 -.2 0 0 0 0 -2 -.1 -2 .1], [.03 .3 .06], 50000, 250, 'controlled.gif')
% Chaotic followed by near-singularity: orbit3b([0 35],[-.97 -.466 .243 -.433 .97 -.466 .243 -.433 0 .932 0 .866],[1 1 1], 17500, 100,'chaos.gif')
function z=orbit3b(inter,ic,m,n,p,fname)
h=(inter(2)-inter(1))/n;    % plot n points
x1=ic(1);vx1=ic(2);y1=ic(3);vy1=ic(4);  % grab initial conds for body 1
x2=ic(5);vx2=ic(6);y2=ic(7);vy2=ic(8);  % grab initial conds for body 2
x3=ic(9);vx3=ic(10);y3=ic(11);vy3=ic(12);   % grab initial conds for body 3
y(1,:)=[x1 vx1 y1 vy1 x2 vx2 y2 vy2 x3 vx3 y3 vy3];t(1)=inter(1);  % build y vector    
set(gca,'XLim',[-5 5],'YLim',[-5 5],'XTick',[-5 0 5],'YTick',...
  [-5 0 5],'Drawmode','fast','Visible','on');
cla;
%sun=line('color','y','Marker','.','markersize',25,...
%  'xdata',0,'ydata',0);
drawnow;
tail1=line('color',[1,0,0],'LineStyle','-','erase','normal',...
  'xdata',[],'ydata',[],'LineWidth',2);
hist1=line('color',[1,0,0],'LineStyle',':','erase','normal',...
  'xdata',[],'ydata',[]);
tail2=line('color',[0,1,0],'LineStyle','-','erase','normal',...
  'xdata',[],'ydata',[],'LineWidth',2);
hist2=line('color',[0,1,0],'LineStyle',':','erase','normal',...
  'xdata',[],'ydata',[]);
tail3=line('color',[0,0,1],'LineStyle','-','erase','normal',...
  'xdata',[],'ydata',[],'LineWidth',2);
hist3=line('color',[0,0,1],'LineStyle',':','erase','normal',...
  'xdata',[],'ydata',[]);
head1=line('color',[1,0,0],'Marker','.','markersize',4 * log(m(1)) + 35,...
  'erase','normal','xdata',[],'ydata',[]);
head2=line('color',[0,1,0],'Marker','.','markersize',4 * log(m(2)) + 35,...
  'erase','normal','xdata',[],'ydata',[]);
head3=line('color',[0,0,1],'Marker','.','markersize',4 * log(m(3)) + 35,...
  'erase','normal','xdata',[],'ydata',[]);

% initialize recency line and full history data
line_len = round(10/h/p);
linedata = zeros(line_len,12);
fullhist = zeros(1,12);
for k=1:line_len
    linedata(k,:) = y(1,:);
end
fullhist(1,:) = y(1,:);

% main loop - simulate for p steps, then sa
for k=1:n/p
  for i=1:(line_len-1)
      linedata(i, :) = linedata(i+1, :);
  end
  fullhist(size(fullhist,1)+1,:)=linedata(1,:);
  for i=1:p
    t(i+1)=t(i)+h;
    y(i+1,:)=rungekutta4step(t(i),y(i,:),h, m); % update method for euler/trapezoid/rungekutta
  end
  y(1,:)=y(p+1,:);t(1)=t(p+1);
  linedata(line_len,:) = y(1,:);
  set(head1,'xdata',y(1,1),'ydata',y(1,3))
  set(head2,'xdata',y(1,5),'ydata',y(1,7))
  set(head3,'xdata',y(1,9),'ydata',y(1,11))
  set(tail1,'xdata',linedata(:,1),'ydata',linedata(:,3))
  set(tail2,'xdata',linedata(:,5),'ydata',linedata(:,7))
  set(tail3,'xdata',linedata(:,9),'ydata',linedata(:,11))
  if k == n/p
    set(hist1,'xdata',fullhist(:,1),'ydata',fullhist(:,3))
    set(hist2,'xdata',fullhist(:,5),'ydata',fullhist(:,7))
    set(hist3,'xdata',fullhist(:,9),'ydata',fullhist(:,11))
  end
  drawnow;
 
  %gif code
  if strcmp(fname, '') == 0
    frm = getframe(1);
    img = frame2im(frm);
    [A,map] = rgb2ind(img,256);
    if k == 1 %first frame
        imwrite(A,map,fname,'gif','LoopCount',Inf,'DelayTime',1);
    elseif k == n/p %final frame
        imwrite(A,map,fname,'gif','WriteMode','append','DelayTime',1);
    else
        imwrite(A,map,fname,'gif','WriteMode','append','DelayTime', h * p / 20);
    end
  end
end

% function y=eulerstep(t,x,h)
% % one step of the Euler method
% y=x+h*ydot(t,x);

% function y=trapezoidstep(t,x,h)
% %one step of trapzeoid method
% y=x+h/2*(ydot(t,x) + ydot(t+h, x+h*ydot(t,x)));

function y=rungekutta4step(t,x,h,m)
% one step of RK4 method
k1 = ydot(t, x, m);
k2 = ydot(t + h/2, x + h/2 * k1, m);
k3 = ydot(t + h/2, x + h/2 * k2, m);
k4 = ydot(t + h, x + h * k3, m);
y = x + h/6 * (k1 + 2*k2 + 2*k3 + k4);

function z = ydot(t,x,m)
m1=m(1);m2=m(2);m3=m(3);
g=1;mg1=m1*g;mg2=m2*g;mg3=m3*g;
px1=x(1);py1=x(3);vx1=x(2);vy1=x(4);
px2=x(5);py2=x(7);vx2=x(6);vy2=x(8);
px3=x(9);py3=x(11);vx3=x(10);vy3=x(12);
dist12=sqrt((px2-px1)^2+(py2-py1)^2);
dist13=sqrt((px3-px1)^2+(py3-py1)^2);
dist23=sqrt((px3-px2)^2+(py3-py2)^2);
z=zeros(1,12);
z(1)  = vx1;
z(2)  = ((mg2*(px2-px1))/(dist12^3))+((mg3*(px3-px1))/(dist13^3));
z(3)  = vy1;
z(4)  = ((mg2*(py2-py1))/(dist12^3))+((mg3*(py3-py1))/(dist13^3));
z(5)  = vx2;
z(6)  = ((mg1*(px1-px2))/(dist12^3))+((mg3*(px3-px2))/(dist23^3));
z(7)  = vy2;
z(8)  = ((mg1*(py1-py2))/(dist12^3))+((mg3*(py3-py2))/(dist23^3));
z(9)  = vx3;
z(10) = ((mg1*(px1-px3))/(dist13^3))+((mg2*(px2-px3))/(dist23^3));
z(11) = vy3;
z(12) = ((mg1*(py1-py3))/(dist13^3))+((mg2*(py2-py3))/(dist23^3));
