qwk.m

基于GUI界面铅垂平面弹道仿真。matlab编写。

function varargout = qwk(varargin)
% QWK M-file for qwk.fig
%      QWK, by itself, creates a new QWK or raises the existing
%      singleton*.
%
%      H = QWK returns the handle to a new QWK or the handle to
%      the existing singleton*.
%
%      QWK('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in QWK.M with the given input arguments.
%
%      QWK('Property','Value',...) creates a new QWK or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before qwk_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to qwk_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help qwk

% Last Modified by GUIDE v2.5 31-Aug-2007 21:45:29

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @qwk_OpeningFcn, ...
                   'gui_OutputFcn',  @qwk_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before qwk is made visible.
function qwk_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to qwk (see VARARGIN)
uicontrol(gcf,'Style','Frame',...
   'Units','normalized',...
   'Position',[.84,.018,.002,.923],...
   'Back',[.5 .5 .5]);

% Choose default command line output for qwk
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes qwk wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = qwk_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in push_start.
function push_start_Callback(hObject, eventdata, handles)
% hObject    handle to push_start (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


%Fig=findobj('Type','Figure','Tag','acy');
%hand=getappdata(acy);

%ACY=str2double(get(hand,'string'));
%global H1;
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;
load -ascii ACX.DAT;    %读入数据
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;

%H=get('untitled1',0,'UserData');

%hm=get(untitled1(1),'string');


AE=0.01745329;   		%角度化弧度
h=str2num(get(handles.dt,'string'));;

%fo1=fopen('result.dat','wt');
%fprintf(fo1,'y8        y1         y2         y3         y4         y5         y6         y7      ALPHA  \n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%赋初值
y(1)=str2num(get(handles.vm,'string'));				           %初速
y(2)=str2num(get(handles.thetam,'string'))*AE;                        %初始弹道倾角
y(3)=str2num(get(handles.om,'string'));				           %初始转速
y(4)=str2num(get(handles.phim,'string'))*AE;                           %初始俯仰角
y(5)=str2num(get(handles.xm,'string'));				              %初始x坐标
y(6)=str2num(get(handles.ym,'string'));				           %初始y坐标
y(7)=52.38;	           %起飞质量
y(8)=0   ;                         %初始时间


%fprintf(fo1,'%10.10f %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f\n',...
%             y(8),y(1),y(2)*57.3,y(3)*57.3/360,y(4)*57.3,y(5),y(6),y(7),y(4)-y(2));

while y(6)>=0
   
    rk(8,h);      
    %fprintf(fo1,'%10.10f %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f\n',y(8),y(1),y(2)*57.3,y(3)*57.3/360,y(4)*57.3,y(5),y(6),y(7),y(4)-y(2));

    h1=subplot(handles.h1);plot(y(5),y(6)),hold on,
    h2=subplot(handles.h2);plot(y(8),y(2)*57.3),hold on,
    h3=subplot(handles.h3);plot(y(8),y(1)),hold on,
    h4=subplot(handles.h4);plot(y(8),(y(4)-y(2))*57.3),hold on,
end
reset(h1);
reset(h2);
reset(h3);
reset(h4);

set(get(h1,'xlabel'),'string','x (m)');
set(get(h1,'ylabel'),'string','y (m)');
set(get(h1,'title'),'string','弹道曲线');

set(get(h2,'xlabel'),'string','时间t (s)');
set(get(h2,'ylabel'),'string','弹道倾角\theta (\circ)');
set(get(h2,'title'),'string','弹道倾角随时间的变化');

set(get(h3,'xlabel'),'string','时间t (s)');
set(get(h3,'ylabel'),'string','导弹速度vm (m/s)');
set(get(h3,'title'),'string','导弹速度随时间的变化');

set(get(h4,'xlabel'),'string','时间t (s)');
set(get(h4,'ylabel'),'string','攻角\alpha (\circ)');
set(get(h4,'title'),'string','攻角随时间的变化');
grid(h1);
grid(h2);
grid(h3);
grid(h4);

%fclose(fo1);


%龙格库塔右端函数，dy:导数值；y:积分值

function dy=grkt2f(y)
global H1;
global AJZ;
global AJZT;
global AMa;
global Aa;
global ACX;
global ACY;
global AP;
global APT;
global AXG;
global AXGT;
global AMZWZ1;
global AMZWZ2;
global AMZAF;

RHO=1.225;
SONIC=343.13;

s=0.0227;             %特征面积
q=0.5*RHO*y(1)*y(1);          %动压头
Ma=y(1)/SONIC;          %计算马赫数
L=1.8;                  %特征长度
B=[2.1126 44.0832];      %两级发动机结束工作时间      
g=9.81;


alpha=y(4)-y(2);
ABS_alpha=abs(alpha);
%插值部分开始
P=interp1(APT,AP,y(8))*g;


CX=twointerp(Aa,AMa,ACX,ABS_alpha,Ma);

CY=twointerp(Aa,AMa,ACY,ABS_alpha,Ma);

if alpha<0
    CY=-CY;
end
XF=CX*q*s;
YF=CY*q*s;
wzt=y(3)*L/y(1);

if y(8)<=B(1) 
    MC=2.362;
else  if y(8)>B(2) 
         MC=0;
         P=0;
      else 
         MC=0.21059;
      end
  end
  
if y(8)<B(2)
XG=interp1(AXGT,AXG,y(8));  %插值重心
end

if y(8)<B(2)
  aa=twointerp(Aa,AMa,AMZWZ1,ABS_alpha,Ma);
  bb=twointerp(Aa,AMa,AMZWZ2,ABS_alpha,Ma);
  ab=[aa,bb];
  a=[.9831,.8896];
  MZWZ=interp1(a,ab,XG);
else
    MZWZ=twointerp(Aa,AMa,AMZWZ2,ABS_alpha,Ma);
end

MZAF=twointerp(Aa,AMa,AMZAF,ABS_alpha,Ma); 
if ABS_alpha>0
    MZAF=MZAF*57.3+CY*(XG-0.9381)/(L*alpha); %由于重心的改变，对MZAF进行修正
end
%插值部分结束

JZ=interp1(AJZT,AJZ,y(8));
dy(1)=(P*cos(alpha)-XF-9.81*y(7)*sin(y(2)))/y(7);
dy(2)=(P*sin(alpha)+YF-9.81*y(7)*cos(y(2)))/(y(7)*y(1));
dy(3)=(MZAF*alpha+MZWZ*wzt)*q*s*L/JZ;
dy(4)=y(3);
dy(5)=y(1)*cos(y(2));
dy(6)=y(1)*sin(y(2));
dy(7)=-MC;
dy(8)=1;

%一维插值函数，用于质心、转动惯量及推力的插值，输入变量依次为：待插参数数组，插值点时间值
function W = oneinterp1(t1,t2)
[x y]=size(t1);

if t2<=t1(1,1)%时间小于或者等于最小给定时间值
    j=1;
end    

if t2>t1(x,1)%时间大于最大给定时间值
    j=x-1;
else
     j=1;
     while t2>t1(j+1,1)
          j=j+1;
     end    
end     

W=t1(j,2)+(t1(j,2)-t1(j+1,2))*(t2-t1(j,1))/(t1(j,1)-t1(j+1,1));

%一维插值函数，用于气动参数插值，输入变量依次为：攻角数组，二维待插参数数组，插值点攻角值，二维待插参数数组中要插值的行数
function W = oneinterp2(t1,t3,t4,t6)
[x y]=size(t3);

if t4<=t1(1)%攻角小于或者等于最小给定攻角值
    j=1;
end    

if t4>t1(y)%攻角大于最大给定攻角值
    j=y-1;
else
     j=1;
     while t4>t1(j+1)
          j=j+1;
     end    
end     

W=t3(t6,j)+(t3(t6,j)-t3(t6,j+1))*(t4-t1(j))/(t1(j)-t1(j+1));

%龙格库塔子函数，n:微分方程的个数；h:积分步长

function rk(n,h)

global y;

a(1)=h/2;
a(2)=a(1);
a(3)=h;
a(4)=h;

dy=grkt2f(y);
 
for i=1:n;
	old_y(i)=y(i);
end

for j=1:3;
	for i=1:n;
        y1(i)=old_y(i)+a(j)*dy(i);
		y(i)=y(i)+a(j+1)*dy(i)/3.;
    end
    dy=grkt2f(y1);
end
     
for i=1:n;
    y(i)=y(i)+h*dy(i)/6.;
end

%二维插值函数，输入变量依次为：攻角数组，马赫数数组，待插参数数组，插值点攻角值，插值点马赫数值
%调用了一维插值函数oneinterp2

function F = twointerp(arg1,arg2,arg3,arg4,arg5)
[m n]=size(arg3);
%确定要用到哪二行马赫数
if arg5<=arg2(1)%马赫数小于或者等于最小给定马赫数 
    i=1;
end    

if arg5>arg2(m)%马赫数大于最大给定马赫数
    i=m-1;
else
     i=1;
     while arg5>arg2(i+1)
          i=i+1;
     end    
end     

fy1=oneinterp2(arg1,arg3,arg4,i);
fy2=oneinterp2(arg1,arg3,arg4,i+1);
F=fy1+(fy1-fy2)*(arg5-arg2(i))/(arg2(i)-arg2(i+1));

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double


% --- Executes during object creation, after setting all properties.
function edit2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function edit2_Callback(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double


% --- Executes during object creation, after setting all properties.
function edit3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



function edit3_Callback(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double


% --- Executes during object creation, after setting all properties.
function edit4_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.