Matlab Code for Penguin Kinematic Analysis

%Derek Hornung
%"Ready to Fly" Kinematic Analysis

theta_crank = [0:2*pi/360:2*pi];
theta_crank_deg = theta_crank * 180/pi;
omega_crank = 2*pi; %1 revolution per second CCW
%Deriving the vertical path of the penguin's body based on cam shape
%Cubic regression terms
a = (2.150E-2).*theta_crank(1:270).^3;
b = (2.026E-1).*theta_crank(1:270).^2;
c = (1.061E-2).*theta_crank(1:270);
d = 1.5;
%Assembling y_cam
y_cam(1:270) = -a + b - c + d;
y_cam(271:337) = 3.7*cos(theta_crank(271:337)-3*pi/2);
y_cam(338:361) = 1.5;

y_body = y_cam - 1.5; %Top of base as reference plane for penguin

%Calculating velocity of the penguin's body by taking dy/dt
%Terms of time derivative of cubic regression line of cam radius wrt angle
dadt = 3*(2.150E-2)*omega_crank.*theta_crank(1:270).^2;
dbdt = 2*(2.026E-1)*omega_crank.*theta_crank(1:270);
dcdt = (1.061E-2)*omega_crank;
dddt = 0;
%Assembling v_body
v_body(1:270) = -dadt + dbdt - dcdt + dddt;
v_body(271:337) = -3.7*omega_crank*sin(theta_crank(271:337)-3*pi/2);
v_body(338:361) = 0;

%Calculating acceleration of the penguin's body by taking dv/dt
%Terms of second derivative of cubic regression of cam radius wrt angle
d2adt2 = 3*(2.150E-2)*omega_crank^2.*theta_crank(1:270);
d2bdt2 = 2*(2.026E-1)*omega_crank^2;
d2cdt2 = 0;
%Assembling a_body
a_body(1:270) = -d2adt2 + d2bdt2 - d2cdt2;
a_body(271:337) = -3.7*omega_crank^2*cos(theta_crank(271:337)-3*pi/2);
a_body(338:361) = 0;

%Plotting kinematic analysis results for the motion of the penguin's body
figure()
plot(theta_crank_deg,y_body)
title('Vertical Path of Body wrt Angle of Crank')
xlabel('Angle of crank (degrees)')
ylabel('Position (cm)')

figure()
plot(theta_crank_deg,v_body)
title('Velocity of Body wrt Angle of Crank')
xlabel('Angle of crank (degrees)')
ylabel('Velocity (cm/s)')

figure()
plot(theta_crank_deg,a_body)
title('Acceleration of Body wrt Angle of Crank')
xlabel('Angle of crank (degrees)')
ylabel('Acceleration (cm/s^2)')

%Calculating angular position of wings wrt angle of crank
h_wing = 8.4; %Height of glued base of wings from surface of glacier
h_slot = 7.2; %Height of shoulder slots from bottom of penguin body
y_slot = y_body + h_slot; %Height of slot from surface of glacier
x_slot = 0.9; %Horizontal distance from base of wing to slot
theta_wing = atan2d(x_slot, (h_wing - y_slot));
omega_wing = v_body ./ x_slot;
alpha_wing = a_body ./ x_slot;

figure()
plot(theta_crank_deg,theta_wing)
title('Angle of Wing wrt Angle of Crank')
xlabel('Angle of crank (degrees)')
ylabel('Angle of wing (degrees)')

figure()
plot(theta_crank_deg,omega_wing)
title('Angular Velocity of Wing wrt Angle of Crank')
xlabel('Angle of crank (degrees)')
ylabel('Angular velocity of wing (rad/s)')

figure()
plot(theta_crank_deg,alpha_wing)
title('Angular Acceleration of Wing wrt Angle of Crank')
xlabel('Angle of crank (degrees)')
ylabel('Angular acceleration of wing (rad/s^2)')

%Calculating the vertical position of the wingtips
L_wing = 1.4 + 6.4; %Length of wing from base to tip
L_hidden = sqrt(x_slot^2 + (h_wing - y_slot).^2); %Length of wing in body
L_expose = L_wing - L_hidden; %Length of wing outside body
y_tip = h_wing - L_wing * cosd(theta_wing);
v_tip = L_wing .* omega_wing;
a_tip = L_wing .* alpha_wing;

figure()
plot(theta_crank_deg,y_tip)
title('Vertical Path of Wingtip wrt Angle of Crank')
xlabel('Angle of crank (degrees)')
ylabel('Position (cm)')

figure()
plot(theta_crank_deg,v_tip)
title('Velocity of Wingtip wrt Angle of Crank')
xlabel('Angle of crank (degrees)')
ylabel('Velocity (cm/s)')

figure()
plot(theta_crank_deg,a_tip)
title('Acceleration of Wingtip wrt Angle of Crank')
xlabel('Angle of crank (degrees)')
ylabel('Acceleration (cm/s^2)')

%Actual shape of cam (different range of theta for plotting purposes)
r_cam(92:361) = -a + b - c + d;
r_cam(1:91) = 1.5; %Radius of actual cam over interval 0 < theta < 90 
%r_cam(1:10:361)

x_cam_model = r_cam .* cos(theta_crank);
y_cam_model = r_cam .* sin(theta_crank);

figure()
plot(x_cam_model,y_cam_model)
xlabel('x (cm)')
ylabel('y (cm)')
title('Cam Shape')