title: Projectile Motion with Drag Force Matlab
body: <p>I'm trying to model projectile motion with air resistance.</p>
  
  <p>This is my code:</p>
  
  <pre><code>function [ time , x , y ] = shellflightsimulator(m,D,Ve,Cd,ElAng)
  
  % input parameters are:
  % m mass of shell, kg
  % D caliber (diameter)
  % Ve escape velocity (initial velocity of trajectory)
  % Cd drag coefficient
  % ElAng angle in RADIANS
  
  A = pi.*(D./2).^2; % m^2, shells cross-sectional area (area of circle)
  rho = 1.2 ; % kg/m^3, density of air at ground level
  h0 = 6800; % meters, height at which density drops by factor of 2
  g = 9.8; % m/s^2, gravity
  dt = .1; % time step
  
  % define initial conditions
  x0 = 0; % m
  y0 = 0; % m
  vx0 = Ve.*cos(ElAng); % m/s
  vy0 = Ve.*sin(ElAng); % m/s
  
  N = 100; % iterations
  
  % define data array
  x = zeros(1,N + 1); % x-position,
  x(1) = x0;
  y = zeros(1,N + 1); % y-position,
  y(1) = y0;
  vx = zeros(1,N + 1); % x-velocity,
  vx(1) = vx0;
  vy = zeros(1,N + 1); % y-velocity,
  vy(1) = vy0;
  
  
  i = 1;
  j = 1;
  
  while i &lt; N
      ax = -Cd*.5*rho*A*(vx(i)^2 + vy(i)^2)/m*cos(ElAng); % acceleration in x
      vx(i+1) = vx(i) + ax*dt; % Find new x velocity
      x(i+1) = x(i) + vx(i)*dt + .5*ax*dt^2; % Find new x position
      ay = -g - Cd*.5*rho*A*(vx(i)^2 + vy(i)^2)/m*sin(ElAng); % acceleration in y
      vy(i+1) = vy(i) + ay*dt; % Find new y velocity
      y(i+1) = y(i) + vy(i)*dt + .5*ay*dt^2; % Find new y position
  
      if y(i+1) &lt; 0 % stops when projectile reaches the ground
          i = N;
          j = j+1;
      else
          i = i+1;
          j = j+1;
      end
  end
  
  plot(x,y,'r-')
  end
  </code></pre>
  
  <p>This is what I am putting into Matlab:</p>
  
  <p>shellflightsimulator(94,.238,1600,.8,10*pi/180)</p>
  
  <p>This yields a strange plot, rather than a parabola. Also it appears the positions are negative values. NOTE: ElAng is in radians!</p>
  
  <p>What am I doing wrong? Thanks!</p>