How to plot a second coordinate frame using plot3?

So I am trying to fully understand the intermediate steps of what i did and to do so I would like to plot everything. What I am doing is I have a point rotated by 45 degrees in the z axis of the reference frame, and then rotated by 90 degrees around the new y-axis. The point after the two rotations is at [1; 0;1 ]. I have found what the point looks in the reference frame:

alb = deg2rad(45);

Rz = [ cos(alb) sin(alb) 0;
       -sin(alb) cos(alb) 0
       0 0 1];

beta = deg2rad(90);

Ry = [ cos(beta) 0 sin(beta);
     0 1 0;
     -sin(beta) 0 cos(beta)];

p_body = [1 0 1];

p_ref = Rz*Ry*p'

How can I plot the new coordinate axes after the 45 degrees rotation, and then after the 90 degree rotation? How shoudl I use the following?

Around X-axis:

X = x;
Y = y*cos(theta) - z*sin(theta);
Z = y*sin(theta) + z*cos(theta);

Around Y-axis:

X = x*cos(theta) + z*sin(theta);
Y = y;
Z = z*cos(theta) - x*sin(theta);

Around Z-axis:

X = x*cos(theta) - y*sin(theta);
Y = x*sin(theta) + y*cos(theta);
Z = z;

Answer 1

[Jeffrey Clark]

@Jimmy Neutron, your Z rotation matrix doesn't agree with your Z rotation formula. Either could be correct depending if you want clockwise or counter rotation. I'd go with only matrices in right-hand rule (see wiki page) and skip the fomulas.

To draw the axis after rotation choose points on each of the axis in the original frame, rotate them, plot their lines on a figure using plot3 (or quiver3 if you want arrows) and use the figure's rotation control or using view to view them from some perspective. So:

Rx = @(t) [      1       0       0  ...
          ;      0   cos(t) -sin(t) ...
          ;      0   sin(t)  cos(t) ];

Ry = @(t) [  cos(t)      0   sin(t) ...
          ;      0       1       0  ...
          ; -sin(t)      0   cos(t) ];

Rz = @(t) [  cos(t) -sin(t)      0  ...
          ;  sin(t)  cos(t)      0  ...
          ;      0       0       1  ];

Xold = [1 0 0]; % point [0 0 0] will not move in these rotations
Yold = [0 1 0]; % but would be needed if you also had translatons
Zold = [0 0 1]; % in your matrices

Rxy = Rx(deg2rad(0));
Ryb = Ry(deg2rad(90));
Rza = Rz(deg2rad(45));

Xnew = Rza*Ryb*Rxy*Xold';
Ynew = Rza*Ryb*Rxy*Yold';
Znew = Rza*Ryb*Rxy*Zold';

tiledlayout(2,1,'TileSpacing','compact','Padding','compact')

nexttile
hold on
displayAxis(Xnew,Ynew,Znew)
xlabel('Xnew')
ylabel('Ynew')
zlabel('Znew')
axis([-1.5 1.5 -1.5 1.5 -1.5 1.5])
axis square
axis equal
view(25,45)

nexttile
hold on
displayAxis(Xnew,Ynew,Znew)
xlabel('Xnew')
ylabel('Ynew')
zlabel('Znew')
axis([-1.5 1.5 -1.5 1.5 -1.5 1.5])
axis square
axis equal
camup(Znew)
campos(2*mean([Xnew';-Ynew';Znew']))

function displayAxis(Xnew,Ynew,Znew)
    quiver3(0,0,0,Xnew(1),Xnew(2),Xnew(3))
    quiver3(0,0,0,Ynew(1),Ynew(2),Ynew(3))
    quiver3(0,0,0,Znew(1),Znew(2),Znew(3))
    text(Xnew(1),Xnew(2),Xnew(3),'Xold')
    text(Ynew(1),Ynew(2),Ynew(3),'Yold')
    text(Znew(1),Znew(2),Znew(3),'Zold')
end

[George Abrahams]

Just to note, instead of your displayAxis function, you could use my plotframe function on File Exchange, if you want something more powerful. Hope it helps.