參考文獻:
《幾種增大空地飛彈落角的制導方式比較》
《Strategic and tactical missile guidance》
clear
clc
%-----------飛彈參數---------
V_m=260;%飛行速度
X_m=0;
Y_m=1000; %初始飛行高度
theta_m=0*pi/180; %彈道傾角
%----------目标參數---------
V_t=0; %靜止目标
X_t=5000;
Y_t=0;
theta_t=0;
dtheta_t=0;
n_t=0;
A_t=0;
R=sqrt((X_m-X_t)^2+(Y_m-Y_t)^2);
q=atan((Y_t-Y_m)/(X_t-X_m));
dR=((X_m-X_t)*(V_m*cos(theta_m)-V_t*cos(theta_t))+(Y_m-Y_t)*(V_m*sin(theta_m)-V_t*sin(theta_t)))/sqrt((X_m-X_t)^2+(Y_m-Y_t)^2);
dq=((X_t-X_m)*(V_t*sin(theta_t)-V_m*sin(theta_m))-(Y_t-Y_m)*(V_t*cos(theta_t)-V_m*cos(theta_m)))/((X_m-X_t)^2+(Y_m-Y_t)^2);
n_m=-q+theta_m;
c=3;
qf=-90*pi/180;
g=9.8;
n=1;
t=0;
dt=0.01;
while (dR<0)
if R>500
t_go=R/abs(dR);
Am=4*dq*abs(dR)+2*abs(dR)*(q-qf)/t_go+g*cos(theta_m); %考慮重力補償因素的彈道成型制導律
else
Am=4*dq*abs(dR); %比例導引
end
dtheta_m=Am/V_m; %縱向通道:彈道傾角變化函數
theta_m=theta_m+dtheta_m*dt;
%----------------------------計算坐标----------------------------
X_m=X_m+V_m*cos(theta_m)*dt;
Y_m=Y_m+V_m*sin(theta_m)*dt;
alpha=Am/g/(0.3*g);
n_m=-q+theta_m;
R=sqrt((X_m-X_t)^2+(Y_m-Y_t)^2);
q=atan((Y_t-Y_m)/(X_t-X_m));
dR=((X_m-X_t)*(V_m*cos(theta_m)-V_t*cos(theta_t))+(Y_m-Y_t)*(V_m*sin(theta_m)-V_t*sin(theta_t)))/sqrt((X_m-X_t)^2+(Y_m-Y_t)^2);
dq=((X_t-X_m)*(V_t*sin(theta_t)-V_m*sin(theta_m))-(Y_t-Y_m)*(V_t*cos(theta_t)-V_m*cos(theta_m)))/((X_m-X_t)^2+(Y_m-Y_t)^2);
theta_m_store(n)=theta_m; %儲存彈道傾角
Am_store(n)=Am; %儲存縱向過載
alpha_store(n)=alpha;%儲存攻角
P_m_store(:,n)=[X_m;Y_m]; %儲存攔截彈坐标
n=n+1;
t=t+dt;
end
disp('脫靶量為(m):')
R
disp('飛行時間為(s):')
t
figure(1)
plot(P_m_store(1,:),P_m_store(2,:),X_t,Y_t,'r+')
hold on
xlabel('X/m')
ylabel('Y/m')
figure(2)
plot((1:n-1)*dt,Am_store/g)
hold on
xlabel('time/s')
ylabel('Acceleration/g')
title('加速度')
figure(3)
plot((1:n-1)*dt,theta_m_store*180/pi)
hold on
xlabel('time/s')
ylabel('\theta_m/°')
title('彈道傾角')
figure(4)
plot((1:n-1)*dt,alpha_store)
hold on
xlabel('time/s')
ylabel('\alpha/°')
title('攻角')