Z Skrypty dla studentów Ekonofizyki UPGOW

Aktualna wersja na dzień 21:34, 4 sty 2010

Modelowanie: Armata

Mamy: \[\vec F = m \vec a\] więc: \[\vec F = m \vec \ddot x\] w kartezjańskim układzie współrzędnych mamy: \[\begin{cases} F_x = m \ddot x \\ F_y = m \ddot y \end{cases}\] Siły w ogólności zależą od prędkości, czasu i położenia. W przypadku lotu pocisku możemy założyć, że działa na niego siła tarcia oraz siła ciężkości.

\[\vec F = \vec T + \vec Q \] Gdzie Q to ciężar \[\vec Q = \begin{cases} 0 \\ -mg \end{cases}\]

Siła tarcia zależy od prędkości ruchu posicku względem powietrza. Jeżeli będzie wiał wiatr \[\vec u(t) = \begin{cases} u(t) \\0 \end{cases}\] to trzeba jego prękość odjąć od prędkości ruchu pocisku względem ziemi (wiatr pozorny). \[\vec T = - \gamma (\vec v-\vec{u(t)}) = \begin{cases} -\gamma (v_x-u(t) ) \\ -\gamma v_y \end{cases} \]

Czyli mamy układ równań:

\[\begin{cases} \dot x=v_x \\ \dot y=v_y\\ \dot v_x = -\gamma/m \; (v_x-u(t)) \\ \dot v_y=-\gamma/m\; v_y-g \end{cases}\]

Kod źródłowy

clear all 
close all
disp('This is a gun!')
 
function dx = gun(X,T)
    g=10; m=1; gama=1.5;
	 u=(T>0.1)*10; 
    dx = zeros(4,1);
    dx(1) = X(3);      % x
    dx(2) = X(4);       % y
    dx(3) = -gama/m*( X(3)-u );    % vx
    dx(4) = -gama/m*X(4) - g;	  % vy
    return
end
 
n=125;
t = linspace(0,1.22,n); 
i=1;
subplot(2,1,1)
hold on 
alfas=linspace(0.01,pi/2,10);
v0=1.;
for alfa=alfas
	X0 = [0,0,v0*cos(alfa),v0*sin(alfa)];
	[X,T,MSG]=lsode(@gun,X0,t);
 	plot(X(:,1),X(:,2),'r-')
	i0=find(X(:,2)<0)(1)-1;
	i1=find(X(:,2)<0)(1);
	x_interp(i)=X(i0,1)+( 0 -X(i0,2))/(X(i1,2)-X(i0,2))*(X(i1,1)-X(i0,1));
	plot(x_interp(i), 0 ,'g*')
	i+=1;
endfor
 
hold off
grid on
ylim([0,.05])
xlim([0,.1])
subplot(2,1,2)
plot(alfas*180/pi,x_interp,'bo')
alfa_fine=linspace(0,90,1000);
x_fine=interp1(alfas*180/pi,x_interp,alfa_fine,'cubic');
[zasieg,alfa_i]=max(x_fine);
hold on
plot(alfa_fine,x_fine)
plot(alfa_fine(alfa_i),zasieg,'ro')
printf("Najwiekszy zasieg dla strzalu pod katem %f stopni\n",alfa_fine(alfa_i))