function HW3_2022
close all

solveproblem4
solveproblem5

%% Problem 4
function solveproblem4
    % Properties
    eta = 0.001;   % Viscosity, Ns/m^2
    rhoc = 1050;   % Cell density, kg/m^3
    rhof = 1000;   % Fluid viscosity kg/m^3
    Rc = 10.e-06;  % Cell radius, m

    Rcentrifuge = 36.5e-03;  % Radius of centrifuge
    Omega = 1000*2*pi/60;    % Angular speed, rad/s

    CB = (rhof/rhoc);
    CD = 9*eta/(2*rhoc*Rc^2);

    initial_r = 10.e-03;
    time_int = [0,4000];
    initial_w = [initial_r; 0; 0; Omega];

    options = odeset('Event',@(t,w) centrifuge_event(t,w,Rcentrifuge));

    [times,sols] = ode15s(@(t,w) centrifuge(t,w,CB,CD,Omega),time_int,initial_w,options);

    % Plot the radial velocity
    figure
    plot(times,sols(:,3))
    xlabel('time (s)')
    ylabel('Radial Velocity (mm/s)')
    title('Problem 4.4 (a) ')
    axes1 = gca;
    set(axes1,'FontName','Times New Roman','FontSize',14,'XGrid','on','YGrid',...
        'on');

    % Plot the difference in angular velocity of cell and centrifuge
    figure
    plot(times,sols(:,4)-Omega)
    xlabel('time (s)')
    ylabel('Angular Velocity Difference (rad/s)')
    title('Problem 4.4 (b) ')
    axes1 = gca;
    set(axes1,'FontName','Times New Roman','FontSize',14,'XGrid','on','YGrid',...
        'on');

    % The approximate solution
    rapprox = initial_r*exp( (1-CB)*Omega^2.*times/CD );

    % Plot radial position and compare with approximate solution
    figure
    p1 = plot(times,sols(:,1));
    hold on
    mask = rapprox < Rcentrifuge;
    p2 = plot(times(mask),rapprox(mask));
    xlabel('time (s)')
    ylabel('Radial position (mm)')
    title('Problem 4.6')
    axes1 = gca;
    set(axes1,'FontName','Times New Roman','FontSize',14,'XGrid','on','YGrid',...
        'on');
    legend1 = legend([p1,p2],'Exact Solution','Approximate Solution');
    set(legend1,...
        'Position',[0.21665983606557377 0.708333329529014 0.229166666666668 0.187022895474015],...
        'FontAngle','italic',...
        'FontSize',14);

    % Obtain a new solution with a larger cell
    Rc = 100.e-06;
    CD = 9*eta/(2*rhoc*Rc^2);


    [times,sols] = ode15s(@(t,w) centrifuge(t,w,CB,CD,Omega),time_int,initial_w,options);
    rapprox = initial_r*exp( (1-CB)*Omega^2.*times/CD );
    figure
    p1 = plot(times,sols(:,1));
    hold on
    mask = rapprox < Rcentrifuge;
    p2 = plot(times(mask),rapprox(mask));
    xlabel('time (s)')
    ylabel('Radial position (mm)')
    title('Problem 4.7')
    axes1 = gca;
    set(axes1,'FontName','Times New Roman','FontSize',14,'XGrid','on','YGrid',...
        'on');
    legend1 = legend([p1,p2],'Exact Solution','Approximate Solution');
    set(legend1,...
        'Position',[0.21665983606557377 0.708333329529014 0.229166666666668 0.187022895474015],...
        'FontAngle','italic',...
        'FontSize',14);

end

%% Problem 5
function solveproblem5
    G = 6.67e-11;
    m1 = 7.348e22; m2=m1; m3=m1;
    d = 384.4e06;
    v0 = sqrt(5*G*m1/(4*d));
    T = 2*pi*d/v0;
    initial_w = [-d;0;0;0;d;0;0;v0;0;0;0;-v0];
    time_int = [0,4*T];

    % 5.3
    [times,sols] = ode45(@(t,w) threebody(t,w,G,m1,m2,m3),time_int,initial_w);
    
    figure
    p1 = plot(sols(:,1),sols(:,2));
    hold on
    p2 = plot(sols(:,5),sols(:,6));   
    xlabel('x (m)')
    ylabel('y (m)')
    title('Problem 5.3')
    axes1 = gca;
    set(axes1,'FontName','Times New Roman','FontSize',14,'XGrid','on','YGrid',...
        'on');
    legend1 = legend([p1,p2],'Planet 1','Planet 2');
    set(legend1,...
        'Position',[0.31665983606557377 0.308333329529014 0.229166666666668 0.187022895474015],...
        'FontAngle','italic',...
        'FontSize',14);
    
    
    %5.4
    options = odeset('RelTol',0.000001);
    [times,sols] = ode45(@(t,w) threebody(t,w,G,m1,m2,m3),time_int,initial_w,options);

    figure
    p1 = plot(sols(:,1),sols(:,2));
    hold on
    p2 = plot(sols(:,5),sols(:,6));   
    xlabel('x (m)')
    ylabel('y (m)')
    title('Problem 5.4')
    axes1 = gca;
    set(axes1,'FontName','Times New Roman','FontSize',14,'XGrid','on','YGrid',...
        'on');
    legend1 = legend([p1,p2],'Planet 1','Planet 2');
    set(legend1,...
        'Position',[0.31665983606557377 0.308333329529014 0.229166666666668 0.187022895474015],...
        'FontAngle','italic',...
        'FontSize',14);

     initial_w = [-d;0;0;0;d;0;0;v0;0;0;0;-v0];
         time_int = [0:T/50:4*T];
       [times,sols] = ode45(@(t,w) threebody(t,w,G,m1,m2,m3),time_int,initial_w,options);

    figure
    p1 = plot(sols(:,1),sols(:,2));
    hold on
    p2 = plot(sols(:,5),sols(:,6));   
    xlabel('x (m)')
    ylabel('y (m)')
    title('Problem 5.4')
    axes1 = gca;
    set(axes1,'FontName','Times New Roman','FontSize',14,'XGrid','on','YGrid',...
        'on');
    legend1 = legend([p1,p2],'Planet 1','Planet 2');
    set(legend1,...
        'Position',[0.31665983606557377 0.308333329529014 0.229166666666668 0.187022895474015],...
        'FontAngle','italic',...
        'FontSize',14);
     
    % Problem 5.5
     initial_w = [-d/2;-d;d/10;d;d;-d/2;0;0;0;0;0;0];
     time_int = [0:T/500:40*T];
     [times,sols] = ode45(@(t,w) threebody(t,w,G,m1,m2,m3),time_int,initial_w,options);

    figure
    p1 = plot(sols(:,1),sols(:,2));
    hold on
    p2 = plot(sols(:,3),sols(:,4));   
    p3 = plot(sols(:,5),sols(:,6));   
    xlabel('x (m)')
    ylabel('y (m)')
    title('Problem 5.5')
    axes1 = gca;
    set(axes1,'FontName','Times New Roman','FontSize',14,'XGrid','on','YGrid',...
        'on');
    legend1 = legend([p1,p2,p3],'Planet 1','Planet 2','Planet 3');
    set(legend1,...
        'Position',[0.31665983606557377 0.208333329529014 0.229166666666668 0.187022895474015],...
        'FontAngle','italic',...
        'FontSize',14);   
    
    
    %animate_planets(times,sols)
end


end

function dwdt = centrifuge(t,w,CB,CD,Omega)
  r = w(1); theta = w(2); vr = w(3); omega = w(4);
  
  drdt = vr;  dthetadt = omega;
  dvrdt = r*omega^2 - CB*Omega^2*r - CD*vr;
  domegadt = CD*(Omega-omega) - 2*vr*omega/r;

  dwdt = [drdt;dthetadt;dvrdt;domegadt];
  
end

function [ev,stop,dir] = centrifuge_event(t,w,Rcentrifuge)
  r = w(1);
  ev = Rcentrifuge - r;
  dir = -1;
  stop = 1;
end

function dwdt = threebody(t,w,G,m1,m2,m3)
 allvels = w(7:12);
 p1 = w(1:2); p2 = w(3:4); p3 = w(5:6);  % Position vectors
 
 R12 = norm(p1-p2);  R13 = norm(p1-p3); R23 = norm(p2-p3);
 
 F1 = (G*m1*m2/R12^3)*(p2-p1) + (G*m1*m3/R13^3)*(p3-p1);
 F2 = (G*m1*m2/R12^3)*(p1-p2) + (G*m2*m3/R23^3)*(p3-p2);
 F3 = (G*m3*m1/R13^3)*(p1-p3) + (G*m2*m3/R23^3)*(p2-p3);
 
 dwdt = [allvels;F1/m1;F2/m2;F3/m3];

end

function animate_planets(t,sols)

scrsize = get(0,'ScreenSize');
nn = min(scrsize(3:4));
x0 = scrsize(3)/2-nn/3;
y0 = scrsize(4)/3;
figure1 = figure('Position',[x0,y0,2*nn/3,2*nn/3]);

xmax = max(max(max(sols(:,5)),max(sols(:,3))),max(sols(:,1)));
xmin = min(min(min(sols(:,5)),min(sols(:,3))),min(sols(:,1)));
ymax = max(max(max(sols(:,6)),max(sols(:,4))),max(sols(:,2)));
ymin = min(min(min(min(sols(:,6)),sols(:,4))),min(sols(:,2)));

dx = 0.1*(xmax-xmin)+0.5;
dy = 0.1*(ymax-ymin)+0.5;

for i = 1:length(t)
    clf
    axes1 = gca;

    plot(axes1,sols(1:i,5),sols(1:i,6),'LineWidth',2,'LineStyle',...
    ':','Color',[0 0 1]);
    ylim(axes1,[ymin-dy ymax+dy]);
    xlim(axes1,[xmin-dx xmax+dx]);
    hold on
    plot(axes1,sols(1:i,1),sols(1:i,2),'LineWidth',2,'LineStyle',':',...
    'Color',[1 0 0]);
    plot(axes1,sols(1:i,3),sols(1:i,4),'LineWidth',2,'LineStyle',':',...
    'Color',[0 1 0]);
    plot(axes1,sols(i,1),sols(i,2),'ro','MarkerSize',11,'MarkerFaceColor','r');
    plot(axes1,sols(i,3),sols(i,4),'ro','MarkerSize',11,'MarkerFaceColor','g');
    plot(axes1,sols(i,5),sols(i,6),'ro','MarkerSize',11,'MarkerFaceColor','b');

    pause(0.1);   
    
end

end
