已知初始发射参数，通过设计运载火箭各阶段的弹道程序俯仰角、偏航角等控制变量，将有效载荷准确送入预定轨道，并给出弹道曲线以及各状态变量随时间的变化曲线 据以上条件给出一段matlab程序

以下是一段MATLAB代码，用于模拟火箭发射过程，包括控制变量和状态变量的计算和可视化展示。注意，这只是一个简单的示例，实际上需要更多的参数和控制变量来实现准确的轨道设计和控制。

% 初始化参数
m0 = 10000; % 火箭质量
mf = 5000; % 燃料质量
T = 100000; % 引擎推力
g = 9.81; % 重力加速度
Cd = 0.1; % 空气阻力系数
A = 10; % 火箭横截面积
h0 = 0; % 初始高度
v0 = 0; % 初始速度
theta0 = pi/4; % 初始俯仰角
phi0 = 0; % 初始偏航角

% 计算参数
m = m0 + mf; % 火箭初始总质量
vex = T/mf; % 燃料比冲
dt = 0.1; % 时间步长
t = 0:dt:1000; % 时间数组

% 初始化状态变量
h = zeros(size(t)); % 高度
v = zeros(size(t)); % 速度
theta = zeros(size(t)); % 俯仰角
phi = zeros(size(t)); % 偏航角
a = zeros(size(t)); % 加速度
M = zeros(size(t)); % 火箭质量
D = zeros(size(t)); % 空气阻力

% 初始状态
h(1) = h0;
v(1) = v0;
theta(1) = theta0;
phi(1) = phi0;
a(1) = T/m - g - 0.5*Cd*A/m*v(1)^2;
M(1) = m0 + mf;
D(1) = Cd*A*0.5*rho(h(1))*v(1)^2;

% 迭代计算状态变量
for i = 2:length(t)
    % 计算控制变量
    if h(i-1) < 5000
        theta(i) = theta0;
        phi(i) = phi0;
    elseif h(i-1) < 10000
        theta(i) = theta0 - 0.1*(h(i-1)-5000)/5000;
        phi(i) = phi0 + 0.1*(h(i-1)-5000)/5000;
    else
        theta(i) = 0;
        phi(i) = pi/4;
    end
    
    % 计算状态变量
    M(i) = m0 + mf*exp(-vex/g*(t(i)-t(1)));
    D(i) = Cd*A*0.5*rho(h(i-1))*v(i-1)^2;
    a(i) = T/M(i) - g - D(i)/M(i);
    v(i) = v(i-1) + a(i)*dt;
    h(i) = h(i-1) + v(i)*dt;
    
    % 检查是否燃料用尽
    if M(i) < m0
        mf = 0;
        vex = 0;
        T = 0;
    end
end

% 绘制弹道曲线和状态变量随时间的变化曲线
figure
subplot(2,2,1)
plot(t,h)
xlabel('时间（秒）')
ylabel('高度（米）')
title('火箭高度随时间的变化')

subplot(2,2,2)
plot(t,v)
xlabel('时间（秒）')
ylabel('速度（米/秒）')
title('火箭速度随时间的变化')

subplot(2,2,3)
plot(t,theta)
xlabel('时间（秒）')
ylabel('俯仰角（弧度）')
title('火箭俯仰角随时间的变化')

subplot(2,2,4)
plot(t,phi)
xlabel('时间（秒）')
ylabel('偏航角（弧度）')
title('火箭偏航角随时间的变化')

需要注意的是，这段代码中的控制变量仅仅是一个示例，实际上需要更加精细的控制算法来实现准确的轨道设计和控制。