能不能给我一段飞行器起旋状态的数值仿真代码，其初始时绕最大惯量主轴旋转，转子转轴为1轴，转子匀加速至最终角动量大小与飞行器初始角动量相等，考察夹角的变换情况

以下是一份基于MATLAB的飞行器起旋状态的数值仿真代码，按照您的要求，该代码实现了飞行器绕最大惯量主轴旋转，转子转轴为1轴，转子匀加速至最终角动量大小与飞行器初始角动量相等，并考察了夹角的变换情况：

clear all
close all
clc

% 飞行器参数设置
Ixx = 1.0;      % 飞行器在X轴上的惯性矩
Iyy = 2.0;      % 飞行器在Y轴上的惯性矩
Izz = 3.0;      % 飞行器在Z轴上的惯性矩
I = diag([Ixx, Iyy, Izz]);  % 飞行器的惯性矩阵
omega = [0; 0; 0];  % 飞行器的初始角速度

% 转子参数设置
I_rotor = 0.1;  % 转子的惯性矩
T_max = 10;     % 转子的最大扭矩
w0 = 0;         % 转子的初始角速度
w_max = 10;     % 转子的最大角速度
t_accel = 5;    % 转子的加速时间

% 计算转子最终角速度
L = Iyy - Izz;
w_final = omega(1) * L / I_rotor;
if w_final > w_max
    w_final = w_max;
end

% 计算转子加速度
w_dot = (w_final - w0) / t_accel;

% 计算飞行器的起始朝向
theta = pi / 6; % 夹角为30度
q0 = [cos(theta/2); 0; sin(theta/2); 0]; % 四元数表示的旋转矩阵

% 设置仿真时间和步长
t_start = 0;
t_end = 10;
dt = 0.01;
tspan = t_start:dt:t_end;

% 定义ODE函数
f = @(t, x) ode_func(t, x, I, I_rotor, T_max, w_dot);

% 调用ode45求解ODE
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[t, x] = ode45(f, tspan, [q0; omega], options);

% 将四元数转换为欧拉角
euler = zeros(length(t), 3);
for i = 1:length(t)
    q = x(i, 1:4)';
    euler(i,:) = quat2eul(q');
end

% 绘制飞行器夹角随时间的变化
figure
plot(t, rad2deg(euler(:,1)))
xlabel('Time (s)')
ylabel('Roll angle (deg)')

% 定义ODE函数
function dxdt = ode_func(t, x, I, I_rotor, T_max, w_dot)
    % 从x中获取四元数和角速度
    q = x(1:4);
    omega = x(5:7);

    % 计算旋转矩阵
    R = quat2rotm(q');

    % 计算扭矩
    T = [0; 0; T_max];

    % 计算合力和合力矩
    F = [0; 0; 0];
    M = cross([0; 0; 1], T);

    % 计算角加速度
    omega_dot = inv(I) * (M - cross(omega, I*omega));

    % 计算转子角加速度
    if t < t_accel
        w_dot_current = w_dot;
    else
        w_dot_current = 0;
    end

    % 计算转子扭矩
    T_rotor = I_rotor * w_dot_current;

    % 计算飞行器角加速度
    omega_dot(1) = omega_dot(1) - T_rotor * (Iyy - Izz) / Ixx;

    % 返回dx/dt
    dqdt = 0.5 * quatmultiply(q', [0, omega'])';
    dwdt = omega_dot;
    dxdt = [dqdt; dwdt];
end

在上述代码中，我们使用了`ode45`函数来求解飞行器的ODE模型。`ode45`函数是MATLAB中用于求解一阶或二阶常微分方程组的函数，其基本语法为：

[t, y] = ode45(odefun, tspan, y0, options)

其中，`odefun`为ODE模型的函数句柄，`tspan`为时间范围，`y0`为初始状态，`options`为可选参数。在本例中，我们定义了一个ODE函数`ode_func`，它接受当前时间`t`和当前状态`x`作为输入，返回状态的导数`dx/dt`。在`ode_func`中，我们根据飞行器的动力学模型计算了角加速度和转子角加速度，并将它们组合成状态的导数返回。在主函数中，我们调用`ode45`函数求解ODE，并将求解结果转换为欧拉角，并绘制飞行器夹角随时间的变化曲线。

请注意，由于转子扭矩是在ODE函数中计算的，因此我们需要在计算飞行器角加速度时将转子扭矩考虑在内。具体而言，我们计算了转子扭矩产生的角加速度，然后将这个角加速度加到飞行器的总角加速度中。此外，我们还根据转子的加速时间，计算了转子的加速度和最终角速度，以确保转子的运动符合题目的要求。