r = 1;
theta = 60;
phi = 30;
x1 = r*sind(theta)*cosd(phi);
y1 = r*sind(theta)*sind(phi);
z1 = r* cosd(theta);
[x1,y1,z1]';


z = [0 0 r]';
z_prime = rotz(phi,'deg')*roty(theta,'deg')*z;
% 求出旋转轴 y_prime
y_prime = cross(z,[x1,y1,0]);
%y_prime单位化
y_prime = y_prime/norm(y_prime);
% 绕着y_prime轴旋转-theta角度，将z_prime轴回正为z轴，使用Euler-Rodrigues表示方法
rotAngle = deg2rad(-theta);

b = tan(rotAngle/2)*y_prime;

aform =["ZYX","ZYZ","ZXY","ZXZ",     "YXZ","YXY","YZX", "YZY","XYZ","XYX", "XZY","XZX"]';
arr = [];
for i = 1:length(aform)
    d = aform(i);
    [y,p,r] =  rod2angle(b,d);
    arr = [arr y];
end
% 获得旋转对应的四元数
qt  = rod2quat(b);
% 获得旋转的欧拉角
el = quat2eul(qt,'ZYZ');
% 由欧拉角获得旋转矩阵dcm
R = quat2dcm(qt);
%  旋转当前的z_prime轴回到之前的z轴
R'*z_prime