function [state1, state2] = m_fnShiXu(state)
% 时序函数,提供起飞、级间分离与点火的判据
% 确保state是结构体
if ~isstruct(state)
error('state必须是结构体');
end
% 确保const子结构体存在
if ~isfield(state, 'const') || ~isstruct(state.const)
state.const = struct();
end
% 确保dhconst字段存在
if ~isfield(state.const, 'dhconst')
state.const.dhconst = 0.01; % 设置默认值
end
dh = state.const.dhconst; % 现在安全访问
% 初始化阶段标志和文件指针
state.fpCombined = -1; % 初始化为无效文件标识符
state.fpCombined1 = -1; % 初始化为无效文件标识符
% 载入(气动)参数
state = m_fnLoadData();
state.const.Flag_write=1;
% 存储(弹道)参数
if state.const.Flag_write
% 设置输出数据的保存位置
outputDir = './OUTPUT/';
% 确保输出目录存在
if ~exist(outputDir, 'dir')
mkdir(outputDir);
end
% 尝试打开文件
state.fpCombined = fopen(fullfile(outputDir, 'CombinedData.csv'), 'w');
if state.fpCombined == -1
error('无法打开文件 CombinedData.csv');
end
state.fpCombined1 = fopen(fullfile(outputDir, 'LamdaBWeiHeight.csv'), 'w');
if state.fpCombined1 == -1
fclose(state.fpCombined); % 关闭已打开的文件
error('无法打开文件 LamdaBWeiHeight.csv');
end
% 写入CSV表头
fprintf(state.fpCombined, 'x0[1],x0[2],x0[3],x0[4],x0[5],x0[6],x0[7],phi*RTD,psi*RTD,gamma*RTD,theta*RTD,sigma*RTD,alpha*RTD,beta*RTD,Height,Vnorm,ma,MassRoc,thetaK*RTD,q,F_Engine[1],Lamda*RTD,BWei*RTD,Range,dW22[1],dW22[2],dW22[3],g[1],g[2],g[3],P*101.325,T,midu,p_mathpri.OrElement[1],p_mathpri.OrElement[2],p_mathpri.OrElement[3]*RTD,p_mathpri.OrElement[4]*RTD,p_mathpri.OrElement[1]*(1+p_mathpri.OrElement[2])-6378,p_mathpri.OrElement[1]*(1-p_mathpri.OrElement[2])-6378,nx,nz\n');
end
const= init_rocket_state();
% ====================================================================
% 一级发动机工作阶段
% ====================================================================
state.TimeLeft = const.TimeP1;
state.Flag_PowerOff = 0;
state.phase = 1;
state.iflag = 0;
step = 1;
while ~state.Flag_PowerOff
if state.iflag ~= 0
state.TimeLeft = state.TimeLeft -state.dh;
end
if (state.TimeLeft < state.dh) && (abs(state.TimeLeft - state.dh) > state.const.ErrorTime)
state.dh = state.TimeLeft;
end
% 执行龙格-库塔法积分
state.x0 = m_fnRK(state.x0, state.dh);
% 调用动态计算函数获取所有参数
[state, dx] = m_fnDynamic(state,const);
if state.Flag_write
% 存储结果
state.results(step, :) = [state.x0(1), state.x0(2), state.x0(3), state.x0(4), state.x0(5), state.x0(6), state.x0(7), ...
state.phi*state.const.RTD, state.psi*state.const.RTD, state.gamma*state.const.RTD, ...
state.theta*state.const.RTD, state.sigma*state.const.RTD, state.alpha*state.const.RTD, ...
state.beta*state.const.RTD, state.Height, state.Vnorm, state.ma, state.MassRoc, ...
state.thetaK*state.const.RTD, state.q, state.F_Engine(1), state.Lamda*state.const.RTD, ...
state.BWei*state.const.RTD, state.Range, state.dW22(1), state.dW22(2), state.dW22(3), ...
state.g(1), state.g(2), state.g(3), state.P*101.325, state.T, state.midu, ...
state.OrElement(1), state.OrElement(2), state.OrElement(3)*state.const.RTD, ...
state.OrElement(4)*state.const.RTD, state.OrElement(1)*(1+state.OrElement(2))-6378, ...
state.OrElement(1)*(1-state.OrElement(2))-6378, state.nx, state.nz];
% 写入LamdaBWeiHeight文件
fprintf(state.fpCombined1, '%10.12f,%10.12f,%10.12f\n', state.Height, state.Lamda*state.const.RTD, state.BWei*state.const.RTD);
end
% 检查是否到达关机时间
if abs(state.TimeLeft - state.dh) < state.const.ErrorTime
state.Flag_PowerOff = 1;
end
state.iflag = state.iflag + 1;
step = step + 1;
end
% 保存一级关机时刻的状态
state1 = [state.x0(1), state.x0(2), state.x0(3), state.x0(4), state.x0(5), state.x0(6), state.x0(7), state.phase];
% 输出关键参数到控制台
fprintf(' 时间(s) 高度(m) 绝对速度(m/s) 当地弹道倾角(°) 攻角(°) 俯仰角(°)\n');
fprintf(' 一级关机:%f %f %f %f %f %f\n', ...
state.x0(1), state.Height, state.VEInorm, state.thetaK*state.const.RTD, state.alpha*state.const.RTD, state.phi*state.const.RTD);
% ====================================================================
% 滑翔1阶段
% ====================================================================
state.phase = 2;
state.TimeLeft = state.const.TimeH1;
state.Flag_PowerOff = 0;
state.iflag = 0;
while ~state.Flag_PowerOff
if state.iflag ~= 0
state.TimeLeft = state.TimeLeft - state.dh;
end
if (state.TimeLeft < state.dh) && (abs(state.TimeLeft - state.dh) > state.const.ErrorTime)
state.dh = state.TimeLeft;
end
% 执行龙格-库塔法积分
state.x0 = m_fnRK(state.x0, state.dh);
% 调用动态计算函数获取所有参数
[state, dx] = m_fnDynamic(state);
if state.Flag_write
% 写入CombinedData文件
fprintf(state.fpCombined, '%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f\n', ...
state.x0(1), state.x0(2), state.x0(3), state.x0(4), state.x0(5), state.x0(6), state.x0(7), ...
state.phi*state.const.RTD, state.psi*state.const.RTD, state.gamma*state.const.RTD, ...
state.theta*state.const.RTD, state.sigma*state.const.RTD, state.alpha*state.const.RTD, ...
state.beta*state.const.RTD, state.Height, state.Vnorm, state.ma, state.MassRoc, ...
state.thetaK*state.const.RTD, state.q, state.F_Engine(1), state.Lamda*state.const.RTD, ...
state.BWei*state.const.RTD, state.Range, state.dW22(1), state.dW22(2), state.dW22(3), ...
state.g(1), state.g(2), state.g(3), state.P*101.325, state.T, state.midu, ...
state.OrElement(1), state.OrElement(2), state.OrElement(3)*state.const.RTD, ...
state.OrElement(4)*state.const.RTD, state.OrElement(1)*(1+state.OrElement(2))-6378, ...
state.OrElement(1)*(1-state.OrElement(2))-6378, state.nx, state.nz);
% 写入LamdaBWeiHeight文件
fprintf(state.fpCombined1, '%10.12f,%10.12f,%10.12f\n', state.Height, state.Lamda*state.const.RTD, state.BWei*state.const.RTD);
end
% 检查是否到达关机时间
if abs(state.TimeLeft - state.dh) < state.const.ErrorTime
state.Flag_PowerOff = 1;
end
state.iflag = state.iflag + 1;
end
fprintf('滑翔第一次:%f %f %f %f %f %f\n', ...
state.x0(1), state.Height, state.VEInorm, state.thetaK*state.const.RTD, state.alpha*state.const.RTD, state.phi*state.const.RTD);
% ====================================================================
% 二级第一次开机
% ====================================================================
state.phase = 3;
state.TimeLeft = state.const.TimeP21;
state.Flag_PowerOff = 0;
state.iflag = 0;
% 初始化状态保存变量
state2 = zeros(8, 1);
state.bFirstCapture = true; % 用于记录是否已捕获110km高度状态
while ~state.Flag_PowerOff
if state.iflag ~= 0
state.TimeLeft = state.TimeLeft - state.dh;
end
if (state.TimeLeft < state.dh) && (abs(state.TimeLeft - state.dh) > state.const.ErrorTime)
state.dh = state.TimeLeft;
end
% 执行龙格-库塔法积分
state.x0 = m_fnRK(state.x0, state.dh);
% 调用动态计算函数获取所有参数
[state, dx] = m_fnDynamic(state);
if state.Flag_write
% 写入CombinedData文件
fprintf(state.fpCombined, '%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f\n', ...
state.x0(1), state.x0(2), state.x0(3), state.x0(4), state.x0(5), state.x0(6), state.x0(7), ...
state.phi*state.const.RTD, state.psi*state.const.RTD, state.gamma*state.const.RTD, ...
state.theta*state.const.RTD, state.sigma*state.const.RTD, state.alpha*state.const.RTD, ...
state.beta*state.const.RTD, state.Height, state.Vnorm, state.ma, state.MassRoc, ...
state.thetaK*state.const.RTD, state.q, state.F_Engine(1), state.Lamda*state.const.RTD, ...
state.BWei*state.const.RTD, state.Range, state.dW22(1), state.dW22(2), state.dW22(3), ...
state.g(1), state.g(2), state.g(3), state.P*101.325, state.T, state.midu, ...
state.OrElement(1), state.OrElement(2), state.OrElement(3)*state.const.RTD, ...
state.OrElement(4)*state.const.RTD, state.OrElement(1)*(1+state.OrElement(2))-6378, ...
state.OrElement(1)*(1-state.OrElement(2))-6378, state.nx, state.nz);
% 写入LamdaBWeiHeight文件
fprintf(state.fpCombined1, '%10.12f,%10.12f,%10.12f\n', state.Height, state.Lamda*state.const.RTD, state.BWei*state.const.RTD);
end
% 检查是否到达关机时间
if abs(state.TimeLeft - state.dh) < state.const.ErrorTime
state.Flag_PowerOff = 1;
end
% 捕获110km高度状态
if (state.Height / 1000 >= 110) && state.bFirstCapture
state2(1:7) = state.x0;
state2(8) = state.phase;
state.bFirstCapture = false;
end
state.iflag = state.iflag + 1;
end
fprintf('二级一次关机:%f %f %f %f %f %f\n', ...
state.x0(1), state.Height, state.VEInorm, state.thetaK*state.const.RTD, state.alpha*state.const.RTD, state.phi*state.const.RTD);
% ====================================================================
% 滑翔2阶段
% ====================================================================
state.phase = 4;
state.TimeLeft = state.const.TimeH2;
state.Flag_PowerOff = 0;
state.iflag = 0;
while ~state.Flag_PowerOff
if state.iflag ~= 0
state.TimeLeft = state.TimeLeft - state.dh;
end
if (state.TimeLeft < state.dh) && (abs(state.TimeLeft - state.dh) > state.const.ErrorTime)
state.dh = state.TimeLeft;
end
% 执行龙格-库塔法积分
state.x0 = m_fnRK(state.x0, state.dh);
% 调用动态计算函数获取所有参数
[state, dx] = m_fnDynamic(state);
if state.Flag_write
% 写入CombinedData文件
fprintf(state.fpCombined, '%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f\n', ...
state.x0(1), state.x0(2), state.x0(3), state.x0(4), state.x0(5), state.x0(6), state.x0(7), ...
state.phi*state.const.RTD, state.psi*state.const.RTD, state.gamma*state.const.RTD, ...
state.theta*state.const.RTD, state.sigma*state.const.RTD, state.alpha*state.const.RTD, ...
state.beta*state.const.RTD, state.Height, state.Vnorm, state.ma, state.MassRoc, ...
state.thetaK*state.const.RTD, state.q, state.F_Engine(1), state.Lamda*state.const.RTD, ...
state.BWei*state.const.RTD, state.Range, state.dW22(1), state.dW22(2), state.dW22(3), ...
state.g(1), state.g(2), state.g(3), state.P*101.325, state.T, state.midu, ...
state.OrElement(1), state.OrElement(2), state.OrElement(3)*state.const.RTD, ...
state.OrElement(4)*state.const.RTD, state.OrElement(1)*(1+state.OrElement(2))-6378, ...
state.OrElement(1)*(1-state.OrElement(2))-6378, state.nx, state.nz);
% 写入LamdaBWeiHeight文件
fprintf(state.fpCombined1, '%10.12f,%10.12f,%10.12f\n', state.Height, state.Lamda*state.const.RTD, state.BWei*state.const.RTD);
end
% 检查是否到达关机时间
if abs(state.TimeLeft - state.dh) < state.const.ErrorTime
state.Flag_PowerOff = 1;
end
state.iflag = state.iflag + 1;
end
fprintf('滑翔第二次:%f %f %f %f %f %f\n', ...
state.x0(1), state.Height, state.VEInorm, state.thetaK*state.const.RTD, state.alpha*state.const.RTD, state.phi*state.const.RTD);
% ====================================================================
% 二级第二次开机
% ====================================================================
state.phase = 5;
state.TimeLeft = state.const.TimeP22;
state.Flag_PowerOff = 0;
state.iflag = 0;
while ~state.Flag_PowerOff
if state.iflag ~= 0
state.TimeLeft = state.TimeLeft - state.dh;
end
if (state.TimeLeft < state.dh) && (abs(state.TimeLeft - state.dh) > state.const.ErrorTime)
state.dh = state.TimeLeft;
end
% 执行龙格-库塔法积分
state.x0 = m_fnRK(state.x0, state.dh);
% 调用动态计算函数获取所有参数
[state, dx] = m_fnDynamic(state);
if state.Flag_write
% 写入CombinedData文件
fprintf(state.fpCombined, '%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f,%10.12f\n', ...
state.x0(1), state.x0(2), state.x0(3), state.x0(4), state.x0(5), state.x0(6), state.x0(7), ...
state.phi*state.const.RTD, state.psi*state.const.RTD, state.gamma*state.const.RTD, ...
state.theta*state.const.RTD, state.sigma*state.const.RTD, state.alpha*state.const.RTD, ...
state.beta*state.const.RTD, state.Height, state.Vnorm, state.ma, state.MassRoc, ...
state.thetaK*state.const.RTD, state.q, state.F_Engine(1), state.Lamda*state.const.RTD, ...
state.BWei*state.const.RTD, state.Range, state.dW22(1), state.dW22(2), state.dW22(3), ...
state.g(1), state.g(2), state.g(3), state.P*101.325, state.T, state.midu, ...
state.OrElement(1), state.OrElement(2), state.OrElement(3)*state.const.RTD, ...
state.OrElement(4)*state.const.RTD, state.OrElement(1)*(1+state.OrElement(2))-6378, ...
state.OrElement(1)*(1-state.OrElement(2))-6378, state.nx, state.nz);
% 写入LamdaBWeiHeight文件
fprintf(state.fpCombined1, '%10.12f,%10.12f,%10.12f\n', state.Height, state.Lamda*state.const.RTD, state.BWei*state.const.RTD);
end
% 检查是否到达关机时间
if abs(state.TimeLeft - state.dh) < state.const.ErrorTime
state.Flag_PowerOff = 1;
end
state.iflag = state.iflag + 1;
end
fprintf('二级二次关机:%f %f %f %f %f %f\n', ...
state.x0(1), state.Height, state.VEInorm, state.thetaK*state.const.RTD, state.alpha*state.const.RTD, state.phi*state.const.RTD);
% 输出轨道参数
fprintf('半长轴:%f km\n', state.OrElement(1));
fprintf('偏心率:%f\n', state.OrElement(2));
fprintf('轨道倾角:%f °\n', state.OrElement(3)*state.const.RTD);
fprintf('升交点赤经:%f °\n', state.OrElement(4)*state.const.RTD);
fprintf('远地点高度:%f km\n', state.OrElement(1)*(1+state.OrElement(2))-6378);
fprintf('近地点高度:%f km\n', state.OrElement(1)*(1-state.OrElement(2))-6378);
% 关闭文件
if state.Flag_write && state.fpCombined ~= -1
fclose(state.fpCombined);
end
if state.Flag_write && state.fpCombined1 ~= -1
fclose(state.fpCombined1);
end
end
最新发布