clear;

global y;
global AMa;
global Aa;
global ACX;
global ACY;
global AP;
global APT;
global AXG;
global AXGT;
global AMZWZ1;
global AMZWZ2;
global AMZAF;
global AJZ;
global AJZT;
global xt0;
global yt0;
global vt;
global q0;
global s0;
global K;
global ammax;


load -ascii ACX.txt;    %讀入數據
load -ascii ACY.txt;
load -ascii AP.txt;
load -ascii APT.txt;
load -ascii AJZ.txt;
load -ascii AJZT.txt;
load -ascii AXG.txt;
load -ascii AXGT.txt;
load -ascii Aa.txt;
load -ascii AMa.txt;
load -ascii AMZAF.txt;
load -ascii AMZWZ1.txt;
load -ascii AMZWZ2.txt;

AE=0.01745329;   		%角度化弧度
h=0.05;
vt=-100;
ammax=10;      %最大過載
thetam0=30*AE;        %初始彈道傾角
rmt0=1000;           %彈目距離
q0=20*AE;            %初始彈目視線角
fo1=fopen('result.dat','wt');
%fprintf(fo1,'y1         y2         y3         y4         y5         y6         y7      y8      y9       y10      x      y      xt      yt  \n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%賦初值
y(1)=500.0;				           %初速
y(2)=thetam0;                        %初始彈道傾角
y(3)=0;				           %初始轉速
y(4)=thetam0;                           %初始俯仰角
y(5)=0;				           %初始x坐標
y(6)=0;				           %初始y坐標
y(7)=52.38;	           %起飛質量
y(8)=0   ;                         %初始時間
y(9)=rmt0;                         %r

y(10)=q0;                           %q

xt0=rmt0*cos(q0);
yt0=rmt0*sin(q0);
s0=thetam0;
K=4;

%fprintf(fo1,'%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',...
            % y(1),y(2)*57.3,y(3)*57.3/360,y(4)*57.3,y(5),y(6),y(7),y(8),y(9),y(10)*57.3,xt0+vt*y(8)-y(9)*cos(y(10)),yt0-y(9)*sin(y(10)),xt0+vt*y(8),yt0);
 fprintf(fo1,'%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',...
             y(8)/h,xt0+vt*y(8)-y(9)*cos(y(10)),yt0-y(9)*sin(y(10)),0,0,0,xt0+vt*y(8),yt0,0,0,0);

while y(9)>=5
   
    rk(10,h);    
    y(11)=xt0+vt*y(8)-y(9)*cos(y(10));
    y(12)=yt0-y(9)*sin(y(10));
    y(13)=xt0+vt*y(8);
    y(14)=yt0;
   % fprintf(fo1,'%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',y(1),y(2)*57.3,y(3)*57.3/360,y(4)*57.3,y(5),y(6),y(7),y(8),y(9),y(10)*57.3,y(11),y(12),y(13),y(14));
    fprintf(fo1,'%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n',y(8)/h,xt0+vt*y(8)-y(9)*cos(y(10)),yt0-y(9)*sin(y(10)),0,0,0,xt0+vt*y(8),yt0,0,0,0);

    plot(y(11),y(12),y(13),y(14)),hold on,
end
%hold on,legend('導彈','目標'),grid on,
fclose(fo1);