% Kast med liten boll med luftmotstånd.
%
% Dynamiken ges av ODEn
%
% r'' = -g*ez-sigma*r'*|r'|/m
%

m=0.1;       % Massan (kg)

g=9.81;      % Tyngdaccelerationen (m/s^2)

rho=1.2;     % Luftens densitet (kg/m^3)
A=0.05^2*pi; % Kroppens tvärsnittsarea (m^2)
Cd=0.47;     % Luftmotståndskoefficient, 
             % "drag coefficient" (dimensionslös)
             % Läs mer på http://en.wikipedia.org/wiki/Drag_coefficient

sigma = rho*A*Cd/2;

T = 4;       % Sluttid

% Begynnelsevärden

r0 = [0 0 2]'; % Position  
r1 = [20 0 20]'; % Hastighet, kastar i xz-planet -> y(t) konstant

f = @(t,u) [u(4:6); -u(4:6)*norm(u(4:6),2)*sigma/m - [0;0;g]];

[t,Y] = ode45(f,[0 T],[r0;r1]);

subplot(2,1,1)
plot(t,Y(:,1),t,Y(:,3)); grid on
xlabel('t','FontSize',18)
h=legend('x(t)','z(t)','Location','SouthWest'); set(h,'FontSize',18)
axis([0 T -15 50])
title(['\sigma = ' num2str(sigma)],'FontSize',18)

subplot(2,1,2)
plot(Y(:,1),Y(:,3)); grid on
xlabel('x','FontSize',18)
ylabel('z','FontSize',18)
h=legend('z(x)','Location','SouthWest'); set(h,'FontSize',18)
axis([0 45 -15 20])

disp(['sigma = ' num2str(sigma)])
