5. MATLAB Code

% This script studies the kinematics of a dinosaur automata
clc;clear;
% links
a = 0.5; % crank
[b,c,e,f,g,h,i,j,k,l,m,n] = deal(1); % neck links
d = 1.332; % ground link

%% Angles (deg)
offset = -227; % change based on ground
thetaA_initial = 90;
thetaA = thetaA_initial+offset:1:thetaA_initial+offset+360;

[thetaD,thetaA,thetaB,thetaC] = fourbarpos(a,b,c,d,offset,thetaA,1);
[thetaG,thetaK,thetaE,thetaI,thetaM] = deal(thetaC);
[thetaF,thetaJ,thetaN,thetaL,thetaH] = deal(thetaB);

thetaA = thetaA-offset;
thetaB = thetaB-offset;
thetaC = thetaC-offset;
thetaD = thetaD-offset;
thetaE = thetaE-offset;
thetaF = thetaF-offset;
thetaG = thetaG-offset;
thetaH = thetaH-offset;
thetaI = thetaI-offset;
thetaJ = thetaJ-offset;
thetaK = thetaK-offset;
thetaL = thetaL-offset;
thetaM = thetaM-offset;
thetaN = thetaN-offset;

%% Position
pointP2x = a.*cosd(thetaA)+b.*cosd(thetaB)+e.*cosd(thetaE)+f.*cosd(thetaF)+i.*cosd(thetaI)+j.*cosd(thetaJ)+m.*cosd(thetaM);
pointP2y = a.*sind(thetaA)+b.*sind(thetaB)+e.*sind(thetaE)+f.*sind(thetaF)+i.*sind(thetaI)+j.*sind(thetaJ)+m.*sind(thetaM);

pointP1x = a.*cosd(thetaA)+b.*cosd(thetaB)+e.*cosd(thetaE)+f.*cosd(thetaF)+i.*cosd(thetaI)+j.*cosd(thetaJ)+n.*cosd(thetaN);
pointP1y = a.*sind(thetaA)+b.*sind(thetaB)+e.*sind(thetaE)+f.*sind(thetaF)+i.*sind(thetaI)+j.*sind(thetaJ)+n.*sind(thetaN);

pointHx = a.*cosd(thetaA)+b.*cosd(thetaB)+e.*cosd(thetaE)+f.*cosd(thetaF)+i.*cosd(thetaI)+j.*cosd(thetaJ);
pointHy = a.*sind(thetaA)+b.*sind(thetaB)+e.*sind(thetaE)+f.*sind(thetaF)+i.*sind(thetaI)+j.*sind(thetaJ);

pointP1x = real(pointP1x);
pointP1y = real(pointP1y);
pointP2x = real(pointP2x);
pointP2y = real(pointP2y);
pointHx = real(pointHx);
pointHy = real(pointHy);

% Position Plots 
figure(1)
hold on
plot(pointP1x, pointP1y);
plot(pointP2x, pointP2y);
plot(pointHx, pointHy);
axis('equal');
title('Movement of Dino Head');
xlabel('X Position of Dino Head');
ylabel('Y Position of Dino Head');
legend('P1', 'P2', 'Point H');
hold off

%% Angular Velocity
omegaA = 10; %rad/sec
omegaP1 = a/m.*sind(thetaC-thetaA)./sind(thetaB-thetaC);
omegaP2 = a/n.*sind(thetaA-thetaB)./sind(thetaC-thetaB);

figure(2)
hold on
plot(thetaA, omegaP1);
plot(thetaA, omegaP2);
title('Angular Velocity of Dino Head');
xlabel('Angle of Crank (deg)');
ylabel('Angular Velocity of Dino Head (rad/sec)');
legend('P1', 'P2');
hold off

%% Force Relationships
Fcrank = 1; % assume force applied to crank = 1 lb
Fp1 = omegaA.*a./omegaP1./m.*Fcrank;
Fp2 = omegaA.*a./omegaP2./n.*Fcrank;

TotalBiteForce = Fp2-Fp1;

figure(3)
subplot(3,1,1)
hold on
plot(thetaA, Fp2);
set(gca, 'XLim', [90 450],'YLim',[-50 50]);
title('Force of Top Jaw vs Crank Angle');
ylabel('Force (lb)');
legend('Top Jaw');
subplot(3,1,2)
plot(thetaA, Fp1);
set(gca, 'XLim', [90 450],'YLim',[-50 50]);
title('Force of Bottom Jaw vs Crank Angle');
ylabel('Force (lb)');
legend('Bottom Jaw');
subplot(3,1,3)
plot(thetaA,TotalBiteForce);
set(gca, 'XLim', [90 450],'YLim',[-50 50]);
title('Force of Dino Bite vs Crank Angle');
xlabel('Angle of Crank (deg)');
ylabel('Force (lb)');
legend('Total Bite Force');
hold off

mechAdvTopJaw = abs(Fp2/Fcrank);
mechAdvBotJaw = abs(Fp1/Fcrank);

figure(5)
hold on
plot(thetaA, mechAdvTopJaw);
plot(thetaA, mechAdvBotJaw);
set(gca, 'XLim', [90 450],'YLim',[0 50]);
title('Mechanical Advantage of Dino Bite vs Crank Angle');
xlabel('Angle of Crank (deg)');
ylabel('Mech Adv');
legend('Top Jaw', 'Bottom Jaw');
hold off

%% Animation of Dino Head Position

for i = i:length(thetaA)
 
 figure(4) 
 P1_circle = viscircles([pointP1x(i) pointP1y(i)], 0.05);
 P2_circle = viscircles([pointP2x(i) pointP2y(i)], 0.05);
 H_circle = viscircles([pointHx(i) pointHy(i)], 0.05);
 
 P1_bar = line([pointP1x(i) pointHx(i)], [pointP1y(i) pointHy(i)]);
 P2_bar = line([pointP2x(i) pointHx(i)], [pointP2y(i) pointHy(i)]);
 
 axis('equal');
 set(gca, 'XLim', [-7 0],'YLim',[0 7]);
 title('Dino Head Movement');
 
 labP1 = text(pointP1x(i), pointP1y(i)+0.2,'P1');
 labP2 = text(pointP2x(i), pointP2y(i)+0.2,'P2');
 labH = text(pointHx(i), pointHy(i)+0.2,'H');
 F(i)=getframe(gcf);
 pause(0.01);
 if i<length(thetaA)
 delete(P1_circle);
 delete(P2_circle);
 delete(H_circle);
 delete(P1_bar);
 delete(P2_bar);
 delete(labP1);
 delete(labP2);
 delete(labH);
 end
end
obj = VideoWriter('animation.mp4', 'MPEG-4');
open(obj);
writeVideo(obj,F);
close(obj)