What is wrong with my code?

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]')

[Amit]

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.

Answer 1

[the cyclist]

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.

[Amit]

Thanks a lot. That section of my code seems to be working fine.

I am mainly having problem with the part where my program keeps updating its angle, theta, until the target is overshot. For some reason, it doesn't work properly.

I would be really grateful if you could provide any assistance

[the cyclist]

The test

if x(i+1) ...

is outside your loop over i. Is that what you intended?

[Amit]

yes, it is. Thanks. I just have a problem with the angle. would i need a for loop for the angle to keep increasing until the point at which the target is overshot?

[the cyclist]

Your test of whether x(i+1)>D is embedded inside the test of whether x(i+1)<D, and therefore will never be satisfied. Did you mean this instead?

if x(i+1)<D
    theta=theta+dtheta;
end
if x(i+1)>D
    break
end

That also be written as

if x(i+1)<D
    theta=theta+dtheta;
else
    break
end

[Amit]

unfortunately it still doesn't work. I've made the code slightly simpler and clearer but still i cant find a way of increasing the angle (theta) until x>D.

 %constants
D=10000;
u=600;
m=50;
t_pchute=15;
g=9.81;
a_proj=0.01;
a_pchute=0.05;
C_proj=0.4;
C_pchute=1.2;
p0=1.207;

 %initial conditions
theta=pi/30;
x=0;
i=0;
t=0;
y=0;
dtheta=0.1;
dt=0.1;

 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));

    p=p0*(1-2.333e-5*y).^5;

    Fg=-m*g;
    Fa=-0.5*p*C_proj*a_proj;
    Fp=-0.5*p*C_pchute*a_pchute;

    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

 plot(x,y),grid
 xlabel('Distance/[m]')
 ylabel('Height/[m]')

[the cyclist]

I think this does what you want. I plot once per iteration of the while loop (with a brief pause), so you can see what is happening as theta changes.

%constants
D=10000;
u=600;
m=50;
t_pchute=15;
g=9.81;
a_proj=0.01;
a_pchute=0.05;
C_proj=0.4;
C_pchute=1.2;
p0=1.207;

theta=pi/30;
dtheta = theta/100;  
x=0;

figure
while x(end) < D

theta = theta+dtheta;    

   %initial conditions
  x=0;
  i=0;
  t=0;
  y=0;
  dt=0.1;

   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));
      p=p0*(1-2.333e-5*y).^5;
      Fg=-m*g;
      Fa=-0.5*p*C_proj*a_proj;
      Fp=-0.5*p*C_pchute*a_pchute;
      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

   plot(x,y),grid
   xlabel('Distance/[m]')
   ylabel('Height/[m]')
   pause(0.25)

end

[Amit]

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]

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]

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