%==========================================================================
%
% MOS21
%
% TD : modélisation du frottement de l'air sur un projectile ponctuel
%
%==========================================================================


%=== Nettoyage de Matlab
close all;
clear all;
clc;


%=== Propriétés
m=0.1;      % kg
g=10;       % m/s/s

kS=0.01;     % N/m*s
kN=0.01;     % N/m/m*s*s

x1_0=-2;     % m
x2_0=0.5;   % m
v0=20;      % m/s
alpha=pi/3;
x1p_0=v0*cos(alpha);
x2p_0=v0*sin(alpha);


%=== Solution sans frottement (sf)
dt=0.001;   % s
temps_max=5;% s

n=temps_max/dt;

x1_sf=zeros(n+1,1);
x2_sf=zeros(n+1,1);
for i=1:n+1
    t=(i-1)*dt;
    x1_sf(i,1)=x1p_0*t+x1_0;
    x2_sf(i,1)=-0.5*g*t*t+x2p_0*t+x2_0;
end


%=== Solution avec frottement : modele de Stokes (S)
dt=0.001;   % s
temps_max=5;% s

n=temps_max/dt;

x1_S=zeros(n+1,1);
x2_S=zeros(n+1,1);
for i=1:n+1
    t=(i-1)*dt;
    x1_S(i,1)=-x1p_0*m/kS*exp(-kS/m*t)+x1_0+x1p_0*m/kS;
    x2_S(i,1)=-(x2p_0+m*g/kS)*m/kS*exp(-kS/m*t)-m*g/kS*t+x2_0+m/kS*(x2p_0+m*g/kS);
end


%=== Solution avec frottement : modele de Newton (N)
dt=0.001;   % s
temps_max=5;% s

n=temps_max/dt;

x1_N=zeros(n+1,1);
x2_N=zeros(n+1,1);
x1p_N=zeros(n+1,1);
x2p_N=zeros(n+1,1);

x1_N(1,1)=x1_0;
x2_N(1,1)=x2_0;
x1p_N(1,1)=x1p_0;
x2p_N(1,1)=x2p_0;

for i=1:n
    v=sqrt(x1p_N(i,1)*x1p_N(i,1)+x2p_N(i,1)*x2p_N(i,1));
    x1_N(i+1,1)=x1_N(i,1)+dt*( x1p_N(i,1) + dt/m*(-kN*v*x1p_N(i,1)) );
    x2_N(i+1,1)=x2_N(i,1)+dt*( x2p_N(i,1) + dt/m*(-kN*v*x2p_N(i,1)-m*g) );
    x1p_N(i+1,1)=(x1_N(i+1,1)-x1_N(i,1))/dt;
    x2p_N(i+1,1)=(x2_N(i+1,1)-x2_N(i,1))/dt;
end


%=== Representations des trajectoires
figure(1)
hold on
plot(x1_sf,x2_sf,'-r');
plot(x1_S,x2_S,'-g');
plot(x1_N,x2_N,'-b');
axis equal;
grid;
title('Trajectoire d un projectile sous l action de la pesanteur');
legend('sans frottement','frottement Stokes','frottement Newton');
xlabel('x_1 (m)');
ylabel('x_2 (m)');
hold off