Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

What is wrong with my code?

Asked by Amit on 2 Apr 2013

I am currently doing a school project where i have to plot the trajectory of a projectile launched from the ground with initial speed V0, and angle theta above the horizontal. The projectile has to hit a target distance D=10000m away once it has reached the ground. I have used an initial guess of theta=pi/30 as the angle, so that the projectile does not reach D. The program should then automatically pick a new value of theta. I am doing this by increasing theta from pi/30 in steps dtheta=0.1 until the target is overshot. However, after plotting, I keep getting a wrong and wild graph. Here is my code:

    %contants
    D=10000;
    u=600;
    m=50;
    g=9.81;
    %initial conditions
    theta=pi/30;
    x=0;
    i=0;
    t=0;
    y=0;
    dtheta=0.1;
    dt=0.1;
    while x(i+1)<D
    for i=1:1000
    xdot(i)=u*cos(theta);
    ydot(i)=u*sin(theta);
    v(i)=sqrt(xdot(i).^2+ydot(i).^2);
    x(i+1)=(t.*xdot(i));
    y(i+1)=(t.*ydot(i))-(0.5*g*(t.^2));
    xdbldot(i)=0; 
    ydbldot(i)=-g;
    vel_x(i+1)=xdot(i)+(t.*xdbldot(i));
    vel_y(i+1)=ydot(i)+(t.*ydbldot(i));
    t=t+dt;
    if y(i+1)<0
        break
    end
    end
    if x(i+1)<D
    theta=theta+dtheta;
    if x(i+1)>D
        break
    end
    end
    end
    plot(x,y),grid
    xlabel('Distance/[m]')
    ylabel('Height/[m]')

1 Comment

Amit on 2 Apr 2013

Basically, what I am trying to do is increase the angle, theta, in steps dtheta, until the target, D=10000m, is overshot. I have tried doing this by implementing an if loop, as can be seen above.

However, it just doesn't work, any help will be greatly appreciated.

Amit

Tags

Products

No products are associated with this question.

1 Answer

Answer by the cyclist on 2 Apr 2013
Accepted answer

I think maybe you want this:

x(i+1)=x(i) + (t.*xdot(i));
y(i+1)=y(i) + (t.*ydot(i))-(0.5*g*(t.^2));

when you update x and y.

9 Comments

Amit on 11 Apr 2013

Thank you so much! That worked perfectly! Also allowed me to actually visualise how the trajectory changes as the angle changes. I really appreciate it.

Amit on 13 Apr 2013

If you have some time, could you please help with one final thing? Now, I have incorporated the part where a parachute is deployed at t_pchute=15 seconds. My initial guess for theta is now pi/12 and dtheta=theta/100. At the moment, the trajectory falls short of the D=10000m target. Once again, I would like to know how I can increase the angle theta iteratively in steps dtheta?

I was told to use theta(i+1)=atan(vy(i)/vx(i)) but I have no idea why.

 %Constants
 D=10000;%m
 u=600;%m/s
 m=50;%kg
 t_pchute=15;%s
 g=9.81;%m/s^2
 a_proj=0.01;%m^2
 a_pchute=0.05;%m^2
 C_proj=0.4;
 C_pchute=1.2;
 p0=1.207;%kg/m^3
 theta=pi/12;
 dtheta=theta/100;
 %initial conditions
 x=0;
 i=0;
 t=0;
 y=0;
 dt=0.1;
 p=p0;
 vx=u*cos(theta);
 vy=u*sin(theta);
 v=sqrt(vx.^2+vy.^2);
 for i=1:1000
    Fa(i+1)=0.5*p(i)*C_proj*a_proj*v(i);
    Fp(i+1)=0.5*p(i)*C_pchute*a_pchute*v(i);
    if t<t_pchute
     ax(i+1)=-((Fa(i)*cos(theta(i)))/m)*vx(i); 
     ay(i+1)=-g-(((Fa(i)*sin(theta(i)))/m)*vy(i)); 
    else
     ax(i+1)=-((Fp(i)*cos(theta(i)))/m)*vx(i);   
     ay(i+1)=-g-(((Fp(i)*sin(theta(i)))/m)*vy(i));
    end
    vx(i+1)=vx(i)+(dt.*ax(i));
    vy(i+1)=vy(i)+(dt.*ay(i));
    v(i+1)=sqrt(vx(i).^2+vy(i).^2);
    x(i+1)=x(i)+(dt.*vx(i));
    y(i+1)=y(i)+(dt.*vy(i));
    p(i+1)=p0*(1-2.333e-5*y(i+1)).^5;
    theta(i+1)=atan(vy(i)/vx(i));
    t=t+dt;
    if y(i+1)<0
        break
    end
end
 plot(x,y),grid
 xlabel('Distance/[m]')
 ylabel('Height/[m]')
 title('Projectile Trajectory')
Amit on 13 Apr 2013

I believe that, at the moment, theta changes at each iteration. I would like it to increase by dtheta from the beginning, where x=0 and plot the graph until x(end) and y=0. Whilst x(end) is less than D (10000 metres), theta should increase by dtheta, and the graph should be re-plotted from the the start where x=0. When x(end) reaches D, theta should stop increasing. Not sure how I do this though? Thanks

the cyclist

Contact us