Lab 3: Optimizing mechanisms
%%
close all; % closes all figures
clear all; % clears all variables from memory
clc;    % clears all calculations from the Matlab workspace

% Plot Parameters: these will be used to set the axis limits on the figure
xmin= -40; % leftmost window edge
xmax= 75; % rightmost window edge
ymin= -25; % bottom window edge
ymax= 80; % top window edge

%% Link Parameters
increments = 40; % number of theta2 configuration steps to calculate along mechanism rotation %%%% YOU MAY WANT TO CHANGE THIS

max_rotation_theta2 = 360 *pi/180; % rotation limit of theta2, radians
theta2 = linspace(0,max_rotation_theta2,increments); % link 2 rotation into 'increments' number of angles

r1 = 19;       % link 1 length, cm
r2 = 10;       % link 2 length, cm
r3 = 25;       % link 3 length, cm
r4 = 25;       % link 4 length, cm

a = 1.5 ; % path straightness tolerance distance
%% Rotation time
theta2_dot = 60;                      % rotation rate, rpm
t_rev = 60/theta2_dot;                % rotation time limit, seconds
time = linspace(0,t_rev,increments); % link 2 rotation into 'increments' number of times

%% set up figure
figure(1);                      %create new figure
% set(1,'WindowStyle','Docked')    %dock the figure

%% 4-bar mechanism geometric constants: See textbook section 4.3.3 for derivation. (Eq. 4.3-54)
h1 = r1/r2;
h2 = r1/r3;
h3 = r1/r4;
h4 =(-r1^2-r2^2-r3^2+r4^2)/(2*r2*r3);
h5 =(r1^2+r2^2-r3^2+r4^2)/(2*r2*r4);

%% calculate mechanism motion and plot it
num_extensions = 100;
extensionVector = linspace(5,50,num_extensions); % link 3 extension length

for j=1:length(extensionVector)          % cycle through each extension value
   for i=1:increments                      % step through motion of the mechanism
   hold off;

   % geometric calculations (book eq. 4.3-56 to 4.3-62):
   d = -h1 +(1-h3)*cos(theta2(i)) +h5;
   b = -2*sin(theta2(i));
   e = h1 -(1+h3)*cos(theta2(i)) +h5;
   a_a = -h1 +(1+h2)*cos(theta2(i)) +h4;
   c = h1 -(1-h2)*cos(theta2(i)) +h4;

   theta3_1(i) = 2*atan(((-b-(b^2-4*a_a*c)^0.5)/(2*a_a))); %calculate angle of link 3 (eq. 4.3-64)
   theta4_1(i) = 2*atan(((-b-(b^2-4*d*e)^0.5)/(2*d)));

   % Link Coordinates calculations:
   Ax(i) = 0;                                  % pivot point of link 2 position
   Ay(i) = 0;                                  % pivot point of link 2 position
   Bx(i) = r2*cos( theta2(i) );                % point B position
   By(i) = r2*sin( theta2(i) );                % point B position
   Cx(i,j) = Bx(i) + (r3+extensionVector(j))*cos( theta3_1(i) );          % point C position. We're going to store these for each extension value (indexed by j)
   Cy(i,j) = By(i) + (r3+extensionVector(j))*sin( theta3_1(i) );          % point C position
   Dx(i) = r1 + r4*cos( theta4_1(i) );          % point D position
   Dy(i) = r4*sin( theta4_1(i

4.1 The spreading of a toxic waste that decays in time is often mo...