ml98.m

%% Constants
segments = ["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"] % replace 'j' with the sub-assignment letter
                 % lowercase only
%% A
% sun
P = [0, 0];
M = 10;

% planet
p = [1, 0];
v = [0, 2];

% solve the planetary motions using ode45
n = 200;
% time_range = linspace(0,1,n);
time_range = [0,1];


theta0 = [p(1), p(2), v(1), v(2)];
[ts, thetas] = ode45(@diffeq_a, time_range, theta0);


% plotting
if ismember("a", segments)
    hold on
    axis([-1, 1, -0.8, 0.8])

    % plot sun
    plot(P(1), P(2),"Marker", "O")

    cpos = plot(thetas(:,1),thetas(:,2),"Marker", "O", "Color","red");


    for i=1:length(ts(:,1))
        plot(thetas(1:i,1),thetas(1:i,2))
        set(cpos,"XData",thetas(i,1),"YData",thetas(i,2))
        pause(0.1)
    end
end


%% B
% solve the planetary motions using ode45
n = 200;
% time_range = linspace(0,1,n);
time_range = [0,1];


theta0 = [p(1), p(2), v(1), v(2)];
events = odeset('Events', @event); % add the callback
[ts, thetas] = ode45(@diffeq_a, time_range, theta0, events);

% plotting
if ismember("b", segments)
    clf; hold on
    axis([-1, 1, -0.8, 0.8])

    % plot sun
    plot(P(1), P(2),"Marker", "O")

    cpos = plot(thetas(:,1),thetas(:,2),"Marker", "O", "Color","red");

    for i=1:length(ts(:,1))
        plot(thetas(1:i,1),thetas(1:i,2))
        set(cpos,"XData",thetas(i,1),"YData",thetas(i,2))
        pause(0.1)
    end
end

%% C

time_range = [0,1];

p2 = [-1,0];
v2 = [0,-2];

theta0 = [p(1), p(2), v(1), v(2), p2(1), p2(2), v2(1), v2(2)];
[ts, thetas] = ode45(@diffeq_c, time_range, theta0);


% plotting
if ismember("c", segments)
    clf; hold on
    axis([-1, 1, -0.8, 0.8])

    % plot sun
    plot(P(1), P(2),"Marker", "O")

    cpos1 = plot(thetas(:,1),thetas(:,2),"Marker", "O", "Color","red");
    cpos2 = plot(thetas(:,5),thetas(:,6),"Marker", "O", "Color","green");

    for i=1:length(ts(:,1))
        plot(thetas(1:i,1),thetas(1:i,2))
        plot(thetas(1:i,5),thetas(1:i,6))
        
        set(cpos1,"XData",thetas(i,1),"YData",thetas(i,2))
        set(cpos2,"XData",thetas(i,5),"YData",thetas(i,6))
        
        pause(0.1)
    end
end


%% D

time_range = [0,2];

theta0 = [1; 0; 0; 3.158; 1.05; 0; 0; 7.64];
opts = odeset("RelTol", 1e-12);
[ts, thetas] = ode45(@diffeq_d, time_range, theta0, opts);


% plotting
if ismember("d", segments)
    clf; hold on
    axis([-2, 2, -2, 2])

    % plot sun
    plot(P(1), P(2),"Marker", "O")

    cpos1 = plot(thetas(:,1),thetas(:,2),"Marker", "O", "Color","red");
    cpos2 = plot(thetas(:,5),thetas(:,6),"Marker", "O", "Color","green");
    
    
    for i=1:25:length(ts(:,1))
        plot(thetas(1:i,1),thetas(1:i,2), "Color", "red")
        plot(thetas(1:i,5),thetas(1:i,6), "Color", "green")
        
        set(cpos1,"XData",thetas(i,1),"YData",thetas(i,2))
        set(cpos2,"XData",thetas(i,5),"YData",thetas(i,6))
        
        pause(0.1)
    end
end


%% E

time_range = [0,2];

theta0 = [1.05;0;0;7.64];


opts = odeset("RelTol", 1e-12);
[ts, thetas] = ode45(@diffeq_e, time_range, theta0, opts);



% plotting
if ismember("e", segments)
    clf; hold on
    axis([-2, 2, -2, 2])

    % plot sun
    plot(P(1), P(2),"Marker", "O")
    
    omega = 3.158;
    cpos1 = plot(cos(0),sin(0),"Marker", "O", "Color","red");
    cpos2 = plot(thetas(:,1),thetas(:,2),"Marker", "O", "Color","green");
    
    
    for i=1:25:length(ts(:,1))
        plot(cos(omega*ts(1:i)),sin(omega*ts(1:i)), "Color", "red")
        plot(thetas(1:i,1),thetas(1:i,2), "Color", "green")
        
        set(cpos1,"XData",cos(omega*ts(i)),"YData",sin(omega*ts(i)))
        set(cpos2,"XData",thetas(i,1),"YData",thetas(i,2))
        
        pause(0.1)
    end
end

%% f
time_range = [0,2];

theta0 = [1.05;0;0;7.64];


opts = odeset("RelTol", 1e-12);
[ts, thetas] = ode45(@diffeq_e, time_range, theta0, opts);


% plotting
if ismember("f", segments)
    clf; hold on
    axis([-0.5, 0.5, -0.5, 0.5])

    % plot planet in forced [0,0]
    plot(0, 0,"Marker", "O", "Color", "Red")
    
    omega = 3.158;
    
    cpos = plot(thetas(:,1)-cos(0),thetas(:,2)-sin(0),"Marker", "O", "Color","green");
    
    
    for i=1:25:length(ts(:,1))
        
        % pos of planet
        x = cos(omega*t(i));
        y = sin(omega*t(i));
        
        xs = cos(omega*t(1:i));
        ys = sin(omega*t(1:i)); 
        
  
        plot(thetas(1:i,1)-xs,thetas(1:i,2)-ys, "Color", "green")
        

        set(cpos,"XData",thetas(i,1)-x,"YData",thetas(i,2)-y)
        
        pause(0.1)
    end
end

%% g

time_range = [0,2];

theta0 = [1; 0; 0; 3.158; 1.05; 0; 0; 7.64; 0; 0; 0; -0.3158];
opts = odeset("RelTol", 1e-12);
[ts, thetas] = ode45(@diffeq_g, time_range, theta0, opts);




% plotting
if ismember("g", segments)
    clf; hold on
    axis([-2, 2, -2, 2])


    cpos1 = plot(thetas(:,1),thetas(:,2),"Marker", "O", "Color","red");
    cpos2 = plot(thetas(:,5),thetas(:,6),"Marker", "O", "Color","green");
    cpos3 = plot(thetas(:,9),thetas(:,10),"Marker", "O", "Color","blue");
    
    for i=1:25:length(ts(:,1))
        plot(thetas(1:i,1),thetas(1:i,2), "Color", "red")
        plot(thetas(1:i,5),thetas(1:i,6), "Color", "green")
        plot(thetas(1:i,9),thetas(1:i,10), "Color", "blue")
        
        set(cpos1,"XData",thetas(i,1),"YData",thetas(i,2))
        set(cpos2,"XData",thetas(i,5),"YData",thetas(i,6))
        set(cpos3,"XData",thetas(i,9),"YData",thetas(i,10))
        
        pause(0.1)
    end
end


%% h
time_range = [0,2];

theta0 = [1; 0; 0; 3.158; 1.05; 0; 0; 7.64; 0; 0; 0; -0.3158];
opts = odeset("RelTol", 1e-12);
[ts, thetas] = ode45(@diffeq_g, time_range, theta0, opts);




% plotting
if ismember("h", segments)
    clf; hold on
    axis([-2, 2, -2, 2])


    cpos1 = plot(thetas(:,1),thetas(:,2),"Marker", "O", "Color","red");
    cpos2 = plot(thetas(:,5),thetas(:,6),"Marker", "O", "Color","green");
    
    
    % force sun to [0,0]
    cpos3 = plot(0,0,"Marker", "O", "Color","blue");
    
    for i=1:25:length(ts(:,1))
        % real pos of sun
        x = thetas(i,9);
        y = thetas(i,10);
        
        xs = thetas(1:i,9);
        ys = thetas(1:i,10); 
        
        
        plot(thetas(1:i,1)-xs,thetas(1:i,2)-ys, "Color", "red")
        plot(thetas(1:i,5)-xs,thetas(1:i,6)-ys, "Color", "green")
        
        set(cpos1,"XData",thetas(i,1)-x,"YData",thetas(i,2)-y)
        set(cpos2,"XData",thetas(i,5)-x,"YData",thetas(i,6)-y)
        
        pause(0.1)
    end
end

%% i


time_range = [0,10];

theta0 = [1; -0.25; 0.41; 0.41; -1; -0.25; 0.41; 0.41; 0; 0; -0.82; -0.82];
opts = odeset("RelTol", 1e-12);
[ts, thetas] = ode45(@diffeq_i, time_range, theta0, opts);

% plotting
if ismember("i", segments)
    clf; hold on
    axis([-2, 2, -2, 2])


    cpos1 = plot(thetas(:,1),thetas(:,2),"Marker", "O", "Color","red");
    cpos2 = plot(thetas(:,5),thetas(:,6),"Marker", "O", "Color","green");
    cpos3 = plot(thetas(:,9),thetas(:,10),"Marker", "O", "Color","blue");
    
    for i=1:25:length(ts(:,1))
        plot(thetas(1:i,1),thetas(1:i,2), "Color", "red")
        plot(thetas(1:i,5),thetas(1:i,6), "Color", "green")
        plot(thetas(1:i,9),thetas(1:i,10), "Color", "blue")
        
        set(cpos1,"XData",thetas(i,1),"YData",thetas(i,2))
        set(cpos2,"XData",thetas(i,5),"YData",thetas(i,6))
        set(cpos3,"XData",thetas(i,9),"YData",thetas(i,10))
        
        pause(0.1)
    end
end


%% j


time_range = [0,10];


r = [0.9700436,-0.24308753];
v = [0.466203685,0.43236573];
theta0 = [r(1);r(2);v(1);v(2);-r(1);-r(2);v(1);v(2);0;0;-2*v(1);-2*v(2)];
opts = odeset("RelTol", 1e-12);
[ts, thetas] = ode45(@diffeq_i, time_range, theta0, opts);

% plotting
if ismember("j", segments)
    clf; hold on
    axis([-2, 2, -2, 2])


    cpos1 = plot(thetas(:,1),thetas(:,2),"Marker", "O", "Color","red");
    cpos2 = plot(thetas(:,5),thetas(:,6),"Marker", "O", "Color","green");
    cpos3 = plot(thetas(:,9),thetas(:,10),"Marker", "O", "Color","blue");
    
    for i=1:25:length(ts(:,1))
        plot(thetas(1:i,1),thetas(1:i,2), "Color", "red")
        plot(thetas(1:i,5),thetas(1:i,6), "Color", "green")
        plot(thetas(1:i,9),thetas(1:i,10), "Color", "blue")
        
        set(cpos1,"XData",thetas(i,1),"YData",thetas(i,2))
        set(cpos2,"XData",thetas(i,5),"YData",thetas(i,6))
        set(cpos3,"XData",thetas(i,9),"YData",thetas(i,10))
        
        pause(0.1)
    end
end

%% Diff eqs
function dthetadt = diffeq_a(t, theta)
       
    % position and velocity of planet
    p = [theta(1), theta(2)];
    v = [theta(3), theta(4)];
    
    % mass of planet
    m = 1;
    
    % mass of sun and pos of sun
    % is there a way such that I don't have to declare these twice
    % which is also compatible with ode45?
    P = [0,0];
    M = 10;
    
    % set grav constant to 1 because we're lazy
    G = 1;
    
    % from grav eq in asgnment
    a = (G*M)/norm(P-p)^3 *(P-p);
    
    dpdt = v;
    dvdt = a;
    
    dthetadt = [dpdt(1); dpdt(2); dvdt(1); dvdt(2)];

end

function [val, stop, dir] = event(t,theta)
    val = dot([theta(1), theta(2)],[theta(3), theta(4)]);
    stop = 1;
    dir = 0;
end

function dthetadt = diffeq_c(t,theta)
    % position and velocity of planets
    p1 = [theta(1), theta(2)];
    v1 = [theta(3), theta(4)];
    
    p2 = [theta(5), theta(6)];
    v2 = [theta(7), theta(8)];
    
    % mass of planets
    m = 1;
    
    % yes, still lazy
    G = 1;
    
    P = [0,0];
    M = 10;
    
    
    a1 = (G*M)/norm(P-p1)^3*(P-p1) + (G*m)/norm(p2-p1)^3*(p2-p1);
    a2 = (G*M)/norm(P-p2)^3*(P-p2) + (G*m)/norm(p1-p2)^3*(p1-p2);
    
    dp1dt = v1;
    dv1dt = a1;
    dp2dt = v2;
    dv2dt = a2;
    
    dthetadt = transpose([dp1dt,dv1dt,dp2dt,dv2dt]);
    

end

function dthetadt = diffeq_d(t,theta)
    % position and velocity of planets
    p1 = [theta(1), theta(2)];
    v1 = [theta(3), theta(4)];
    
    p2 = [theta(5), theta(6)];
    v2 = [theta(7), theta(8)];
    
    % mass of planets
    m = [1,0.001];
    
    % yes, still lazy
    G = 1;
    
    P = [0,0];
    M = 10;
    
    
    a1 = (G*M)/norm(P-p1)^3*(P-p1) + (G*m(2))/norm(p2-p1)^3*(p2-p1);
    a2 = (G*M)/norm(P-p2)^3*(P-p2) + (G*m(1))/norm(p1-p2)^3*(p1-p2);
    
    dp1dt = v1;
    dv1dt = a1;
    dp2dt = v2;
    dv2dt = a2;
    
    dthetadt = transpose([dp1dt,dv1dt,dp2dt,dv2dt]);
    

end

function dthetadt = diffeq_e(t,theta)
    % position and velocity of object
    p = [theta(1), theta(2)];
    v = [theta(3), theta(4)];
    
    % mass of planets
    m = 1;
    
    % yes, still lazy
    G = 1;
    
    P = [0,0];
    M = 10;
    
    % planet motion
    omega = 3.158;
    
    p_static = [cos(omega*t), sin(omega*t)];
    
    a = (G*M)/norm(P-p)^3*(P-p) + (G*m)/norm(p_static-p)^3*(p_static-p);
    
    dpdt = v;
    dvdt = a;
    
    dthetadt = transpose([dpdt,dvdt]);
    

end


function dthetadt = diffeq_g(t,theta)
    % position and velocity of objects
    
    % planet
    p1 = [theta(1), theta(2)];
    v1 = [theta(3), theta(4)];
    
    % moon
    p2 = [theta(5), theta(6)];
    v2 = [theta(7), theta(8)];
    
    % sun
    p3 = [theta(9), theta(10)];
    v3 = [theta(11), theta(12)];
    
    % mass of planets
    m = [1,0.001,10];
    
    % yes, still lazy
    G = 1;
    
    a1 = (G*m(2))/norm(p2-p1)^3*(p2-p1) + (G*m(3))/norm(p3-p1)^3*(p3-p1);
    a2 = (G*m(1))/norm(p1-p2)^3*(p1-p2) + (G*m(3))/norm(p3-p2)^3*(p3-p2);
    a3 = (G*m(1))/norm(p1-p3)^3*(p1-p3) + (G*m(2))/norm(p2-p3)^3*(p2-p3);
    
    dp1dt = v1;
    dv1dt = a1;
    dp2dt = v2;
    dv2dt = a2;
    dp3dt = v3;
    dv3dt = a3;
    
    dthetadt = transpose([dp1dt,dv1dt,dp2dt,dv2dt,dp3dt,dv3dt]);
    

end


function dthetadt = diffeq_i(t,theta)
    % position and velocity of objects
    
    % planet
    p1 = [theta(1), theta(2)];
    v1 = [theta(3), theta(4)];
    
    % moon
    p2 = [theta(5), theta(6)];
    v2 = [theta(7), theta(8)];
    
    % sun
    p3 = [theta(9), theta(10)];
    v3 = [theta(11), theta(12)];
    
    % mass of planets
    m = [1,1,1];
    
    % yes, still lazy
    G = 1;
    
    a1 = (G*m(2))/norm(p2-p1)^3*(p2-p1) + (G*m(3))/norm(p3-p1)^3*(p3-p1);
    a2 = (G*m(1))/norm(p1-p2)^3*(p1-p2) + (G*m(3))/norm(p3-p2)^3*(p3-p2);
    a3 = (G*m(1))/norm(p1-p3)^3*(p1-p3) + (G*m(2))/norm(p2-p3)^3*(p2-p3);
    
    dp1dt = v1;
    dv1dt = a1;
    dp2dt = v2;
    dv2dt = a2;
    dp3dt = v3;
    dv3dt = a3;
    
    dthetadt = transpose([dp1dt,dv1dt,dp2dt,dv2dt,dp3dt,dv3dt]);
    

end