Skip to content

Instantly share code, notes, and snippets.

@kasperfred kasperfred/ml98.m
Created Feb 27, 2019

Embed
What would you like to do?
%% 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.