Create local coordinate system and update throughout flight

I am trying to monitor the position of a rocket throughout its trajectory by creating a local coordinate system. I would like it to update throughout the flight. My vector calculus is a bit rusty, so I'm having trouble getting the right angles to get the new axes. Right now, I'm trying to plot the new axes in a crude way, but I would appreciate any help in figuring this out or advice for a better way to do this.

close;
clear;
clc;

t = 0:pi/16:pi;
x = sin(t);
y = cos(t);
z = t;

[m,n] = size(t);

radius = zeros(1,n);
theta = zeros(1,n);   % angle from +z axis to vector (zenith)
phi = zeros(1,n);     % angle from +x axis to vector projection on xy-plane (azimuth)
rocket_angle = zeros(1,n);

lsm = 0.25; % line size modifier

for i=1:length(x)
  plot3(x(i),y(i),z(i),'*r');
  %plot3(x(i),y(i),0,'*b'); % vector projection
  hold on;
  pause(0.01);
  radius(i) = sqrt(x(i)^2 + y(i)^2 + z(i)^2);
  theta(i) = acosd(z(i)/radius(i)); % global coordinate system
  phi(i) = atand(y(i)/x(i)); % global coordinate system
  
  if i > 1
    rocket_angle(i) = atand((y(i)-y(i-1))/(x(i)-x(i-1)));
    % show rocket x-axis
    line([x(i)-lsm*sind(rocket_angle(i)) x(i)+lsm*sind(rocket_angle(i))], [y(i)+lsm*sind(rocket_angle(i)) y(i)-lsm*sind(rocket_angle(i))], [z(i)-lsm*sind(rocket_angle(i)) z(i)+lsm*sind(rocket_angle(i))]);
    % show rocket y-axis
    % show rocket z-axis (cross product of x and y)
  end
  
end

xlabel('X');
ylabel('Y');
zlabel('Z');
grid on;

Answer 1

[darova]

Just find tangent of a trajectory

x = 0:.1:pi;
y = sin(x);
u0 = diff(x);               % component of a tangent
v0 = diff(y);               % component of a tangent
u1 = u0./hypot(u0,v0);      % normalization
v1 = v0./hypot(u0,v0);
line(x,y)
h = line(0,0,'color', 'r');
axis equal
for i = 1:length(x)-1
    set(h,'xdata',x(i)+[u1(i) 0 -v1(i)]/2);     % tangent and normal
    set(h,'ydata',y(i)+[v1(i) 0  u1(i)]/2);     % tangent and normal
    pause(0.5)
end