Projectile Motion Simulation

The vacuum parabola

With gravity as the only force, the horizontal velocity never changes and the vertical velocity falls by g every second. For launch speed v₀ at angle θ from height h₀:

From ground level (h₀ = 0) this gives the familiar results:

sin 2θ peaks at θ = 45°, so that is the best angle — but only from ground level. Launched from a height, the projectile has extra time to travel while it drops, and the best angle falls to tan θ* = v₀ / √(v₀² + 2gh₀). Energy conservation fixes the impact speed at √(v₀² + 2gh₀) whatever the angle. From ground level, switching to the Moon multiplies every distance and time by 9.81/1.62 ≈ 6.1 and leaves the shape unchanged.

Adding quadratic air drag

For a ball moving fast through air, drag grows with the square of speed and points against the velocity:

Here ρ is the air density, C_d the drag coefficient, A the cross-section and m the mass. For a baseball-sized ball (ρ = 1.2 kg/m³, C_d = 0.35, A = 43 cm², m = 145 g) this gives k ≈ 0.006 m⁻¹. Because |v| couples the horizontal and vertical motion, there is no closed-form solution; the lab integrates the equations numerically.

What changes, and why:

• Lower and shorter. Drag removes energy all the way, so the peak is lower and the ball lands much closer; compare the two curves.
• Lopsided. The ball is fast on the way up and slow on the way down, so drag bends the descent more: it falls more steeply than it rose (the dots bunch up on the way down), and the impact angle exceeds the launch angle.
• Terminal speed. A falling ball stops accelerating when drag balances gravity, at v_t = √(g/k). A ball that falls far enough approaches that speed however hard it was thrown.
• A flatter best angle. Time spent high up is time spent losing speed, so a lower trajectory pays off. The angle scan finds the optimum; it keeps falling as k or v₀ increases.

Keep the setting physical: the Moon has no air, so k = 0 there, and the air on Mars is less than 2% as dense as Earth's at the surface, so k is correspondingly smaller.

How the simulation computes it

The vacuum trajectory is drawn from the closed-form solution. The drag trajectory is integrated with the classical fourth-order Runge–Kutta method (RK4) on the state (x, y, v_x, v_y), using at least 1,000 steps per vacuum flight time and smaller steps when drag is strong. The landing is found by shrinking the last step with the secant method until the height is zero, so the range is not rounded to a whole step. With k = 0 the integrator reproduces the closed-form numbers, a built-in check.

The best angle with drag is found by computing the range at every whole degree from 1° to 89°, then refining around the best one with a golden-section search.

Example result

A ball launched at 40 m/s and 45° from the ground on Earth would travel 163.1 m in a vacuum, reaching 40.8 m and landing after 5.77 s at 40.0 m/s. With k = 0.005 m⁻¹ it reaches only 31.2 m, lands after 5.03 s at 103.3 m, and hits the ground at 27.0 m/s, descending at 55.8°. The best angle with drag is 41.5° instead of 45.0°, for a range of 104.0 m.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator. The excerpts show the core of each.

// Quadratic air drag: a = (0, −g) − k·|v|·v, with k = ρ·C_d·A / (2m) in 1/m.

export function accel(vx, vy, g, k) {
  const s = Math.sqrt(vx * vx + vy * vy);
  return [-k * s * vx, -g - k * s * vy];
}

// One RK4 step for the state [x, y, vx, vy].

export function rk4Step(st, dt, g, k) {
  const [x, y, vx, vy] = st;
  const a1 = accel(vx, vy, g, k);
  const a2 = accel(vx + (dt / 2) * a1[0], vy + (dt / 2) * a1[1], g, k);
  const a3 = accel(vx + (dt / 2) * a2[0], vy + (dt / 2) * a2[1], g, k);
  const a4 = accel(vx + dt * a3[0], vy + dt * a3[1], g, k);
  return [
    x + (dt / 6) * (vx + 2 * (vx + (dt / 2) * a1[0]) + 2 * (vx + (dt / 2) * a2[0]) + (vx + dt * a3[0])),
    y + (dt / 6) * (vy + 2 * (vy + (dt / 2) * a1[1]) + 2 * (vy + (dt / 2) * a2[1]) + (vy + dt * a3[1])),
    vx + (dt / 6) * (a1[0] + 2 * a2[0] + 2 * a3[0] + a4[0]),
    vy + (dt / 6) * (a1[1] + 2 * a2[1] + 2 * a3[1] + a4[1]),
  ];
}

// Integrate until the projectile comes back to y = 0, with at least `steps` steps per vacuum

// flight time and k·v0·dt ≤ 0.05 so strong drag is resolved. The last, partial step is

// found by the secant method on its length, so the landing point is as accurate as RK4.

export function fly(v0, th, h0, g, k, steps, path) {
  const T0 = vacuum(v0, th, h0, g).T;
  const n = Math.max(steps, Math.ceil((T0 * k * v0) / 0.05));
  const dt = T0 / n;
  let st = [0, h0, v0 * Math.cos(th), v0 * Math.sin(th)], t = 0, top = h0;
  if (path) path.push(st);
  for (let i = 0; i < 100 * n; i++) {
    const next = rk4Step(st, dt, g, k);
    if (next[1] < 0) {
      let a = 0, ya = st[1], b = dt, yb = next[1], land = next;
      for (let j = 0; j < 6 && yb !== ya; j++) {
        const c = b - (yb * (b - a)) / (yb - ya);
        land = rk4Step(st, c, g, k);
        a = b; ya = yb; b = c; yb = land[1];
      }
      t += b; st = [land[0], 0, land[2], land[3]];
      break;
    }
    st = next; t += dt;
    if (st[1] > top) top = st[1];
    if (path) path.push(st);
  }
  if (path) path.push(st);
  return { T: t, range: st[0], height: top, speed: Math.sqrt(st[2] * st[2] + st[3] * st[3]), angle: Math.atan2(-st[3], st[2]) / DEG, dt };
}

def accel(vx, vy, g, k):
    """Quadratic air drag: a = (0, -g) - k*|v|*v, with k = rho*C_d*A/(2m) in 1/m."""
    s = math.sqrt(vx * vx + vy * vy)
    return (-k * s * vx, -g - k * s * vy)


def rk4_step(st, dt, g, k):
    x, y, vx, vy = st
    a1 = accel(vx, vy, g, k)
    a2 = accel(vx + dt / 2 * a1[0], vy + dt / 2 * a1[1], g, k)
    a3 = accel(vx + dt / 2 * a2[0], vy + dt / 2 * a2[1], g, k)
    a4 = accel(vx + dt * a3[0], vy + dt * a3[1], g, k)
    return (
        x + dt / 6 * (vx + 2 * (vx + dt / 2 * a1[0]) + 2 * (vx + dt / 2 * a2[0]) + (vx + dt * a3[0])),
        y + dt / 6 * (vy + 2 * (vy + dt / 2 * a1[1]) + 2 * (vy + dt / 2 * a2[1]) + (vy + dt * a3[1])),
        vx + dt / 6 * (a1[0] + 2 * a2[0] + 2 * a3[0] + a4[0]),
        vy + dt / 6 * (a1[1] + 2 * a2[1] + 2 * a3[1] + a4[1]),
    )


def fly(v0, th, h0, g, k, steps, path=None):
    """Integrate to y = 0; the final partial step is found by the secant method."""
    T0 = vacuum(v0, th, h0, g)["T"]
    n = max(steps, math.ceil(T0 * k * v0 / 0.05))
    dt = T0 / n
    st, t, top = (0, h0, v0 * math.cos(th), v0 * math.sin(th)), 0, h0
    if path is not None:
        path.append(st)
    for _ in range(100 * n):
        nxt = rk4_step(st, dt, g, k)
        if nxt[1] < 0:
            a, ya, b, yb, land = 0, st[1], dt, nxt[1], nxt
            j = 0
            while j < 6 and yb != ya:
                c = b - yb * (b - a) / (yb - ya)
                land = rk4_step(st, c, g, k)
                a, ya, b, yb = b, yb, c, land[1]
                j += 1
            t += b
            st = (land[0], 0, land[2], land[3])
            break
        st = nxt
        t += dt
        if st[1] > top:
            top = st[1]
        if path is not None:
            path.append(st)
    if path is not None:
        path.append(st)
    return {"T": t, "range": st[0], "height": top, "speed": math.sqrt(st[2] * st[2] + st[3] * st[3]),
            "angle": math.atan2(-st[3], st[2]) / DEG, "dt": dt}

Run the full model from a terminal:

node simulation.js --v0 40 --angle 45
python simulation.py --v0 40 --angle 45

Download the model

Assumptions & limitations

• Point mass. No spin, so no Magnus lift or curving; no tumbling.
• Constant drag coefficient. Real C_d changes with speed (for example, the drag crisis of a baseball or golf ball), and linear drag dominates for very small or slow objects.
• Still, uniform air. No wind, and air density does not change with altitude.
• Flat ground, uniform gravity. Landing is at y = 0; the Earth's curvature and rotation are ignored, which is fine for ranges of a few kilometers or less.

Related simulations

• OS-035Harmonic oscillatorA mass on a spring with friction and a periodic push: damping regimes, phase portraits and resonance.
• OS-036Double pendulumTwo nearly identical double pendulums that part ways: sensitive dependence on initial conditions.
• OS-012Random walkCoin-flip steps on a line or grid: distance grows like √n, and in 1-D and 2-D every walker comes back.

Background: Euler vs Runge–Kutta: Time Stepping in Simulations

References

1. Taylor, J. R. (2005). Classical Mechanics, chapter 2 (projectiles and air resistance). University Science Books.
2. Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007). Numerical Recipes, 3rd ed., chapter 17. Cambridge University Press.