Why is my projectile motion code only working at certain input angles.

This projectile motion is supposed to project the distance of a baseball with drag consdered. It is not working at some input angles but it is at others. For some reason i can get a good output reading when my input is 51 m/s and 15 deg. Ive gotten a few other random ggod readings but i'm not sure why it's like that so any help would be greatly appreciated.

clear 
 close all
 V = input('input exit velocity [m/s]: '); 
 launch_angle = input('input launch angle [deg]: ');
 Vx = V * cos(launch_angle);
 Vy = V * sin(launch_angle);
 %intial conditions 
  G = 9.80665;             %m/s^2 Acceleration due to Gravity       
  DC = 0.4;                %Drag coefficient
  Area = 0.00426;          %m^2 cross section area of ball
  Mass = 0.149;            %aKg mass of ball
  x(1) = 0;                %intial x postion
  y(1) = 1.2;              %inital y postion
  xf(1) = 0;               %inital xf postion
  yf(1) = 0;               %intial yf postion
  AP = 1.2;                %kg/m^3 Air Density @ Sea  Level
  D = AP*DC*Area/2;        %constant for drag calculations    
  t(1) = 0;                %sets intial time
  dt = 0.01;               %s set the intervals at which time will be evalutated
  i = 1;                   %sets counter/index   
  while min(y)> -0.01;                                          
       t = t + dt;                                               
       i = i + 1;                                                
       xf(i) = xf(i-1)+ Vx.*dt;                                  
       AxD = - ( D / Mass ) * V * Vx                            
       AyD = -G - ( D / Mass ) * V * Vy                         
       Vx_new = Vx + AxD * dt;                                      
       Vy_new = Vy + AyD * dt;                                       
       x(i) = x(i-1) + Vx_new * dt + 0.5 * AxD * dt^2;               
       y(i) = y(i-1) + Vy_new * dt + 0.5 * AyD * dt^2;         
       Vx = Vx_new;
       Vy = Vy_new;
   end;
   proj_distance = -xf(i).*3.28;
 disp(proj_distance)

[Voss]

Use cosd and sind when the angle is in degrees:

launch_angle = input('input launch angle [deg]: ');
Vx = V * cosd(launch_angle);
Vy = V * sind(launch_angle);

And I think you want to convert all xf at the end, not only the last element of xf:

proj_distance = -xf.*3.28;

(This is not an exhaustive list of things that might need to be fixed.)

[Nicola Pugliese]

I think I want the last reading because I only want the output to be the final distance. And yes that definitely had to be fixed with my sin and cos but the outputs are still off for some reason.

[Voss]

That makes sense.

Maybe see if you can get it to look right with no drag (D = 0) first.

Answer 1

[Mathieu NOE]

hello

made a for loop on angles to test the code (for one given V velocity) and it worked as soon as I put the "clear x y" line in the code

clear 
 close all
%  V = input('input exit velocity [m/s]: '); 
%  launch_angle = input('input launch angle [deg]: ');
  V = 10; 
  
  figure(1), hold on
  launch_angle = (5:5:85);
  
  for ci = 1:numel(launch_angle)
     clear x y % here <=
     Vx = V * cosd(launch_angle(ci));
     Vy = V * sind(launch_angle(ci));
     %intial conditions 
      G = 9.80665;             %m/s^2 Acceleration due to Gravity       
      DC = 0.4;                %Drag coefficient
      Area = 0.00426;          %m^2 cross section area of ball
      Mass = 0.149;            %aKg mass of ball
      x(1) = 0;                %intial x postion
      y(1) = 1.2;              %inital y postion
      xf(1) = 0;               %inital xf postion
      yf(1) = 0;               %intial yf postion
      AP = 1.2;                %kg/m^3 Air Density @ Sea  Level
      D = AP*DC*Area/2;        %constant for drag calculations    
      t(1) = 0;                %sets intial time
      dt = 0.01;               %s set the intervals at which time will be evalutated
      i = 1;                   %sets counter/index   
      while min(y)> -0.01                                          
           t = t + dt;                                               
           i = i + 1;                                                
           xf(i) = xf(i-1)+ Vx.*dt;                                  
           AxD = - ( D / Mass ) * V * Vx;                            
           AyD = -G - ( D / Mass ) * V * Vy;                         
           Vx_new = Vx + AxD * dt;                                      
           Vy_new = Vy + AyD * dt;                                       
           x(i) = x(i-1) + Vx_new * dt + 0.5 * AxD * dt^2;               
           y(i) = y(i-1) + Vy_new * dt + 0.5 * AyD * dt^2;         
           Vx = Vx_new;
           Vy = Vy_new;
       end;
       proj_distance = -xf(i).*3.28;
     disp(proj_distance)
     legstr{ci} = (['Angle :' num2str(launch_angle(ci)) ]);
      plot(x,y)
 
  end
  legend(legstr);
  hold off

[Nicola Pugliese]

Thank you this works much better but my results are still a bit inacurrate. Thats most likely a mathematical issue though.

Answer 2

[James Tursa]

Drag depends on current velocity, not initial velocity. So you need to recalculate V at each step. E.g.,

Vx = Vx_new;
Vy = Vy_new;
V = norm([Vx Vy]);

Also it looks like you are double banging the acceleration into the position solution. The acceleration gets into velocity, which you then integrate into delta position via the term. But you also integrate it directly into position via the term. So it is getting into the position solution twice. For now, I would suggest doing simple Euler integration for everything. E.g.,

x(i) = x(i-1) + Vx * dt;
y(i) = y(i-1) + Vy * dt;
Vx = Vx + AxD * dt;                                      
Vy = Vy + AyD * dt;                                       
V = norm([Vx Vy]);

No need for the Vx_new and Vy_new variables.

Once you get the simple Euler scheme working, you can explore more advanced schemes such as Modified Euler, RK4, or even ode45( ).

[Nicola Pugliese]

This works better but the results that I am getting are still very innacurate. Could it be something wrong with the code mathematically? This is where I am at right now with it.

clear 
 close all
 Vmph = input('input exit velocity [mph]: '); 
 V = Vmph./2.237; %convert mph to m/s
 launch_angle = input('input launch angle [deg]: '); 
  
  figure(1), hold on
  
  for ci = 1:numel(launch_angle)
     clear x y 
     Vx = V * cosd(launch_angle(ci));
     Vy = V * sind(launch_angle(ci));
     %intial conditions 
      G = 9.80665;             %m/s^2 Acceleration due to Gravity       
      DC = 0.4;                %Drag coefficient
      Area = 0.00426;          %m^2 cross section area of ball
      Mass = 0.149;            %aKg mass of ball
      x(1) = 0;                %intial x postion
      y(1) = 1.2;              %inital y postion
      xf(1) = 0;               %inital xf postion
      yf(1) = 0;               %intial yf postion
      AP = 1.2;                %kg/m^3 Air Density @ Sea  Level
      D = AP*DC*Area/2;        %constant for drag calculations    
      t(1) = 0;                %sets intial time
      dt = 0.01;               %s set the intervals
      i = 1;                   %sets counter/index   
      while min(y)> -0.01                                          
           t = t + dt;                                               
           i = i + 1;                                                
           xf(i) = xf(i-1)+ Vx.*dt;                                  
           AxD = - ( D / Mass ) * V * Vx;                            
           AyD = -G - ( D / Mass ) * V * Vy;                                                                
           x(i) = x(i-1) + Vx * dt;
           y(i) = y(i-1) + Vy * dt;
           Vx = Vx + AxD * dt;                                      
           Vy = Vy + AyD * dt;                                       
           V = norm([Vx Vy]);
       end;
       xfinal = xf(i).*3.28;
     disp(xfinal)
     legstr{ci} = (['Angle :' num2str(launch_angle(ci)) ]);
      plot(x,y)
 
  end
  legend(legstr);
  hold off

[Mathieu NOE]

hello

what makes you believe there is something wrong mathematically ? the longest range is obtained with an initial angle of 45° which is the correct answer.

I looked again at the drag force equation and the projection on x and y axes are correct (IMO)

beside that , minor upgrade of the code as constants don't need to be reinitialized at every loop

clear 
 close all
%  V = input('input exit velocity [m/s]: '); 
%  launch_angle = input('input launch angle [deg]: ');
  V = 10; 
  
  figure(1), hold on
  launch_angle = (5:5:85);
  
  % constants
      G = 9.80665;             %m/s^2 Acceleration due to Gravity       
      Cd = 0.4;                %Drag coefficient
      Area = 0.00426;          %m^2 cross section area of ball
      Mass = 0.149;            %aKg mass of ball
      rho_air = 1.2;                %kg/m^3 Air Density @ Sea  Level
      D = rho_air*Cd*Area/2;        %constant for drag calculations    
      dt = 0.01;               %s set the intervals at which time will be evalutated

  for ci = 1:numel(launch_angle)
     clear x y
     Vx = V * cosd(launch_angle(ci));
     Vy = V * sind(launch_angle(ci));
     %intial conditions 
      x(1) = 0;                %intial x postion
      y(1) = 1.2;              %inital y postion
      t(1) = 0;                %sets intial time
      i = 1;                   %sets counter/index   
      while min(y)> -0.01                                          
           t = t + dt;                                               
           i = i + 1;                                                
           AxD = - ( D / Mass ) * V * Vx;                            
           AyD = -G - ( D / Mass ) * V * Vy;                         
           Vx = Vx + AxD * dt;                                      
           Vy = Vy + AyD * dt;     
           V = sqrt(Vx.^2 + Vy.^2);
           x(i) = x(i-1) + Vx * dt;               
           y(i) = y(i-1) + Vy * dt;         

      end
      proj_distance = x(i).*3.28;
     disp(proj_distance)
     legstr{ci} = (['Angle :' num2str(launch_angle(ci)) ]);
      plot(x,y)
 
  end
  legend(legstr);
  hold off