%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                    航跡發生器
%
%  輸入參數:
%  t          仿真時間
%  T          仿真步長
%  atti       橫滾、俯仰、航向（單位：度）
%  atti_rate  橫滾速率、俯仰速率、航向速率（單位：度/秒）
%  veloB      飛機運動速度——X右翼、Y機頭、Z天向（單位：米/秒）
%  acceB      飛機運動加速度——X右翼、Y機頭、Z天向（單位：米/秒/秒）
%  posi       航跡發生器初始位置經度、緯度、高度（單位：度、度、米）
%  veloN      飛機運動速度——X東向、Y北向、Z天向（單位：米/秒）
%  posiN      經度、緯度、高度（單位：度、度、米）        
%                           程序設計：吳玲  日期：2007/12/05
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%仿真時間設置%%%%%%
t=0;
T=1;
TraceData=[];
t_stop=3600.0;

%%%%%%%%%%%%%%%%航跡發生器%%%%%%%%%%%%%%%%%
atti=zeros(3,1);     %橫滾、俯仰、航向（單位：度）
atti_rate=zeros(3,1);%橫滾速率、俯仰速率、航向速率（單位：度/秒）
veloB=zeros(3,1);    %飛機運動速度——X右翼、Y機頭、Z天向（單位：米/秒）
acceB=zeros(3,1);    %飛機運動加速度——X右翼、Y機頭、Z天向（單位：米/秒/秒）
posi=zeros(3,1);     %航跡發生器初始位置經度、緯度、高度（單位：度、度、米）
posi=[118;32;300.0];  

Re=6378137.0;  %地球半徑（米） 
f=1/298.257;   %地球的橢圓率
posiN=posi;    %初始位置與航跡位置一致
long=posiN(1,1)*pi/180.0;lati=posiN(2,1)*pi/180.0;heig=posiN(3,1);
    %飛行器位置

Rm=Re*(1-2*f+3*f*sin(lati)*sin(lati));
Rn=Re*(1+f*sin(lati)*sin(lati));
    %地球曲率半徑求解
atti(1,1)=0.0;
atti(2,1)=0.0;
atti(3,1)=90.0;  %初始航向角度（單位：度）

veloB(2,1)=0.0;  %飛機初始運動速度——機頭（單位：米/秒）
while t<=t_stop
    %
    if( t==0 )
      acceB(2,1) = 0;               %初始對準時間    
   elseif(t<=20)                    %滑跑
      acceB(2,1)=4.0;
    elseif(t<=24)                   %加速拉起
      atti_rate(2,1)=7.5;
      acceB(2,1)=5.0;            
   elseif (t<=64)                   %爬高
      atti_rate(2,1)=0.0;
   elseif (t<=68)                   %改平
      atti_rate(2,1)=-7.5;
      acceB(2,1)=0.0;
   elseif (t<=668)                  %平直飛行
      atti_rate(2,1)=0.0;
   elseif (t<=671)                  %傾斜預轉彎
      atti_rate(1,1)=10.0;
   elseif (t<=731)                  %轉彎
      atti_rate(3,1)=1.5;
      atti_rate(1,1)=0.0;
   elseif (t<=734)                  %改平
      atti_rate(3,1)=0.0;
      atti_rate(1,1)=-10.0;
   elseif (t<=1334)                 %平直飛行
      atti_rate(1,1)=0.0;
   elseif (t<=1342)                 %加速拉起
      atti_rate(2,1)=7.5;
      acceB(2,1)=2.0;
   elseif (t<=1374)                 %加速爬高
      atti_rate(2,1)=0.0;
   elseif (t<=1382)                 %改平
      atti_rate(2,1)=-7.5;
      acceB(2,1)=0.0;
   elseif (t<=1862)                 %平直飛行
      atti_rate(2,1)=0.0;
   elseif (t<=1902)                 %減速飛行
      acceB(2,1)=-2.5;
   elseif (t<=1905)                 %傾斜預轉彎
      atti_rate(1,1)=10.0;
      acceB(2,1)=0.0;
   elseif (t<=1965)                 %轉彎
      atti_rate(3,1)=1.5;
      atti_rate(1,1)=0.0;
   elseif (t<=1968)                 %改平
      atti_rate(3,1)=0.0;
      atti_rate(1,1)=-10.0;
   elseif (t<=2568)                 %平直飛行
      atti_rate(1,1)=0.0;
   elseif (t<=2574)                 %低頭
      atti_rate(2,1)=-7.5;
   elseif (t<=2594)                 %俯沖
      atti_rate(2,1)=0.0;
   elseif (t<=2600)                 %改平
      atti_rate(2,1)=7.5;
   elseif (t<=2603)                 %傾斜預轉彎
      atti_rate(2,1)=0.0;
      atti_rate(1,1)=10.0;
   elseif (t<=2663)                 %轉彎
      atti_rate(3,1)=1.5;
      atti_rate(1,1)=0.0;
   elseif (t<=2666)                 %改平
      atti_rate(3,1)=0.0;
      atti_rate(1,1)=-10.0;
   elseif (t<=3600)                 %平直飛行
      atti_rate(1,1)=0.0;            
   end
   %
   t=t+T;
   t
   
   veloB(2,1)=veloB(2,1)+acceB(2,1)*T;    
   
   atti(1,1)=atti(1,1)+atti_rate(1,1)*T;
   atti(2,1)=atti(2,1)+atti_rate(2,1)*T;
   atti(3,1)=atti(3,1)+atti_rate(3,1)*T;
  
   roll=atti(1,1)*pi/180.0;pitch=atti(2,1)*pi/180.0;head=atti(3,1)*pi/180.0;

  %坐標系N(地理系)-->B(機體系)
  Cbn=[cos(roll)*cos(head)+sin(roll)*sin(pitch)*sin(head), -cos(roll)*sin(head)+sin(roll)*sin(pitch)*cos(head), -sin(roll)*cos(pitch);
       cos(pitch)*sin(head),                               cos(pitch)*cos(head),                                sin(pitch);
       sin(roll)*cos(head)-cos(roll)*sin(pitch)*sin(head), -sin(roll)*sin(head)-cos(roll)*sin(pitch)*cos(head), cos(roll)*cos(pitch)];

  veloN=Cbn'*veloB;
  Ve=veloN(1,1);Vn=veloN(2,1);Vu=veloN(3,1);

  heig=heig+T*Vu;
  lati=lati+T*(Vn/(Rm+heig));
  long=long+T*(Ve/((Rn+heig)*cos(lati)));

  posiN(1,1)=long*180.0/pi;       %單位：度
  posiN(2,1)=lati*180.0/pi;       %單位：度
  posiN(3,1)=heig;

  
  
  TraceData=[TraceData;t,posiN',veloB(2,1),veloN'];
end
figure(1);
plot3(TraceData(:,2),TraceData(:,3),TraceData(:,4),'m');
title('飛行航跡仿真');
xlabel('經度（^{\circ}）');ylabel('緯度（^{\circ}）');zlabel('高度（m）');
grid;
figure(2);
subplot(3,1,1);plot(TraceData(:,1),TraceData(:,2));title('飛行器經、緯、高度仿真');ylabel('經度（^{\circ}）');grid;
subplot(3,1,2);plot(TraceData(:,1),TraceData(:,3));ylabel('緯度（^{\circ}）');grid;
subplot(3,1,3);plot(TraceData(:,1),TraceData(:,4));xlabel('t(sec)');ylabel('高度（m）');grid;

figure(4);
subplot(3,1,1);plot(TraceData(:,1),TraceData(:,6));title('飛行器東、北、天速度仿真');ylabel('東向速度(m/s)');grid;
subplot(3,1,2);plot(TraceData(:,1),TraceData(:,7));ylabel('北向速度(m/s)');grid;
subplot(3,1,3);plot(TraceData(:,1),TraceData(:,8));xlabel('t(sec)');ylabel('天向速度(m/s)');grid;
%%%%%%%%%%%%%存儲仿真數據%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
save trace1.dat TraceData -ASCII;    %存儲航跡數據