c24l3.m

这是本战略战术导弹制导的书中的matlab程序,比书中的 forchan程序简单易懂

count=0;
FR=3.;
DIAM=1.;
XL=20.;
CTW=0.;
CRW=6.;
HW=2.;
CTT=0.;
CRT=2.;
HT=2.;
XN=4.;
XCG=10.;
XHL=19.5;
WACT=150.;
ZACT=.7;
TF=5.;
VM=3000.;
XNCG=10.;
WCR=50.;
ZETA=.7;
TAU=.3;
ALT=0.;
XNT=64.4;
XNP=3.;
VC=4000.;
T1=.1;
T2=.15;
PHIRN=.000002;
PHIGL=20.;
PHIFN=.00000008;
RA=30000.;
R=0.;
MISS=1;
TINT=0.;
A=1000.;
if ALT<=30000.
	RHO=.002378*exp(-ALT/30000.);
else
	RHO=.0034*exp(-ALT/22000.);
end
WGT=1000.;
XNLLIN=0.;
XACC=XCG;
SWING=.5*HW*(CTW+CRW);
STAIL=.5*HT*(CTT+CRT);
SREF=3.1416*DIAM*DIAM/4.;
XLP=FR*DIAM;
SPLAN=(XL-XLP)*DIAM+1.33*XLP*DIAM/2.;
XCPN=2*XLP/3;
AN=.67*XLP*DIAM;
AB=(XL-XLP)*DIAM;
XCPB=(.67*AN*XLP+AB*(XLP+.5*(XL-XLP)))/(AN+AB);
XCPW=XLP+XN+.7*CRW-.2*CTW;
XMACH=VM/A;
XIYY=WGT*(3*((DIAM/2)^2)+XL*XL)/(12*32.2);
TMP1=(XCG-XCPW)/DIAM;
TMP2=(XCG-XHL)/DIAM;
TMP3=(XCG-XCPB)/DIAM;
TMP4=(XCG-XCPN)/DIAM;
B=sqrt(XMACH^2-1);
Q=.5*RHO*VM*VM;
P1=WGT*XNCG/(Q*SREF);
Y1=2+8*SWING/(B*SREF)+8*STAIL/(B*SREF);
Y2=1.5*SPLAN/SREF;
Y3=8*STAIL/(B*SREF);
Y4=2*TMP4+8*SWING*TMP1/(B*SREF)+8*STAIL*TMP2/(B*SREF);
Y5=1.5*SPLAN*TMP3/SREF;
Y6=8*STAIL*TMP2/(B*SREF);
P2=Y2-Y3*Y5/Y6;
P3=Y1-Y3*Y4/Y6;
ALFTR=(-P3+sqrt(P3*P3+4.*P2*P1))/(2.*P2);
DELTR=-Y4*ALFTR/Y6-Y5*ALFTR*ALFTR/Y6;
CNA=2+1.5*SPLAN*ALFTR/SREF+8*SWING/(B*SREF)+8*STAIL/(B*SREF);
CND=8*STAIL/(B*SREF);
CMAP=2*TMP4+1.5*SPLAN*ALFTR*TMP3/SREF+8*SWING*TMP1/(B*SREF);
CMA=CMAP+8*STAIL*TMP2/(B*SREF);
CMD=8*STAIL*TMP2/(B*SREF);
XMA=Q*SREF*DIAM*CMA/XIYY;
XMD=Q*SREF*DIAM*CMD/XIYY;
ZA=-32.2*Q*SREF*CNA/(WGT*VM);
ZD=-32.2*Q*SREF*CND/(WGT*VM);
WZ=sqrt((XMA*ZD-ZA*XMD)/ZD);
WAF=sqrt(-XMA);
ZAF=.5*WAF*ZA/XMA;
XK1=-VM*(XMA*ZD-XMD*ZA)/(1845*XMA);
XK2=XK1;
TA=XMD/(XMA*ZD-XMD*ZA);
XK3=1845*XK1/VM;
W=(TAU*WCR*(1+2.*ZAF*WAF/WCR)-1)/(2*ZETA*TAU);
W0=W/sqrt(TAU*WCR);
Z0=.5*W0*(2*ZETA/W+TAU-WAF^2/(W0*W0*WCR));
XKC=(-W0^2/WZ^2-1.+2.*Z0*W0*TA)/(1.-2.*Z0*W0*TA+W0*W0*TA*TA);
XKA=XK3/(XK1*XKC);
XK0=-W*W/(TAU*WAF*WAF);
XK=XK0/(XK1*(1+XKC));
WI=XKC*TA*W0*W0/(1+XKC+W0^2/WZ^2);
XKR=XK/(XKA*WI);
XKDC=1.+1845./(XKA*VM);
%for R=-.06:.01:.06
%Number of case reduced by factor of 2 to reduce run time
for R=-.06:.02:.06
	X1=0.;
	X2=0.;
	X3=0.;
	X4=0.;
	X5=0.;
	X6=0.;
	X7=0.;
	X8=0.;
	X9=0.;
	X10=0.;
	X11=0.;
	X12=0.;
	X13=0.;
	X14=0.;
	X15=0.;
	if MISS==1
		X3=1.;
	elseif MISS==2
		X8=1.;
	elseif MISS==3
		X6=XNP*VC/32.2;
	end
	T=0;
   %H=.0001;
   %Integration interval increased by factor of 100 to reduce running time
	H=.01;
	S=0;
	TP=T+.00001+TINT;
	while TP<=(TF - 1e-5)
		X1OLD=X1;
		X2OLD=X2;
		X3OLD=X3;
		X4OLD=X4;
		X5OLD=X5;
		X6OLD=X6;
		X7OLD=X7;
		X8OLD=X8;
		X9OLD=X9;
		X10OLD=X10;
		X11OLD=X11;
		X12OLD=X12;
		X13OLD=X13;
		X14OLD=X14;
		X15OLD=X15;
		STEP=1;
		FLAG=0;
		while STEP <=1
			if FLAG==1
				X1=X1+H*X1D;
				X2=X2+H*X2D;
				X3=X3+H*X3D;
				X4=X4+H*X4D;
				X5=X5+H*X5D;
				X6=X6+H*X6D;
				X7=X7+H*X7D;
				X8=X8+H*X8D;
				X9=X9+H*X9D;
				X10=X10+H*X10D;
				X11=X11+H*X11D;
				X12=X12+H*X12D;
				X13=X13+H*X13D;
				X14=X14+H*X14D;
				X15=X15+H*X15D;
				TP=TP+H;
				STEP=2;
			end
			TGO=TP+.00001;
 			X1D=X2;
			X2D=X3;
			Y1PZ=(X6/T2+X5)/T1;
			X3D=Y1PZ/(VC*TGO);
			X4D=-Y1PZ;
			X5D=-Y1PZ+R*Y1PZ;
			Y2PZ=-XKA*WI*X7;
			X6D=-X6/T2+XNP*VC*XKDC*Y2PZ/32.2;
			X7D=XKR*WACT*WACT*X8;
			X8D=X9-2.*ZACT*WACT*X8;
			Y4PZ=XK1*(-32.2*X2-Y2PZ);
			X9D=-WACT*WACT*X8+X10*WAF*WAF-WAF*WAF*Y4PZ/WZ^2;
			Y3PZ=XK3*(X7D+WI*X7+(X4-X5)/57.3);
			X10D=-2.*ZAF*WAF*(X10-Y4PZ/WZ^2)+X11+TA*Y3PZ;
			X11D=-WAF*WAF*(X10-Y4PZ/WZ^2)+Y4PZ+Y3PZ;
			X12D=X1*X1;
			X13D=(Y1PZ/(VC*TGO))^2;
			X14D=(Y1PZ*VC*TGO/RA)^2;
			X15D=Y1PZ^2;
			FLAG=1;
		end;
		FLAG=0;
		X1=.5*(X1OLD+X1+H*X1D);
		X2=.5*(X2OLD+X2+H*X2D);
		X3=.5*(X3OLD+X3+H*X3D);
		X4=.5*(X4OLD+X4+H*X4D);
		X5=.5*(X5OLD+X5+H*X5D);
		X6=.5*(X6OLD+X6+H*X6D);
		X7=.5*(X7OLD+X7+H*X7D);
		X8=.5*(X8OLD+X8+H*X8D);
		X9=.5*(X9OLD+X9+H*X9D);
		X10=.5*(X10OLD+X10+H*X10D);
		X11=.5*(X11OLD+X11+H*X11D);
		X12=.5*(X12OLD+X12+H*X12D);
		X13=.5*(X13OLD+X13+H*X13D);
		X14=.5*(X14OLD+X14+H*X14D);
		X15=.5*(X15OLD+X15+H*X15D);
	end
	XMFN=sqrt(X15*PHIFN);
	XMRN=sqrt(X14*PHIRN);
	XMGL=sqrt(X13*PHIGL);
	XMUDNT=XNT*sqrt(X12/TGO);
	RMS=sqrt(XMFN^2+XMRN^2+XMGL^2+XMUDNT^2);
	count=count+1;
	ArrayR(count)=R;
	ArrayXMFN(count)=XMFN;
	ArrayXMRN(count)=XMRN;
	ArrayXMGL(count)=XMGL;
	ArrayXMUDNT(count)=XMUDNT;
	ArrayRMS(count)=RMS;
end
figure
plot(ArrayR,ArrayXMFN,ArrayR,ArrayXMRN,ArrayR,ArrayXMGL,ArrayR,...
   ArrayXMUDNT,ArrayR,ArrayRMS),grid
xlabel('Radome Slope')
ylabel('Standard Deviation of Miss (Ft) ')
clc
output=[ArrayR',ArrayXMFN',ArrayXMRN',ArrayXMGL',ArrayXMUDNT',ArrayRMS'];
save datfil.txt output /ascii
disp 'simulation finished'