From 0629fbf4f2ad9c379fe7c59f5e47df34eabe0575 Mon Sep 17 00:00:00 2001
From: =?utf8?q?Nils=20Forss=C3=A9n?= <forssennils@gmail.com>
Date: Wed, 7 Dec 2022 01:37:30 +0100
Subject: [PATCH] pushing

---
 brick_ekv.m    |   6 +++
 brick_events.m |   6 +++
 mek_upg.asv    |  90 +++++++++++++++++++++++++++++++++++
 mek_upg.m      | 125 ++++++++++++++++++++++++++++++++-----------------
 4 files changed, 185 insertions(+), 42 deletions(-)
 create mode 100644 brick_ekv.m
 create mode 100644 brick_events.m
 create mode 100644 mek_upg.asv

diff --git a/brick_ekv.m b/brick_ekv.m
new file mode 100644
index 0000000..713a49e
--- /dev/null
+++ b/brick_ekv.m
@@ -0,0 +1,6 @@
+function y_dot=brick_ekv(t, y, g)
+z = y(1);
+z_dot = y(2);
+y_dot = zeros(2,1);
+y_dot(1) = z_dot;
+y_dot(2) = -g;
\ No newline at end of file
diff --git a/brick_events.m b/brick_events.m
new file mode 100644
index 0000000..2b61cb9
--- /dev/null
+++ b/brick_events.m
@@ -0,0 +1,6 @@
+function [value,isterminal,direction] = brick_events(t, y, stopheight)
+
+value = y(1) - stopheight;
+isterminal=1;
+direction=0;
+end
\ No newline at end of file
diff --git a/mek_upg.asv b/mek_upg.asv
new file mode 100644
index 0000000..17fec0d
--- /dev/null
+++ b/mek_upg.asv
@@ -0,0 +1,90 @@
+%% init
+clear
+clc
+
+set(gcf, 'Position', [50 50 1200 900]);
+
+stopheight=0;
+g = 9.82;
+k = 100;
+m = 0.5;
+b = 0.3;
+
+t_max = 10;
+z_0 = sqrt(3)*b;
+
+hold on;
+
+%% calcs 
+options_collar = odeset("Events", @(t,y) collar_events(t, y, stopheight), ...
+    'RelTol', 1e-3,'AbsTol', 1e-6, 'MaxStep', 0.01);
+
+[t_vek_collar,Y_collar]=ode45(@(t,y) collar_ekv(t,y,g,k,m,b), [0 t_max], [z_0 0], options_collar);
+
+options_brick = odeset("Events", @(t,y) brick_events(t, y, stopheight), ...
+    'RelTol', 1e-3,'AbsTol', 1e-6, 'MaxStep', 0.01);
+
+[t_vek_brick,Y_brick]=ode45(@(t,y) brick_ekv(t,y,g), [0 t_max], [z_0 0], options_brick);
+
+z_collar = Y_collar(:, 1);
+z_dot_collar = Y_collar(:,2);
+s_collar = z_0 - z_collar;
+
+z_brick = Y_brick(:,1);
+z_dot_brick = Y_brick(:,2);
+s_brick = z_0 - z_brick; 
+
+%% plot 1 - hastighet
+subplot(2,2,1);
+
+p1 = plot(s_collar, z_dot_collar);
+p2 = plot(s_brick, z_dot_brick);
+xlim([0 max(s_brick)]);
+xlabel("Fallsträcka (m)");
+ylabel("Hastighet (m/s)");
+legend([p1 p2], {"Krage", "Kloss"});
+title("Hastighet");
+
+%% plot 2 - acceleration
+zdotdot_collar = zeros(length(t_vek_collar), 1);
+for i=1:length(t_vek_collar)
+    zdotdot_collar(i) = -g - ((k/m)*z_collar(i)*(1-((3*b)/(2*sqrt((z_collar(i)^2) + (b.^2))))));
+end
+subplot(2,2,2)
+p3 = plot(s_collar, zdotdot_collar);
+p4 = plot(s_brick, -g.*ones(length(s_brick), 1));
+xlim([0 max(s_brick)]);
+xlabel("Fallsträcka (m)");
+ylabel("Acceleration (m/s^2)");
+legend([p3 p4], {"Krage", "Kloss"}, "Location", "southeast");
+title("Acceleration")
+
+%% plot 3 - normalkraft
+subplot(2,2,3)
+normal = b.*k.*((1-(3.*b))./(2.*sqrt((s_brick.^2) + (b.^2))));
+p1 = plot(s_brick, normal);
+xlim([0 max(s_brick)]);
+xlabel("Fallsträcka (m)");
+ylabel("Normalkraft (N)");
+title("Normalkraft")
+
+%% plot 4 - energi
+T = (m.*(z_dot_collar.^2)./2;
+V_g = m*g*z_collar;
+V_e = k.*(sqrt((z_collar.^2) + (b.^2)) - (3.*b./2))./2;
+E_tot = T + V_g + V_e;
+subplot(2,2,4);
+p5 = plot(s_collar, T);
+p6 = plot(s_collar, V_g);
+p7 = plot(s_collar, V_e);
+p8 = plot(s_collar, E_tot);
+xlim([0 max(s_collar)]);
+xlabel("Fallsträcka (m)");
+ylabel("Energi (J)");
+legend([p5 p6 p7 p8], {"Rörelseenergi", "Potentiell energi", "Fjäderenergi", "Total energi"});
+title("Energi");
+
+%% print
+min(z_dot_collar)
+min(z_dot_brick)
+
diff --git a/mek_upg.m b/mek_upg.m
index dec2aaf..86532ce 100644
--- a/mek_upg.m
+++ b/mek_upg.m
@@ -1,59 +1,100 @@
-%BANAN
+%% init
+clear
+clc
+clf
+
+set(gcf, 'Position', [50 50 1200 900]);
+
 stopheight=0;
 g = 9.82;
 k = 100;
 m = 0.5;
 b = 0.3;
 
-
 t_max = 10;
 z_0 = sqrt(3)*b;
+%% calcs 
+options_collar = odeset('Events', @(t,y) collar_events(t, y, stopheight), ...
+    'RelTol', 1e-3,'AbsTol', 1e-6, 'MaxStep', 0.01);
 
+[t_vek_collar,Y_collar]=ode45(@(t,y) collar_ekv(t,y,g,k,m,b), [0 t_max], [z_0 0], options_collar);
 
-options = odeset("Events", @(t,y) collar_events(t, y, stopheight), ...
+options_brick = odeset('Events', @(t,y) brick_events(t, y, stopheight), ...
     'RelTol', 1e-3,'AbsTol', 1e-6, 'MaxStep', 0.01);
 
-[t_vek,Y]=ode45(@(t,y) collar_ekv(t,y,g,k,m,b), [0 t_max], [z_0 0], options);
+[t_vek_brick,Y_brick]=ode45(@(t,y) brick_ekv(t,y,g), [0 t_max], [z_0 0], options_brick);
 
-z = Y(:, 1);
-z_dot = Y(:,2);
-s = z_0 - z;
+z_collar = Y_collar(:, 1);
+z_dot_collar = Y_collar(:,2);
+s_collar = z_0 - z_collar;
 
-subplot(1,3,1);
-p1 = plot(s, z_dot);
-xlabel("s");
-ylabel("zdot");
-title("Hastighet krage");
+z_brick = Y_brick(:,1);
+z_dot_brick = Y_brick(:,2);
+s_brick = z_0 - z_brick; 
 
-zdotdot = zeros(length(t_vek), 1);
-for i=1:length(t_vek)
-    zdotdot(i) = -g - ((k/m)*z(i)*(1-((3*b)/(2*sqrt((z(i)^2) + (b.^2))))));
-end
-subplot(1,3,2)
-p2 = plot(s, zdotdot);
-xlabel("s");
-ylabel("ydotdot");
-title("Acceleration krage")
-
-step = 0.01;
-y = zeros(t_max/step, 1);
-v = zeros(t_max/step, 1);
-for i=0:length(new_t)
-
-    v(i+1) = g*step*i;
-    y(i+1) = z_0 - (v(i+1)*step*i/2);
-    if z_0 - (g*step*i*step*i/2) < 0
-        v = v(1:i, 1);
-        y = y(1:i, 1);
-        v = flip(v,1);
-
-        break
-    end
+%% plot 1 - hastighet
+subplot(2,2,1);
+p1 = plot(s_collar, z_dot_collar);
+hold on;
+p2 = plot(s_brick, z_dot_brick);
+xlim([0 max(s_collar) + 0.05*max(s_collar)]);
+xlabel('Fallsträcka (m)');
+ylabel('Hastighet (m/s)');
+hold off;
+legend([p1 p2], {'Krage', 'Kloss'});
+title('Hastighet');
+
+%% plot 2 - acceleration
+zdotdot_collar = zeros(length(t_vek_collar), 1);
+for i=1:length(t_vek_collar)
+    zdotdot_collar(i) = -g - ((k/m)*z_collar(i)*(1-((3*b)/(2*sqrt((z_collar(i)^2) + (b.^2))))));
 end
-subplot(1,3,3)
-p3 = plot(y, v);
-xlabel("s");
-ylabel("v");
+subplot(2,2,2)
+p3 = plot(s_collar, zdotdot_collar);
+hold on;
+p4 = plot(s_brick, -g.*ones(length(s_brick), 1));
+xlim([0 max(s_collar) + 0.05*max(s_collar)]);
+xlabel('Fallsträcka (m)');
+ylabel('Acceleration (m/s^2)');
+hold off;
+legend([p3 p4], {'Krage', 'Kloss'}, 'Location', 'southeast');
+title('Acceleration')
+
+%% plot 3 - normalkraft
+
+subplot(2,2,3)
+normal = b.*k.*((1-(3.*b))./(2.*sqrt((s_brick.^2) + (b.^2))));
+p1 = plot(s_brick, normal);
+hold on;
+xlim([0 max(s_brick) + 0.05*max(s_brick)]);
+xlabel('Fallsträcka (m)');
+ylabel('Normalkraft (N)');
+title('Normalkraft')
+hold off;
+
+%% plot 4 - energi
+z_collar(end + 1) = z_collar(end);
+z_dot_collar(end + 1) = 0;
+s_collar(end + 1) = s_collar(end);
+
+T = abs((m.*(z_dot_collar.^2))./2);
+V_g = abs(m*g*z_collar);
+V_e = abs(k.*(sqrt((z_collar.^2) + (b.^2)) - (3.*b./2))./2);
+E_tot = T + V_g + V_e;
+subplot(2,2,4);
+p5 = plot(s_collar, T);
+hold on;
+p6 = plot(s_collar, V_g);
+p7 = plot(s_collar, V_e);
+p8 = plot(s_collar, E_tot);
+xlim([0 max(s_collar) + 0.05*max(s_collar)]);
+xlabel('Fallsträcka (m)');
+ylabel('Energi (J)');
+hold off;
+legend([p5 p6 p7 p8], {'Rörelseenergi', 'Potentiell energi', 'Fjäderenergi', 'Total energi'}, 'Location', 'north');
+title('Energi');
+
+%% print
+min(z_dot_collar)
+min(z_dot_brick)
 
-max(v)
-z_dot(end)
-- 
2.30.2