[MATLAB] Javelot

Salut, notre professeur nous propose un projet sous Matlab MAIS LE problème que j'ai affronté c'est que la fonction ne s'affiche pas.

Svp est ce que vous pouvez le courigé si c'etais possible:

% boucle i de balayage des valeurs de dθ/dt(t=0)

for i=-2:2

% boucle m de balayage des valeurs de θ(τ0)

for m=1:17;

y=zeros(100,6);

final=0;

y(1,5)=m*5/180*pi; % angle θ initial (de 5° à 85 ° par pas de 5°)

options = odeset('RelTol',1e-4);

[t,y] = ode45('equations',[0 10],[0 25*cos(y(1,5)) 2 25*sin(y(1,5)) y(1,5) i],options);

e=size(y);

% boucle k enregistre à quel pas de temps k le javelot atteint le sol

for k=1:e(1);

if (y(k,3)>0) & (y(k+1,3)<0)

final=k;

break;

end

end

% fin boucle k

% Le vecteur à 3 dimensions range contiendra les résultats finaux

range(1,m,i+3)=m*5; % angle θ initial en degrés

range(2,m,i+3)=-y(final,3)*(y(final+1,1)-y(final,1))/(y(final+1,3)-y(final,3))+y(final,1);

range(3,m,i+3)=i; % vitesse angulaire dθ/dt(t=0) initiale

axis([0 150 0 50]);

title(['Trajectoire pour une vitesse angulaire \theta/dt=,int2str(i),'rad/sec]);

xlabel(['Portée [m]= ',int2str(range(2,m,i+3))]);

ylabel('Hauteur [m]');

comet(y(1:final+1,1),y(1:final+1,3));

plot(y(1:final+1,1),y(1:final+1,3));

hold on % force chacune des trajectoires à s'imprimer sur la même figure

end

% fin boucle m

hold off

end

% fin boucle i

for i=-2:2

plot(range(1,:,i+3),range(2,:,i+3),'LineWidth',2);

hold on;

end

title('Diagramme Portée en fonction de Angle Attaque \theta ');

xlabel('Angle \theta [°]');

ylabel('Portée[m]');

%Sauvegarde dans le fichier javlot.txt de la portée du javelot en fonction de θ(t=0) pour

chacune des 5 valeurs initiales de dθ/dt(t=0)

fid = fopen('c:\javelot.txt','w');

fprintf(fid,'angle deg; portée ; dteta/dt(t=0)\n');

fprintf(fid,'% 6.2f ; %6.2f ; %6.2f \n',range);

fclose(fid);

hold off;

Le fichier equation.m dont le module de résolution ode45.m a besoin est donné par:

function z = equation(t,y)

z = zeros(6,1);

Delta=0.255;

Speed = sqrtm(y(2)*y(2) + y(4)*y(4));

Phi = atan(y(4)/y(2));

Alpha = sin(y(5) - Phi);

Drag = 0.00024*exp(5.157*abs(Alpha));

Lift = 0.0;

if Alpha > 0.001

Lift = 0.0127*exp(1.34*log(abs(Alpha)));

end

z(1) = y(2);

z(2) = - (Drag*y(2)*Speed + Lift*y(4)*Speed)/0.80625;

z(3) = y(4);

z(4) = - (Drag*y(4)*Speed - Lift*y(2)*Speed)/0.80625 - 9.81;

z(5) = y(6);

z(6) = - Delta*Speed*Speed*(Drag*sin(Alpha) + Lift*cos(Alpha))/0.42;

J'att vos répense