软件功能

卫星平面轨道入轨仿真。即根据入轨点的参数，对其后的运动轨迹进行计算。可以模拟出关机点参数与轨道形状的关系。第一、第二宇宙速度等。

火箭自由飞行段弹道仿真。根据主动段结束时的关机点参数，推算自由飞行段轨迹。可以观察出射程、飞行时间、再入速度与关机点参数的关系。

近地面的抛体仿真。当关机点参数很小时，模型表现为水平地面的斜抛运动。

软件进行数值计算之后，会以动画效果演示飞行轨迹。

软件带有绘图保持功能，可以比较不同的初始参数的仿真结果。

软件运行效果

软件开始界面。

左下角的开始仿真按钮被水印盖住了。

软件帮助及说明。

弹道导弹自由飞行段的仿真。

通过“绘图保持”按钮，可以比较不同关机点参数下的弹道。每次改完参数，然后点击“绘图保持”，再点击“开始仿真”。

在图形的下方有个文本框，列出本次弹道参数，包括飞行时间s，射程km，弹头再入大气层的速度mk/s。

当弹道与地面没有交点时，得到卫星轨道。

这里要注意设置仿真时间，不能太长，不然动画会一直持续很久。

可以比较不同的速度、弹道倾角对应的轨道情况。速度达到一定程度，弹头逃逸，离开地球。

近地面抛体运动。

需要点击“曲线缩放”按钮，放大画面。这样才能看到这么小的抛体轨迹。当然也可以仿真，比较多个参数的效果。

动画演示。

仿真软件附件：

ballistic missile.mlappinstall 42.60KB MLAPPINSTALL 420次下载

运行次软件，需要提前下载Matlab可以安装APP的版本。我用的是Matlab2018b版本。下载之后双击自动安装，然后就可以在我的APP里找到这个应用。下图第二行的第一个。

软件使用的说明

动力学模型：

目前的模型还比较简单，但是基本的估算是没问题的，尤其是弹头自由飞行段和卫星较高的轨道。以后再考虑加上地球引力摄动模型和大气阻力模型，这样的话计算结果会更精确些，尤其是低轨道和弹道的再入段部分。但是相应的，需要输入的参数也会变多。

这是一个中心引力场下的二体运动模型，这个是有解析解的，解析解的部分可以参考我以前写的文章。当然也可以通过数值计算求解。

因为平面看图比较直观，而且圆锥曲线的轨道是个平面。于是投影到平面下进行计算仿真。

数值计算方法：

采用变步长的四阶龙格库塔法进行求解。具体的表达式如下

在Matlab里实现很方便，只需要就一句代码 :

Other
开启自动换行复制
 [T_out,Y_out]=ode45(@daodan,[t0 T],Y0,options)

Matlab独立仿真代码源码：

贴到Matlab里就能运行，无需安装。想改的话自己改吧~

Other
开启自动换行复制
%平面弹道导弹被动段弹道仿真（中心引力体运动）
clc;clear;close all;

T=2e4;%最长仿真时间s
higth_k = 10e3;%关机点高度m
speed_k = 7.9e3;%关机点速度大小m/s
theta_k = 50*pi/180;%关机点弹道倾角rad
latitude_k = -50*pi/180;%关机点纬度deg

t0=0;%仿真开始时间s
RE=6371e3;%地球半径m
rk=RE+higth_k;
Y0=[speed_k*sin(theta_k+latitude_k),-speed_k*cos(theta_k+latitude_k),rk*cos(latitude_k),rk*sin(latitude_k)]';
options=odeset('Reltol',1e-4,'AbsTol',1e-4,'Events',@eventfun);
[T_out,Y_out]=ode45(@daodan,[t0 T],Y0,options);

figure(1);
%绘图
h1=plot(Y_out(1,3)*1e-3,Y_out(1,4)*1e-3,'-.ob');
a1=gca;
grid on;
xlabel('x/km');
ylabel('z/km');
%绘制地球
q=0:0.01:2*pi;
x_earth = RE*cos(q);
z_earth = RE*sin(q);
hold on;
h2=plot([0,Y0(3)]*1e-3,[0,Y0(4)]*1e-3,'-ok');
h3=plot(x_earth*1e-3,z_earth*1e-3,'k','LineWidth',2);
axis equal;
for i=1:1:length(T_out)
    set(h1,'xdata',Y_out(1:i,3)*1e-3,'ydata',Y_out(1:i,4)*1e-3);
    pause(0.1);
end

latitude_end = atan2(Y_out(end,4),Y_out(end,3));
distance_fly = abs(latitude_end - latitude_k)*RE*1e-3;
fprintf('弹头飞行距离为：%.3f km\n',distance_fly);

%事件检测函数，数值积分终止条件
function [value,isterminal,direction] = eventfun(T_out,Y_out)
RE=6371e3;%地球半径m
value=sqrt(Y_out(3)*Y_out(3)+Y_out(4)*Y_out(4))-RE+T_out*0; %触发时间，当其值为0的时候，时间会触发 到地球表面的距离为0
isterminal=1;   %设为1时会，触发时间会停止求解器，设0时触发不影响工作
direction=-1;    %触发方向设1时是上升触发，设-1是下降触发，设0是双向触发
end

%二体运动方程
function dy=daodan(t,Y)
dy=zeros(4,1);%Y=[Vx,Vz,x,z]'m
GM=3.98603e14;%地球引力常数
R=sqrt(Y(3)*Y(3)+Y(4)*Y(4))+t*0;%到地心距离m
dy([1,2])=-GM/R^3*[Y(3),Y(4)]';%动力学方程m/s^2
dy([3,4])=Y([1,2]);%运动学方程m/s
end