设惯性系为基坐标系Ω1，某运动物体的本体坐标系为Ω2，该物体通过3次欧拉转动使原始本体坐标系变为新的坐标系Ω3，已知Ω2求坐标系Ω3，使用c语言实现

假设欧拉转动的顺序为 Z-Y-X，即先绕 Z 轴转动，再绕 Y 轴转动，最后绕 X 轴转动。

首先需要定义旋转矩阵，以 Z 轴旋转为例，旋转角度为θ，则旋转矩阵为：

Rz = [cos(θ)  -sin(θ)  0]
     [sin(θ)   cos(θ)  0]
     [0        0       1]

同理，Y 轴旋转矩阵为：

Ry = [cos(θ)   0   sin(θ)]
     [0        1   0     ]
     [-sin(θ)  0   cos(θ)]

X 轴旋转矩阵为：

Rx = [1   0        0       ]
     [0   cos(θ)  -sin(θ)]
     [0   sin(θ)   cos(θ)]

将三个旋转矩阵按照顺序相乘即可得到总的旋转矩阵 R：

R = Rx * Ry * Rz

假设初始时物体在 Ω1 坐标系下的位置向量为 P1，物体在 Ω2 坐标系下的位置向量为 P2，则 Ω3 坐标系下的位置向量 P3 可以用旋转矩阵 R 将 P2 转换为 Ω1 坐标系下的位置向量 P1，再用同样的方法将 P1 转换为 Ω3 坐标系下的位置向量 P3：

P1 = R * P2
P3 = inverse(R) * P1

其中 inverse(R) 表示旋转矩阵 R 的逆矩阵，可以使用矩阵求逆的方法得到。

下面是使用 C 语言实现的代码：

#include <stdio.h>
#include <math.h>

typedef struct {
    double x;
    double y;
    double z;
} Vector3;

typedef struct {
    double m[3][3];
} Matrix3x3;

// Z 轴旋转矩阵
Matrix3x3 Rz(double theta) {
    Matrix3x3 m = {
        {cos(theta), -sin(theta), 0},
        {sin(theta),  cos(theta), 0},
        {         0,           0, 1}
    };
    return m;
}

// Y 轴旋转矩阵
Matrix3x3 Ry(double theta) {
    Matrix3x3 m = {
        { cos(theta), 0, sin(theta)},
        {          0, 1,          0},
        {-sin(theta), 0, cos(theta)}
    };
    return m;
}

// X 轴旋转矩阵
Matrix3x3 Rx(double theta) {
    Matrix3x3 m = {
        {1,          0,           0},
        {0, cos(theta), -sin(theta)},
        {0, sin(theta),  cos(theta)}
    };
    return m;
}

// 求矩阵的逆矩阵
Matrix3x3 inverse(Matrix3x3 m) {
    Matrix3x3 inv;
    double det =
        m.m[0][0] * (m.m[1][1] * m.m[2][2] - m.m[1][2] * m.m[2][1]) -
        m.m[0][1] * (m.m[1][0] * m.m[2][2] - m.m[1][2] * m.m[2][0]) +
        m.m[0][2] * (m.m[1][0] * m.m[2][1] - m.m[1][1] * m.m[2][0]);
    double inv_det = 1.0 / det;
    inv.m[0][0] =  (m.m[1][1] * m.m[2][2] - m.m[1][2] * m.m[2][1]) * inv_det;
    inv.m[0][1] = -(m.m[0][1] * m.m[2][2] - m.m[0][2] * m.m[2][1]) * inv_det;
    inv.m[0][2] =  (m.m[0][1] * m.m[1][2] - m.m[0][2] * m.m[1][1]) * inv_det;
    inv.m[1][0] = -(m.m[1][0] * m.m[2][2] - m.m[1][2] * m.m[2][0]) * inv_det;
    inv.m[1][1] =  (m.m[0][0] * m.m[2][2] - m.m[0][2] * m.m[2][0]) * inv_det;
    inv.m[1][2] = -(m.m[0][0] * m.m[1][2] - m.m[0][2] * m.m[1][0]) * inv_det;
    inv.m[2][0] =  (m.m[1][0] * m.m[2][1] - m.m[1][1] * m.m[2][0]) * inv_det;
    inv.m[2][1] = -(m.m[0][0] * m.m[2][1] - m.m[0][1] * m.m[2][0]) * inv_det;
    inv.m[2][2] =  (m.m[0][0] * m.m[1][1] - m.m[0][1] * m.m[1][0]) * inv_det;
    return inv;
}

// 将向量 v 乘以矩阵 m
Vector3 mul(Matrix3x3 m, Vector3 v) {
    Vector3 r;
    r.x = m.m[0][0] * v.x + m.m[0][1] * v.y + m.m[0][2] * v.z;
    r.y = m.m[1][0] * v.x + m.m[1][1] * v.y + m.m[1][2] * v.z;
    r.z = m.m[2][0] * v.x + m.m[2][1] * v.y + m.m[2][2] * v.z;
    return r;
}

int main() {
    // 假设欧拉转动的顺序为 Z-Y-X
    double theta1 = M_PI / 4;   // 绕 Z 轴旋转的角度
    double theta2 = M_PI / 3;   // 绕 Y 轴旋转的角度
    double theta3 = M_PI / 6;   // 绕 X 轴旋转的角度

    // 物体在 Ω2 坐标系下的位置向量
    Vector3 P2 = {1, 2, 3};

    // 计算旋转矩阵 R
    Matrix3x3 R = Rx(theta3);
    R = mul(Ry(theta2), R);
    R = mul(Rz(theta1), R);

    // 计算 P1 和 P3
    Vector3 P1 = mul(R, P2);
    Vector3 P3 = mul(inverse(R), P1);

    printf("P1 = (%.2f, %.2f, %.2f)\n", P1.x, P1.y, P1.z);
    printf("P3 = (%.2f, %.2f, %.2f)\n", P3.x, P3.y, P3.z);

    return 0;
}