\documentclass{ctexart}

\usepackage{listings}
\usepackage{xcolor}
\usepackage{textcomp,url} % 必须加上，否则报错
\usepackage[framed,numbered,autolinebreaks,useliterate]{mcode}

\begin{document}
	\begin{lstlisting}[language=matlab]
		height = 1400;          %天问一号动力减速段的下落距离,约1500m到100m高度动力减速
		t = 0;                  %初始开始动力减速的时间为0，随迭代增加
		T_var = 0.2;            %迭代的时间步长为0.2s
		c = 2940;     %比冲c，用质量描述推进剂的量，比冲以速度量纲表现，国际单位为米每秒（m/s）
		g=3.71;       %火星重力加速度约为3.71m/s
		iteration = 1000:100:7500;     %推力从1000N到7500N迭代
		M_var = [];       %质量变化
		V_var = [];       %速度变化
		t_time = [];      %时间
		X_var = [];       %水平方向位移
		for i = iteration %迭代计算
		%动力减速段为主要消耗燃料的阶段
		%减速初始的运行状态变量数据
		Ft = i;           %推力Ft=迭代变量i，每个推力对应一种状态
		alpha=(11.2/180)*pi;     %初速度与水平面的夹角
		Vx0=94.172;             %水平初速度，进入动力减速段的初速度约为96m/s，水平初速是96*cos(11.2°)
		Vy0=18.645;             %竖直初速度，进入动力减速段的初速度约为96m/s，水平初速是96*sin(11.2°)
		m0 = 1285;          %初始质量，随燃料消耗质量降低
		Ay0 = g-Ft*sin(alpha)/(m0-Ft/c*t);%竖直初加速度，在竖直方向受重力加速度和由发动机推力产生的反向加速度
		Ax0 = -Ft*cos(alpha)/(m0-Ft/c*t);%水平初加速度，在水平方向发动机推力产生的反向加速度
		sum = 0;              %迭代计数，每迭代一次+1
		X_x = Vx0*t+0.5*Ax0*t^2;          %X方向位移分量
		Y_y = Vy0*t+0.5*Ay0*t^2;          %Y方向位移分量
		%% 迭代计算各个方向的速度分量和位移分量
		while(Y_y<height)               %循环迭代，当Y_res>=h时，结束
		sum = sum+1;
		%计算速度
		Vx = Vx0+Ax0*T_var;             %速度v=vo+a*t，每迭代一次速度变化一次
		Vy = Vy0+Ay0*T_var;
		Vx0 = Vx;                       %一次迭代的末状态量等于下一次初状态量
		Vy0 = Vy;
		%计算位移
		X = Vx0*T_var+1/2*Ax0*T_var^2;      %位移，由牛顿定律得x=v*t+1/2*a*t^2
		Y = Vy0*T_var+1/2*Ay0*T_var^2;
		X_x = X_x + X;              %总位移
		Y_y = Y_y + Y;
		%减速时间
		Time = sum*T_var;             %时间等于迭代计数*时间步长
		%加速度
		SIN = Vy/sqrt(Vy^2+Vx^2);               %末速度与水平方向角度正弦值
		COS = Vx/sqrt(Vy^2+Vx^2);               %末速度与水平方向角度余弦值
		Ay = g-Ft*SIN/(m0-Ft/c*Time);           %加速度分量
		Ax = -Ft*COS/(m0-Ft/c*Time);
		Ax0 = Ax;                       %一次迭代的末状态量等于下一次初状态量
		Ay0 = Ay;
		end
		M = m0-Ft/c*Time;               %末质量等于初质量减消耗的燃料质量
		Time = sum*T_var;%运动时间
		V_res = sqrt(Vx^2+Vy^2);            %末速度
		degree = atan(Vy/Vx)*180/pi;        %末速度角度
		M_var = [M_var;Ft/c*Time];          %质量变化
		V_var = [V_var;V_res,Vx];               %速度变化
		t_time = [t_time;Time];             %记录运行时间
		X_var = [X_var;X_x];
		end
		%画出相应的图像
		%图一：燃料-推力图
		figure;
		plot(iteration,M_var,'LineWidth',2);
		title('燃料消耗随推力变化');
		xlabel('推力 Ft ( N )');
		ylabel('燃油消耗质量 △M ( kg )');
		%图二：水平位移-推力图
		figure;
		plot(iteration,X_var,'LineWidth',2);
		title('水平方向位移随推力变化');
		xlabel('推力 Ft( N )');
		ylabel('位移 X ( m )');
		%图三：速度-推力图
		figure;
		plot(iteration,V_var(:,2),'LineWidth',2);
		title('末速度水平分量随推力变化');
		xlabel('推力 Ft(N)');
		ylabel('速度(m/s)');
		%图四：时间-推力图
		figure;
		plot(iteration,t_time,'LineWidth',2);
		title('时间随推力变化');
		xlabel('推力 Ft(N)');
		ylabel('时间 ( s )');
		Answer = [iteration',M_var,V_var,t_time,X_var];     %迭代计算中的计算结果
		xlswrite('data.xls',Answer);                 %结果导出到表格内，方便查询和分析数据
		%在Excel表格文件data.xls中，每列数据分别为
		%1、推力大小（单位：N）
		%2、燃料消耗质量（单位：kg）
		%3、总末速度(单位：m/s）
		%4、水平方向末速度分量（单位：m/s）
		%5、时间（单位：s)
		%通过分析表格内容，可以得到燃料最省的推力大小，绘制相应的轨迹图像
	\end{lstlisting}
\end{document}