9 - ReadyToFly.m

clear all;
theta = 0:.1:360;
w = 2;
s = zeros(1,length(theta));
v = zeros(1,length(theta));
a = zeros(1,length(theta));
vr = zeros(1,length(theta));
ar = zeros(1,length(theta));
alpha = zeros(1,length(theta));
thetarad = theta*pi/180;
for i = 1:length(theta)
    if (theta(i) >= 0 && theta(i) < 50)
        s(i) = 0;
        v(i) = 0;
        vr(i) = 0;
        a(i) = 0;
        ar(i) = 0;
    else
        s(i) = -0.0006406*thetarad(i)^6+0.0103723*thetarad(i)^5-0.0666412*thetarad(i)^4+0.1992861*thetarad(i)^3-0.207613*thetarad(i)^2+0.0764435*thetarad(i);
        v(i) = -0.0038436*w*thetarad(i)^5+0.0518615*w*thetarad(i)^4-0.2665648*w*thetarad(i)^3+0.5978583*w*thetarad(i)^2-0.415226*w*thetarad(i)+0.0764435*w;
        vr(i) = -0.0038436*thetarad(i)^5+0.0518615*thetarad(i)^4-0.2665648*thetarad(i)^3+0.5978583*thetarad(i)^2-0.415226*thetarad(i)+0.0764435;
        a(i) = -0.019218*w^2*thetarad(i)^4+0.207446*w^2*thetarad(i)^3-0.7996944*w^2*thetarad(i)^2+1.1957166*w^2*thetarad(i)-0.415226*w^2;
        ar(i) = -0.019218*thetarad(i)^4+0.207446*thetarad(i)^3-0.7996944*thetarad(i)^2+1.1957166*thetarad(i)-0.415226;
    end
    alpha(i) = atand(-ar(i)/vr(i))-theta(i);
end

thetacom = 296.2:.1:360;
thetarad2 = thetacom*pi/180;
ss = zeros(1,length(thetacom));
scom = zeros(1,length(thetacom));
for k = 1:length(thetacom)
    ss(k) = -0.0006406*thetarad2(k)^6+0.0103723*thetarad2(k)^5-0.0666412*thetarad2(k)^4+0.1992861*thetarad2(k)^3-0.207613*thetarad2(k)^2+0.0764435*thetarad2(k);
    if (thetacom(k) > 296.1 && thetacom(k) <= 349.2)
        scom(k) = sqrt(1.415^2-(1.415*sind(thetacom(k)-296.1))^2)-0.555;
    else
        scom(k) = sqrt((1.1316/sind(thetacom(k)-296.1))^2-1.1316^2)-0.555;
    end
end

d = 0.355;
L = 3.105;
phio = 65;
l = d/sind(phio);
fo = -l*cosd(phio);
thetafin = 110:0.1:296;
thetafinrad = thetafin*pi/180;
x = zeros(1,length(thetafin));
y = zeros(1,length(thetafin));
sfin = zeros(1,length(thetafin));
tfin = zeros(1,length(thetafin));
dphi = zeros(1,length(thetafin));
d2phi = zeros(1,length(thetafin));
phi = zeros(1,length(thetafin));
findis = zeros(1,length(thetafin));
finvel = zeros(1,length(thetafin));
finacc = zeros(1,length(thetafin));
xv = zeros(1,length(thetafin));
yv = zeros(1,length(thetafin));
xa = zeros(1,length(thetafin));
ya = zeros(1,length(thetafin));
for j = 1:length(thetafin)
    sfin(j) = -0.0006406*thetafinrad(j)^6+0.0103723*thetafinrad(j)^5-0.0666412*thetafinrad(j)^4+0.1992861*thetafinrad(j)^3-0.207613*thetafinrad(j)^2+0.0764435*thetafinrad(j)-s(1101);
    tfin(j) = thetafinrad(j)/w;
    dphi(j) = 1.7728*tfin(j)^3-10.1982*tfin(j)^2+17.8268*tfin(j)-8.3598;
    d2phi(j) = 5.3184*tfin(j)^2-20.3964*tfin(j)+17.8268;
    phi(j) = atand((fo+sfin(j))/d);
    x(j) = L*cosd(phi(j));
    y(j) = L*sind(phi(j));
    findis(j) = sqrt(x(j)^2+y(j)^2);

    xv(j) = -L*dphi(j)*sind(phi(j));
    yv(j) = L*dphi(j)*cosd(phi(j));
    finvel(j) = sqrt(xv(j)^2+yv(j)^2);

    xa(j) = -L*(d2phi(j)*sind(phi(j))+(dphi(j))^2*cosd(phi(j)));
    ya(j) = L*(d2phi(j)*cosd(phi(j))-(dphi(j))^2*sind(phi(j)));
    finacc(j) = sqrt(xa(j)^2+ya(j)^2);
end

% filename = 'phi.xlsx';
% writematrix(tfin,filename,'Sheet',1,'Range','A1');
% writematrix(phi,filename,'Sheet',1,'Range','A2');

figure(1);
plot(theta,s,'LineWidth',2);
hold on;
title('Ready To Fly: Penguin Displacement');
xlabel('Angle (deg)');
ylabel('Displacement (in)');

theta2  = [0
10
20
30
40
50
60
70
80
90
100
110
120
130
140
150
160
170
180
190
200
210
220
230
240
250
260
270
280
290
296
300
310
320
330
340
350
360];

thetacom2 = [300
310
320
330
340
350
360];

s2 = [0
0
0
0
0
0
0.0156
0.0313
0.0469
0.0625
0.0938
0.125
0.1563
0.1875
0.25
0.2813
0.3438
0.4063
0.4375
0.5
0.5313
0.5938
0.625
0.6563
0.7188
0.75
0.8125
0.8438
0.85
0.8594
0.86
0.8438
0.8125
0.75
0.625
0.5
0.3125
0];

scom2 = [0.8438
0.8125
0.75
0.625
0.5
0.3125
0];

plot(theta2,s2,'-o');
legend('Trendline','Measurements','Location','northwest')
hold off;

xx = [0
0
0
0
0
0
0.05
0.08
0.12
0.15
0.15
0.15
0.15
0.15
0.15
0.18
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.17
0.15
0.15
0.11
0.07
0.03
0
-0.15
-0.36
-0.7
-0.94
-1.08
-1.13];

xxtheta = [0
0.174532925
0.34906585
0.523598776
0.698131701
0.872664626
1.047197551
1.221730476
1.396263402
1.570796327
1.745329252
1.919862177
2.094395102
2.268928028
2.443460953
2.617993878
2.792526803
2.967059728
3.141592654
3.316125579
3.490658504
3.665191429
3.839724354
4.01425728
4.188790205
4.36332313
4.537856055
4.71238898
4.886921906
5.061454831
5.166174586
5.235987756
5.410520681
5.585053606
5.759586532
5.934119457
6.108652382];

figure(2);
subplot(2,1,1);
plot(theta,v,'LineWidth',2);
hold on;
title('Ready To Fly: Penguin Velocity (in/s)');
xlabel('Angle (deg)');
ylabel('Velocity (in/s)');
hold off;
subplot(2,1,2);
plot(thetarad,vr,'LineWidth',2);
hold on;
plot(xxtheta,xx,'-o');
title('Ready To Fly: Penguin Velocity (in/rad)');
xlabel('Angle (rad)');
ylabel('Velocity (in/rad)');
legend('d/dt(Trendline)','Measurements','Location','southwest');
hold off;

figure(3);
subplot(2,1,1);
plot(theta,a,'LineWidth',2);
hold on;
title('Ready To Fly: Penguin Acceleration (in/s^2)');
xlabel('Angle (deg)');
ylabel('Acceleration (in/s^2)');
hold off;
subplot(2,1,2);
plot(thetarad,ar,'LineWidth',2);
hold on;
title('Ready To Fly: Penguin Acceleration (in/rad^2)');
xlabel('Angle (rad)');
ylabel('Acceleration (in/rad^2)');
hold off;

figure(4);
plot(thetacom2,scom2,'-o','LineWidth',2);
hold on;
plot(thetacom,ss,'LineWidth',2);
plot(thetacom,scom,'LineWidth',2);
xlabel('Angle (deg)');
ylabel('Displacement (in)');
title('Penguin Displacement Equation Comparison');
legend('Measurements','Trendline','Trigonometry Equations');
hold off;

figure(5);
plot(x,y);
hold on;
title('Ready To Fly: Penguin Fin Position');
xlabel('Horizontal Position (in)');
ylabel('Vertical Position (in)');
hold off;

figure(6);
subplot(3,1,1);
plot(tfin,findis);
hold on;
title('Ready To Fly: Penguin Fin Motion Profile');
xlabel('Time (s)');
ylabel('Position Magnitude (in)');
hold off;

subplot(3,1,2);
plot(tfin,finvel);
hold on;
xlabel('Time (s)');
ylabel('Velocity Magnitude (in/s)');
hold off;

subplot(3,1,3);
plot(tfin,finacc);
hold on;
xlabel('Time (s)');
ylabel('Acceleration Magnitude (in/s^2)');
hold off;

figure(7);
subplot(3,2,1);
plot(tfin,x);
hold on;
title('Ready To Fly: Penguin Fin Motion Profile [Real]');
xlabel('Time (s)');
ylabel('Position [Real] (in)');
hold off;

subplot(3,2,3);
plot(tfin,xv);
hold on;
xlabel('Time (s)');
ylabel('Velocity [Real] (in/s)');
hold off;

subplot(3,2,5);
plot(tfin,xa);
hold on;
xlabel('Time (s)');
ylabel('Acceleration [Real] (in/s^2)');
hold off;

subplot(3,2,2);
plot(tfin,y);
hold on;
title('Ready To Fly: Penguin Fin Motion Profile [Imaginary]');
xlabel('Time (s)');
ylabel('Position [Imaginary] (in)');
hold off;

subplot(3,2,4);
plot(tfin,yv);
hold on;
xlabel('Time (s)');
ylabel('Velocity [Imaginary] (in/s)');
hold off;

subplot(3,2,6);
plot(tfin,ya);
hold on;
xlabel('Time (s)');
ylabel('Acceleration [Imaginary] (in/s^2)');
hold off;