马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 kkkttt 于 2010-10-10 22:45 编辑
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % PNGE.m
- % Proportional Navigation Guidance Emluator
- % 导弹拦截的比例制导仿真程序
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function [IT,MD,PT]=PNGE(MT,MTSD,ILP,ILD,IV,IAL,S)
- %用比例导引律仿真拦截器拦截来袭导弹的时间步长模拟算法
- % ChengAihua,PLA Information Engineering University,ZhengZhou,China
- % Email:aihuacheng@gmail.com
- % All rights reserved
- %比例导引律:拦截器速度向量的偏转角速度(法向过载)与目标视线的旋转角速度成正比
- %输入参数列表
- % MT Missile Track 来袭弹道导弹运动轨迹(N×3矩阵)
- % MTSD Missile Track's Sample Disdance 导弹轨迹采样时间间隔,兼仿真时间步长
- % ILP Interceptor's Launch Position 拦截器发射地点(1×3矩阵)
- % ILD Interceptor's Launch Direction 拦截器发射方向归一化向量(1×3矩阵)
- % IV Interceptor's Velocity 拦截器飞行速率
- % IAL Interceptor's Acceleration Limit 拦截器(角)加速度上限
- % S Scale 比例导引律中的比例系数
- %输出参数列表
- % IT Interceptor's Track 拦截器的运行轨迹
- % MD Miss Disdance 脱靶量
- % PT Pursuit Time 追踪时间
- % 程序还将绘出拦截器拦截导弹的空间轨迹曲线·
- %第一步:初始化
- XYZ1=MT(1,:);
- xyz1=ILP;
- XYZ2=MT(2,:);
- v1=ILD;
- counter=2
- IT=zeros(0,3);
- IT=[IT;xyz1];
- L1=xyz1-XYZ1;
- LL1=L1.^2;
- L1=L1./sqrt(sum(LL1));
- L2=xyz1-XYZ2;
- LL2=L2.^2;
- L2=L2./sqrt(sum(LL2));
- v2=v1+S.*(L1-L2);
- IA=(IV/MTSD).*(v2-v1);
- IA=sqrt(sum(IA.^2));
- if IA>IAL
- k=IAL/IA;
- v2=v1+(k*S).*(L1-L2);
- end
- vv2=v2.^2;
- v2=v2./sqrt(sum(vv2));
- xyz2=xyz1+(MTSD*IV).*v2;
- IT=[IT;xyz2];
- XYZ1=MT(2,:);
- XYZ2=MT(3,:);
- xyz1=xyz2;
- %下面是迭代过程
- while 1
- %第二步:计算视线方向向量
- L1=xyz1-XYZ1;
- LL1=L1.^2;
- L1=L1./sqrt(sum(LL1));
- L2=xyz1-XYZ2;
- LL2=L2.^2;
- L2=L2./sqrt(sum(LL2));
- %第三步:利用比例引导律计算v2
- v2=v1+S.*(L1-L2);
- %第四步:计算拦截器需要的偏转加速度
- IA=(IV/MTSD).*(v2-v1);
- IA=sqrt(sum(IA.^2));
- %第五步:如果需要的加速度超过拦截器加速度上限,则修正v2向量
- if IA>IAL
- k=IAL/IA;
- v2=v1+(k*S).*(L1-L2);
- end
- %第六步:计算xyz2
- vv2=v2.^2;
- v2=v2./sqrt(sum(vv2));%对v2做归一化处理
- xyz2=xyz1+(MTSD*IV).*v2;
- %第六步:记录和更新
- IT=[IT;xyz2];
- counter=counter+1;
- XYZ1=MT(counter,:);
- XYZ2=MT(counter+1,:);
- xyz1=xyz2;
- if IT(counter,3)>MT(counter,3)||IT(counter,3)>MT(counter-1,3)||IT(counter,3)>MT(counter-2,3)
- Len=size(IT,1);
- P1=IT(Len-1,:);
- P2=IT(Len,:);
- Q1=MT(Len,:);
- Q2=MT(Len+1,:);
- x1=P1(1);y1=P1(2);z1=P1(3);
- x2=Q1(1);y2=Q1(2);z2=Q1(3);
- l1=P1(1)-P2(1);m1=P1(2)-P2(2);n1=P1(3)-P2(3);
- l2=Q1(1)-Q2(1);m2=Q1(2)-Q2(2);n2=Q1(3)-Q2(3);
- A=[l1,m1;l2,m2];
- B=[m1,n1;m2,n2];
- C=[n1,l1;n2,l2];
- M=[x2-x1,y2-y1,z2-z1;l1,m1,n1;l2,m2,n2];
- MD=det(M)/(sqrt((det(A))^2+(det(B))^2+(det(C))^2));
- PT=(Len-1)*MTS*;
- %绘出导弹和拦截器的轨迹图
- %**p=1:10:971
- %plot3(MT(ppp,1),MT(ppp,2),MT(ppp,3),'ro-')
- plot3(MT(:,1),MT(:,2),MT(:,3),'r.-')
- grid on
- hold on
- %qqq=1:10:351
- %plot3(IT(qqq,1),IT(qqq,2),IT(qqq,3),'bo-')
- plot3(IT(:,1),IT(:,2),IT(:,3),'b-')
- return
- end
- end
复制代码 程序作者:aiwa
Email:cheng_aihua@163.com
Q Q:330264704
MSN:aihuacheng@hotmail.com
Google Talk:aihuacheng@gmail.com
版权归原作者所有
|