MATLAB Code of Analysis

minradius=1.6; %cam min radius
maxradius=3.9; %cam max radius
omega=-1;
winglength=7.6;
bodytowing=1.2; %grounded joint of wing to slot
theta=(0:1:360);
thetastop=90-asind(minradius/maxradius); %point where period of no motion starts
radius=zeros(1,361);
for i=1:1:round(thetastop)+1
radius(i)=maxradius*cosd(i-1); %cosine of max radius, may drop below minradius as thetastop is rounded
end
for i=181:1:361
radius(i)=minradius+(maxradius-minradius)*(1-cosd(i-181))/2;%Simple harmonic motion
end
height=radius-minradius; %may drop below 0 as thetastop is rounded
for i=1:1:361 %if below 0, set to 0
if height(i)<0
height(i)=0;
end
end
velocity=zeros(1,361);
for i=1:1:round(thetastop)+1
velocity(i)=-maxradius*sind(i-1); %derivative of above
end
for i=181:1:361
velocity(i)=(maxradius-minradius)*sind(i-181)/2; %derivative of above
end
for i=1:1:361 %see explanation in position part
if height(i)==0
velocity(i)=0;
end
end
acceleration=zeros(1,361);
for i=1:1:round(thetastop)+1
acceleration(i)=-maxradius*cosd(i-1); %derivative of above
end
for i=181:1:361
acceleration(i)=(maxradius-minradius)*cosd(i-181)/2; %derivative of above
end
for i=1:1:361 %see explanation in position part
if height(i)==0
velocity(i)=0;
end
end
wingpositiony=-winglength*sind(45)+height*2*winglength*sind(45)/(maxradius-minradius); %body position curve scaled to end positions
wingpositionx=(winglength^2-wingpositiony.^2).^(1/2); %Pythagorean
wingvelocityy=winglength/bodytowing*velocity; %body velocity scaled by length
wingvelocityx=zeros(1,361);
for i=1:1:361
angle=atand(wingpositiony(i)/wingpositionx(i))+90;
wingvelocityx(i)=cotd(angle)*wingvelocityy(i); %Pythagorean
end
wingvelocitymag=(wingvelocityx.^2+wingvelocityy.^2).^(1/2);
wingangvelocity=wingvelocitymag/winglength;
for i=1:1:round(thetastop)+1
wingangvelocity(i)=-wingangvelocity(i); %clockwise during drop
end
wingaccelerationy=winglength/bodytowing*acceleration;
wingaccelerationx=zeros(1,361);
for i=2:1:360
wingaccelerationx(i)=(wingvelocityx(i+1)-wingvelocityx(i-1))/2*180/pi; %numerical differentiation
end
wingaccelerationx(1)=wingaccelerationx(2)^2/wingaccelerationx(3); %smoothing
wingaccelerationx(361)=wingaccelerationx(360)^2/wingaccelerationx(359); %smoothing
wingaccelerationx(181)=wingaccelerationx(182)^2/wingaccelerationx(183); %smoothing
wingaccelerationx(66)=wingaccelerationx(65)^2/wingaccelerationx(64); %smoothing
wingaccelerationx(67)=wingaccelerationx(66)*wingaccelerationx(65)/wingaccelerationx(64); %smoothing
wingangaccel=zeros(1,361);
for i=2:1:360
wingangaccel(i)=(wingangvelocity(i+1)-wingangvelocity(i-1))/2*180/pi; %numerical differentiation
end
wingangaccel(1)=wingangaccel(2)^2/wingangaccel(3); %smoothing
wingangaccel(361)=wingangaccel(360)^2/wingangaccel(359); %smoothing
wingangaccel(181)=wingangaccel(182)^2/wingangaccel(183); %smoothing
wingangaccel(66)=wingangaccel(65)^2/wingangaccel(64); %smoothing
wingangaccel(67)=wingangaccel(66)*wingangaccel(65)/wingangaccel(64); %smoothing
figure1=figure;
plot(theta,height);
title('Vertical Position of Penguin vs Theta')
xlabel('Theta (deg)');
ylabel('Vertical Position (cm)');
figure2=figure;
plot(theta,velocity);
title('Velocity of Penguin vs Theta')
xlabel('Theta (deg)');
ylabel('Velocity (cm/s)');
figure3=figure;
plot(theta,acceleration);
title('Acceleration of Penguin vs Theta')
xlabel('Theta (deg)');
ylabel('Acceleration (cm/s^2)');
figure4=figure;
plot(theta,wingpositionx,theta,wingpositiony);
title('Position of the Wing vs Theta')
xlabel('Theta (deg)');
ylabel('Wing Position (cm)');
legend('X','Y');
figure5=figure;
plot(theta,wingvelocityx,theta,wingvelocityy);
title('Velocity of the Wing vs Theta')
xlabel('Theta (deg)');
ylabel('Wing Velocity (cm/s)');
legend('X','Y');
figure6=figure;
plot(theta,wingaccelerationy);
title('Acceleration of the Wing in Y vs Theta')
xlabel('Theta (deg)');
ylabel('Wing Acceleration in Y (cm/s^2)');
figure7=figure;
plot(theta,wingangvelocity);
title('Angular Velocity of Wing vs Theta')
xlabel('Theta (deg)');
ylabel('Angular Velocity of wing (rad/s)');
figure8=figure;
plot(theta,wingaccelerationx);
title('Acceleration of the Wing in X vs Theta')
xlabel('Theta (deg)');
ylabel('Wing Acceleration in X (cm/s^2)');
figure9=figure;
plot(theta,wingangaccel);
title('Angular Acceleration of Wing vs Theta')
xlabel('Theta (deg)');
ylabel('Angular Acceleration of wing (rad/s^2)');
saveas(figure1,'BodyPosition.jpg');
saveas(figure2,'BodyVelocity.jpg');
saveas(figure3,'BodyAcceleration.jpg');
saveas(figure4,'WingPosition.jpg');
saveas(figure5,'WingVelocity.jpg');
saveas(figure6,'WingAccelerationY.jpg');
saveas(figure7,'WingAngularVel.jpg');
saveas(figure8,'WingAccelerationX.jpg');
saveas(figure9,'WingAngAccel.jpg');