无控纵向导弹弹道计算

function myjob（）
clcclear
tic
indata
z=15;
n=8;
h=0.01;
Y=[0 20 z 0 z 0 0 0];  %初始值
fid=fopen（‘result.txt‘‘w‘）;
%将数据写入文本文件
for i=1:n
    fprintf（fid‘%8.4f ‘Y（i））;
end
old4=YFZ（Y）;
Y=RK（nhY）;
%write
for i=1:n
    fprintf（fid‘%8.4f ‘Y（i））;
end
old3=YFZ（Y）;
Y=RK（nhY）;
%write
for i=1:n
    fprintf（fid‘%8.4f ‘Y（i））;
end
old2=YFZ（Y）;
Y=RK（nhY）;
%write
for i=1:n
    fprintf（fid‘%8.4f ‘Y（i））;
end
old1=YFZ（Y）;
while Y（7）>=0
    [Yold1old2old3old4]=ADM（Yhold1old2old3old4）;
    %write
    for i=1:n
        fprintf（fid‘%8.4f ‘Y（i））;
    end
end
fclose（fid）;
Y
Y（1）
%结果处理
load result.txt
s=size（result2）;
da=reshape（result8s/8）;
da=da‘;
plot（da（:6）da（:7））title（‘弹道曲线‘）xlabel（‘x/m‘）ylabel（‘y/m‘）
figure
plot（da（:1）da（:2））title（‘速度曲线‘）xlabel（‘t/s‘）ylabel（‘V/（m/s）‘）
a=da（:5）-da（:3）;
figure
plot（da（:1）a）title（‘攻角曲线‘）xlabel（‘t/s‘）ylabel（‘a/rad‘）
toc

%阿当姆斯预估校正法
function [Yold1old2old3old4]=ADM（Yhold1old2old3old4）
y=Y+h*（55*old1-59*old2+37*old3-9*old4）/24;
f=YFZ（y）;
Y=Y+h*（9*f+19*old1-5*old2+old3）/24;
old4=old3;
old3=old2;
old2=old1;
old1=YFZ（Y）;

%四阶龙格库塔法
function Y=RK（nhY）
b=[h/2 h/2 h h];
dy=YFZ（Y）;
oldy=Y;
Y1=[];
for j=1:3
    for i=1:n
        Y1（i）=oldy（i）+b（j）*dy（i）;
        Y（i）=Y（i）+b（j+1）*dy（i）/3;
    end
    dy=YFZ（Y1）;
end
Y=Y+b（1）*dy/3;

%右端项函数
function dy=YFZ（Y）
ts=3.6;
tf=ts+0.5;
T=273.5+20;
Ma=Y（2）/（20*sqrt（T））;
m0=819.06;
S=0.308;
D=0.36;
RHO=1.225;
g=9.8;
q=RHO*Y（2）*Y（2）/2;
if Y（1）    mro=0;
else
    mro=152;
end
if Y（1）    P=76930;
    mc=28.224;
elseif Y（1）>=ts && Y（1）    P=0;
    mc=0;
elseif Y（1）>=tf && Y（1）<145.15
    P=2270;
    mc=1.094;
else 
    P=0;
    mc=0;
end
m=m0-mro-Y（8）;
a=Y（5）-Y（3）;
%插值子程序的调用:cx、cya、mza、mzw、Jz
[cxcyamzamzwJz]=inter（Y（1）Maatf）;
dy（1）=1;
dy（2）=P*cos（a/57.3）/m-cx*q*S/m-g*sin（Y（3）/57.3）;
dy（3）=57.3*（P*sin（a/57.3）+cya*a*q*S-m*g*cos（Y（3）/57.3））/（m*Y（2））;
dy（4）=57.3*（mza*a+mzw*Y（4）*D/（57.3*2*Y（2）））*q*S*D/Jz;
dy（5）=Y（4）;
dy（6）=Y（2）*cos（Y（3）/57.3）;
dy（7）=Y（2）*sin（Y（3）/57.3）;
dy（8）=mc;

%插值子程序
function [cxcyamzamzwJz]=inter（tMaatf）
load data;
A=[0 1 2 3 4 5 6];
MC=[0.0 0.1 0.3 0.5 0.7 0.8 0.9 0.95];
T=[4.1 28.08 56.29 90.14 145.15];
MM=[0.7 0.75 0.8 0.85 0.9 0.95 1.00];
MCA=[0 1 2 3 4 6 8];
if t    Jz=interpt（da7（1:41）da7（1:42）t）;
    Jz=Jz*9.8;
    cya=interpt（da1（:1）da1（:2）Ma）;
    mza=interpt（da1（:1）da1（:3）Ma）;
    mzw=interpt（da1（:1）da1（:4）Ma）;
    cx=interpt2（AMCda2aMa）;
else
    Jz=interpt（da7（5:101）da7（5:102）t）;
    Jz=Jz*9.8;
    cya=interpt（da3（:1）da3（:2）Ma）;
    mza=interpt2（TMMda4tMa）;
    mzw=interpt2（TMMda5tMa）;
    cx=interpt2（MCAMMda6aMa）;
end

%一维线性插值
function b=interpt（xya）;
if（a>=x（1）&&a<=x（end））
   for i=1:length（x）
       if（a>=x（i）&&a<=x（i+1））
           k=i;