clear
clc
syms t
syms x(t) y(t)

% x, y in meters; t in seconds
x(t) = -1000 + 60*t - (t-100)^3/1000;
y(t) = 20*(.005*t^2-t+50);

% 5.1
fplot(x,y,[0,200]);
title('Trajectory');
xlabel('x (m)');
ylabel('y (m)');

% 5.2
vx = diff(x,t); % m/s
vy = diff(y,t); % m/s

speed = sqrt(vx^2+vy^2); % m/s

fplot(speed,[0,200]);
title('Speed');
xlabel('Time (s)');
ylabel('Speed (m/s)');

% 5.3
% (i) Find accelerations in i,j basis
ax = diff(vx,t);
ay = diff(vy,t);

% (ii) Find a unit vector tangent to path
tn1 = vx/speed;
tn2 = vy/speed;

% (iii) Find a unit vector normal to path
% in 2D, a vector 90 deg ccw from a vector [A,B] is [-B,A]
n1 = -vy/speed;
n2 = vx/speed;

% (iv) Find an and at using dot products
an = ax*n1 + ay*n2;
at = ax*tn1 + ay*tn2;

fplot(at, [0,200]);
hold on;
fplot(an, [0,200]);
title('Acceleration');
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');
legend('Tangential', 'Normal')
hold off
% 5.4
g = 9.8; % m/s^2
% Definition of load factor
lf = (ax*n1 + ay * n2)/g + n2;

fplot(lf, [0,200])
title('Load Factor')
xlabel('Time (s)')
ylabel('Load Factor')