铅垂平面飞行弹道仿真及分析

-0.5738 -0.5871 -0.6014 -0.6153 -0.6287 -0.6415; -0.6272 -0.6407 -0.6553 -0.6694 -0.6830 -0.6960]; %%当xg=.8896时 %其他参数

S=0.0227; L=1.8; SONIC=343.13; RHO=1.225;

%四阶龙格-库塔法子函数 function rk_4(n,h) global y dy; dy=zeros(n,1); old_y=zeros(1,n); y1=zeros(1,n); a(1)=h/2; a(2)=h/2; a(3)=h; a(4)=h;

dery(y); for i=1:n

old_y(i)=y(i); end

for j=1:3

for i=1:n

y1(i)=old_y(i)+a(j)*dy(i); y(i)=y(i)+a(j+1)*dy(i)/3; end

dery(y1); end

for i=1:n

y(i)=y(i)+a(1)*dy(i)/3; end

%右端子函数 function dery(y)

global dy ;

global L S SONIC RHO;

global ma abs_alpha cx cy mzaf mzwz jz alpha p mc; aa=zeros(1,4);

q=RHO*y(2)*y(2)/2; ma=y(2)/SONIC; alpha=y(5)-y(3);

abs_alpha=abs(alpha); interp();

aa(1)=sin(y(3)*pi/180);aa(2)=cos(y(3)*pi/180);

aa(3)=sin(alpha*pi/180);aa(4)=cos(alpha*pi/180); if alpha<0 cy=-cy;end

xf=cx*q*S;yf=cy*q*S; wzt=y(4)*L/y(2); dy(1)=1;

dy(2)=(p*aa(4)-xf-9.81*y(8)*aa(1))/y(8);

dy(3)=(p*aa(3)+yf-9.81*y(8)*aa(2))/(y(2)*y(8))/pi*180;%转换成度 dy(4)=(mzaf*alpha+mzwz*wzt)*q*S*L/jz; dy(5)=y(4)/pi*180; dy(6)=y(2)*aa(2); dy(7)=y(2)*aa(1); dy(8)=-mc;

%插值子函数

function interp()

global acx acy ajz amzaf amzwz axg ap amc agc andm andaf b L ; global ma abs_alpha cx cy mzaf mzwz xg jz p mc; global y; if y(1)

%a=interp11(tt,3,y(1));

xg=interp11(axg,10,y(1)); jz=interp11(ajz,14,y(1)); p=interp11(ap,11,y(1)); mc=interp11(amc,6,y(1)); else xg=axg(2,10); jz=ajz(2,14); p=ap(2,11); mc=amc(2,6);

end

cx=interp33(ma,abs_alpha,9,6,andm,andaf,acx); cy=interp33(ma,abs_alpha,9,6,andm,andaf,acy);

mzaf0=interp33(ma,abs_alpha,9,6,andm,andaf,amzaf);

if abs_alpha~=0.

mzaf=mzaf0+cy*(xg-agc(1))/(abs_alpha*L);%去掉了*RAD else mzaf=mzaf0;

end

mzwz=interp33(ma,abs_alpha,9,6,andm,andaf,amzwz);

%不等距单变元线性插值子函数 function res=interp11(yy,n,x) for j=1:(n-1)

if x<=yy(1,j+1) i=j;break; else i=n-1; end end

res=yy(2,i)+(yy(2,i+1)-yy(2,i))*(x-yy(1,i))/(yy(1,i+1)-yy(1,i));

%等距双变元抛物线线性插值子函数 function res=interp33(x,qq,n1,n2,a,bb,yy) h=(bb(2)-bb(1))/(n2-1); i=fix((qq-bb(1))/h+1); if (i-1)<0 i=1;

elseif (i-n2)>=0 i=n2-1;

end

yy1=interp31(x,n1,i,a,yy); yy2=interp31(x,n1,i+1,a,yy);

res=((qq-(i-1)*h-bb(1))*yy2-(qq-i*h-bb(1))*yy1)/h;

%等距单变元抛物线插值子函数 function res=interp31(x,n1,i,a,yy) h=(a(2)-a(1))/(n1-1); i1=fix((x-a(1))/h+1); if (i1-1)<=0

i1=1;

else if (i1-n1)>=0 i1=n1-2;end end

x0=a(1)+(i1-1)*h; x1=x0+h; x2=x0+2*h;

c0=0.5*(x-x1)*(x-x2)/(h^2); c1=-(x-x0)*(x-x2)/(h^2); c2=0.5*(x-x0)*(x-x1)/(h^2);

res=c0*yy(i1,i)+c1*yy(i1+1,i)+c2*yy(i1+2,i);

%保存节点处数据

function result(k)

global R1 R2 R3 R4 R5 R6 R7 R8 R9; global y;

R1(k)=y(1); R2(k)=y(2); R3(k)=y(3); R4(k)=y(4); R5(k)=y(5); R6(k)=y(6); R7(k)=y(7); R8(k)=y(8);

R9(k)=y(5)-y(3);

%输出数据文件子函数

function savedata(ii)

global R1 R2 R3 R4 R5 R6 R7 R8 R9;

fp=fopen('c:\\\\data.dat','w');

fprintf(fp,'T V SIGOMA WZ E X Y M ALPHA\\n'); for i=1:ii

fprintf(fp,'%f %f %f %f %f %f %f %f %f\\n',R1(i),R2(i),... R3(i),R4(i),R5(i),R6(i),R7(i),R8(i),R9(i)); end fclose(fp);

%绘图子函数

function drawing()

global R1 R2 R3 R4 R5 R6 R7 R8 R9; figure('numbertitle','on','name','弹道曲线') plot(R6,R7); title('弹道曲线')

xlabel('x/m'); ylabel('y/m');

figure('numbertitle','on','name','速度曲线') plot(R1,R2); title('速度曲线') xlabel('T/s');

ylabel('V/(m/s)');

figure('numbertitle','on','name','攻角曲线') plot(R1,R9); title('攻角曲线') xlabel('T/s');

ylabel('alpha/度');