Demo entry 6354929

ss

   

Submitted by k on Apr 11, 2017 at 18:10
Language: Matlab. Code size: 2.7 kB.

function bal_homework4
clear;clc;

%ÉèÖóõÖµÒÔ¼°ÆäËûÒÑÖªÁ¿
K=2;
V0=472.12;
Vm=300;
xm0=28000;
ym0=16000;
x0=0;
y0=1000;
R0=sqrt((ym0-y0)^2+(xm0-x0)^2);
q0=atan((ym0-y0)/(xm0-x0));
theta0=30*pi/180;
m0=2022;
time_in=0:0.001:60;
par_in=[R0;q0;theta0;V0;m0;x0;y0];

%Êýֵ΢·Ö
op=odeset('Events',@eventfun);
[time,out]=ode45(@dao,time_in,par_in,op);

%¼ÆËã¹¥½Ç
rhoA=1.225.*exp(-0.00015.*(out(:,7)));
a=K*(out(:,5).*out(:,4)).*(1./out(:,1).*(out(:,4).*sin((out(:,2)-out(:,3)))-Vm*sin(out(:,2)-pi)))+out(:,5).*cos(out(:,3))*9.8;
b=Push(time)+Ya(out(:,4),rhoA)*57.3;
alphaD=a./b*57.3;
%¼ÆËã·¨Ïò¼ÓËÙ¶È
V_w=out(:,4);V_w(end)=[];
time_w=time;time_w(end)=[];
W=V_w.*(diff(out(:,3))./diff(time));

%Êä³ö·ÂÕæ½á¹û
figure(1);  plot(out(:,6),out(:,7),'r-',xm0-Vm*time,ym0+0*time,'g-','LineWidth',1.25);ylim([0,18000]);
            title('µ¼Òýµ¯µÀ');xlabel('ˮƽ¾àÀë/km');ylabel('¸ß¶È/km');grid on;
            legend('µ¼µ¯','Ä¿±ê','Location','best');
figure(2);  plot(time,alphaD,'b-','LineWidth',1.25);
            title('¹¥½Ç±ä»¯');xlabel('ʱ¼ä/(Ãë)');ylabel('¹¥½Ç/(¶È)');grid on;
figure(3);  plot(time_w,W,'m-','LineWidth',1.25);
            title('·¨Ïò¼ÓËٶȱ仯');xlabel('ʱ¼ä/(s)');ylabel('·¨Ïò¼ÓËÙ¶È/(m/s^2)');grid on;
figure(4);  plot(time,out(:,4)./340,'c-','LineWidth',1.25);
            title('Ëٶȱ仯');xlabel('ʱ¼ä/(Ãë)');ylabel('ËÙ¶È/(ÂíºÕ)');grid on;
% set([1,2,3,4],'Position',get(0,'ScreenSize'));  %ͼÏñ´°¿Ú×î´ó»¯
% saveas(1,'µ¼Òýµ¯µÀ','png');saveas(2,'¹¥½Ç±ä»¯','png');saveas(3,'·¨Ïò¼ÓËٶȱ仯','png');saveas(4,'Ëٶȱ仯','png'); %±£´æÇúÏß            
            
%--------------------------------------------------------------------------
% µ¼µ¯Ô˶¯·½³Ì
function dy=dao(t,y)
rho=1.225*exp(-0.00015*y(7));
eta=y(2)-y(3);
etam=y(2)-pi;
P=Push(t);
G=y(5)*9.8;
alpha=(K*y(5)*y(4)*(1/y(1)*(y(4)*sin(eta)-Vm*sin(etam)))+G*cos(y(3)))/(P+Ya(y(4),rho)*57.3);
X=Drag(y(4),alpha,rho);
Y=Ya(y(4),rho)*alpha*57.3;
dy=[-y(4)*cos(eta)+Vm*cos(etam);
    1/y(1)*(y(4)*sin(eta)-Vm*sin(etam));
    1/(y(5)*y(4))*(P*sin(alpha)+Y-G*cos(y(3)));
    1/y(5)*(P*cos(alpha)-X-G*sin(y(3)));
    -Mc(t);
    y(4)*cos(y(3));
    y(4)*sin(y(3))];
end

%--------------------------------------------------------------------------
%ÉèÖÃÊýֵ΢·Öµü´ú½áÊøÌõ¼þ
function [value,isterminal,direction] = eventfun(~,y)   %ÉèÖõü´ú½áÊøÌõ¼þ
        value=1;
        if y(1)<20
            value=0;       %´¥·¢Ê±¼ä£¬µ±ÆäֵΪ0µÄʱºò£¬Ê±¼ä»á´¥·¢
        end
        isterminal=1;    %ÉèΪ1ʱ»á£¬´¥·¢Ê±¼ä»áÍ£Ö¹Çó½âÆ÷£¬Éè0ʱ´¥·¢²»Ó°Ï칤×÷
        direction=-1;    %´¥·¢·½ÏòÉè1ʱÊÇÉÏÉý´¥·¢£¬Éè-1ÊÇϽµ´¥·¢£¬Éè0ÊÇË«Ïò´¥·¢
end
end

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).