DCM与欧拉角相互转换

“声明”：本程序欧拉角→DCM所得结果与MATLAB自带的转换函数有出入，即为其转置。只是因为不同学科的惯用方法不同。使用者可以修改代码，灵活使用。

一、定义旋转矩阵

分别定义X，Y，Z旋转矩阵的表示方法：

1.1 X的旋转矩阵

function R_x = ROTX(phi_deg)
    phi = deg2rad(phi_deg);
    R_x = [1, 0, 0;
           0, cos(phi), sin(phi);
           0, -sin(phi), cos(phi)];
end

1.2 Y的旋转矩阵

function R_y = ROTY(theta_deg)
    theta = deg2rad(theta_deg);
    R_y = [cos(theta), 0, -sin(theta);
           0, 1, 0;
           sin(theta), 0, cos(theta)];  
end

1.3 Z的旋转矩阵

function R_z = ROTZ(psi_deg)
    psi = deg2rad(psi_deg);
    R_z = [cos(psi), sin(psi), 0;
           -sin(psi), cos(psi), 0;
           0, 0, 1];
end

二、编写调用函数

此处主要是DCM和欧拉角互相转换的函数，包含12种欧拉角的顺序转换。

2.1欧拉角→DCM

function R = eulerToDCM(phi_deg, theta_deg, psi_deg)
    typeDCM = input('请输入DCM类型：', 's');
    switch typeDCM
       case 'XZX'
           R=ROTX(psi_deg)*ROTZ(theta_deg)*ROTX(phi_deg);
       case 'XYX'
           R=ROTX(psi_deg)*ROTY(theta_deg)*ROTX(phi_deg);
       case 'YXY'
           R=ROTY(psi_deg)*ROTX(theta_deg)*ROTY(phi_deg);
       case 'YZY'
           R=ROTY(psi_deg)*ROTZ(theta_deg)*ROTY(phi_deg);
       case 'ZYZ'
           R=ROTZ(psi_deg)*ROTY(theta_deg)*ROTZ(phi_deg);
       case 'ZXZ'
           R=ROTZ(psi_deg)*ROTX(theta_deg)*ROTZ(phi_deg);
       case 'XZY'
           R=ROTY(psi_deg)*ROTZ(theta_deg)*ROTX(phi_deg);
       case 'XYZ'
           R=ROTZ(psi_deg)*ROTY(theta_deg)*ROTX(psi_deg);
       case 'YXZ'
           R=ROTZ(psi_deg)*ROTX(theta_deg)*ROTY(phi_deg);
       case 'YZX'
           R=ROTX(psi_deg)*ROTZ(theta_deg)*ROTY(phi_deg);
       case 'ZYX'
           R=ROTX(psi_deg)*ROTY(theta_deg)*ROTZ(phi_deg);
       case 'ZXY'
           R=ROTX(psi_deg)*ROTY(theta_deg)*ROTZ(pi_deg);
    end
    disp(R)
end

2.2欧拉角→DCM

function DCMrToEuler(DCM)
    typeDCM = input('请输入DCM类型：', 's');    
    switch typeDCM
        case 'XZX'
           phi=atan2(DCM(1,3),DCM(1,2));
           theta=acos(DCM(1,1));
           psi=atan2(DCM(3,1),-DCM(2,1));
           disp(['XZX的欧拉角是：']);
       case 'XYX'
           phi=atan2(DCM(1,2),-DCM(1,3));
           theta=acos(DCM(1,1));
           psi=atan2(DCM(2,1),DCM(3,1));
           disp(['XYX的欧拉角是：']);
       case 'YXY'
           phi=atan2(DCM(2,1),DCM(2,3));
           theta=acos(DCM(2,2));
           psi=atan2(DCM(1,2),-DCM(3,2));
           disp(['YXY的欧拉角是：']);
       case 'YZY'
           phi=atan2(DCM(2,3),-DCM(2,1));
           theta=acos(DCM(2,2));
           psi=atan2(DCM(3,2),DCM(1,2));
           disp(['YZY的欧拉角是：']);
       case 'ZYZ'
           phi=atan2(DCM(3,2),DCM(3,1));
           theta=atan2(sqrt(1-DCM(3,3)^2),DCM(3,3));
           psi=atan2(DCM(2,3),-DCM(1,3));
           disp(['ZYZ的欧拉角是：']);
       case 'ZXZ'
           phi=atan2(DCM(3,1),-DCM(3,2));
           theta=acos(DCM(3,3));
           psi=atan2(DCM(1,3),DCM(2,3));
           disp(['ZXZ的欧拉角是：']);
       case 'XZY'
           phi=atan2(DCM(2,3),DCM(2,2));
           theta=asin(-DCM(2,1));
           psi=atan2(DCM(3,1),DCM(1,1));
           disp(['XZY的欧拉角是：']);
       case 'XYZ'
           phi=atan2(-DCM(3,2),DCM(3,3));
           theta=asin(DCM(3,1));
           psi=atan2(-DCM(2,1),DCM(1,1));
           disp(['XYZ的欧拉角是：']);
       case 'YXZ'
           phi=atan2(DCM(3,1),DCM(3,3));
           theta=acos(-DCM(3,2));
           psi=atan2(DCM(1,2),DCM(2,2));
           disp(['YXZ的欧拉角是：']);
       case 'YZX'
           phi=atan2(-DCM(1,3),DCM(1,1));
           theta=asin(DCM(1,2));
           psi=atan2(-DCM(3,2),DCM(2,2));
           disp(['YZX的欧拉角是：']);
       case 'ZYX'
           phi=atan2(DCM(1,2),DCM(1,1));
           theta=asin(-DCM(1,3));
           psi=atan2(DCM(2,3),DCM(3,3));
           disp(['ZYX的欧拉角是：']);
       case 'ZXY'
           phi=atan2(-DCM(2,1),DCM(2,2));
           theta=acos(DCM(2,3));
           psi=atan2(-DCM(1,3),DCM(3,3));
           disp(['ZXY的欧拉角是：']);
    end
    if phi>pi/2
        phi = phi-pi;
    elseif  phi<-pi/2
        phi = phi + pi;
    end
    if psi>pi/2
        psi = psi-pi;
    elseif  psi<-pi/2
        psi = psi + pi;
    end
    output=[rad2deg(phi),rad2deg(theta),rad2deg(psi)];
    disp(output)
end

三、主应用程序

通过调用“二”中的两个函数进行转换，可以使用MATLAB制作APP。注意，输入欧拉角与DCM时均需要使用数组的方式输入。

type = input('工作类型(eulerToDCM按1，DCMToeuler按2)：');
switch type
    case 1
        A = input('请输入欧拉角：');
        R = eulerToDCM(A(1,1), A(1,2), A(1,3));
    case 2
        DCM = input('请输入DCM：');
        DCMrToEuler(DCM)
end

四、结语

通过输入数学公式，实现12种欧拉角与DCM相互转换。使用者可以使用“三”中的程序进行转换，也可以新建程序对“二”中的函数进行调用，或者摘取所需要的欧拉角顺序的代码进行编程。