Adding air drag to a golf ball trajectory equation

Question

I'm developing a 2D golf game in VB.NET 2005, but I am stuck on how to implement air or wind drag that should affect the ball.

Already I have these equations for projectile:

Vo                            `<-Initial velocity of golfball when hit or fired`

Vx=VoCos(theta)               `<-Vertical component of velocity of golfball`

Vy=VoSin(theta)-gt*           `<-Horizontal component of velocity of golfball`

X=VoCos(theta)t              `<-Vertical dixtance of golfball`

Y=VoSin(theta)t-(0.5)g(t*t)  `<-Horizontal dixtance of golfball`

How do I add air drag to this equation to properly affect the velocity of the golf ball? I don't have any idea how to do it, has anyone worked with similar equations?

Answer 1

I'm not sure if there even exists a closed form for drag or wind, but it is quite easy to simulate in a step-wise fashion (like all the physics libraries do):

1) set your initial condition:

x, y, v_x, v_y , (for t = 0)

2) update position:

x += v_x * dt, y += v_y * dt , (where dt is the time elapsed since the last update)

3) calculate these velocity helpers:

v_squared = v_x * v_x + v_y * v_y and v_length = √(v_squared)

4) calculate drag force:

f_drag = c * v_squared , (where c is the coefficient of friction small!)

5) accumulate forces:

f_x = (-f_drag * v_x / v_length), f_y = (-f_drag * v_y / v_length) + (-g * mass) , (where mass is the mass of your golf ball)

6) update velocity:

v_x += f_x * dt / mass, v_y += f_y * dt / mass

That's basically Euler's Method http://en.wikipedia.org/wiki/Euler_method for approximating those physics.

EDIT

A bit more on how the simulation as requested in the comments:

• The initial condition (t = 0) in your case is

x = 0, y = 0,

v_x = v₀ * cos(θ),

v_y = v₀ * sin(θ)

It's basically the same as in your basic trajectory formula where every occurrence of t is replaced by 0.

• The kinetic energy

KE = 0.5 * m * (V * V) = 0.5 * m * v_squared

is valid for every t (see v_squared as in (3)

• The potential energy

PE = m * g * y

is also always valid.

• If you want to get the current (x,y) for a given t₁ what you need to do is initialize the simulation for t = 0 and do small dt updates until t = t₁

• If you already calculated (x,y) for a t₁ and you want to know their values for a t₂ where t₁ < t₂ all you need to do is calculating those small dt update steps from t₁ to t₂

Pseudo-Code:

simulate(v0, theta, t1)
  dt = 0.1
  x = 0
  y = 0
  vx = v0 * cos(theta)
  vy = v0 * sin(theta)
  for (t = 0; t < t1; t += dt)
    x += vx * dt
    y += vy * dt
    v_squared = vx * vx + vy * vy
    v_length = sqrt(v_squared)
    f_drag = c * v_squared
    f_grav = g * mass
    f_x = (-f_drag * vx / v_length)
    f_y = (-f_drag * vy / v_length) + (-f_grav)
    v_x += f_x * dt / mass
    v_y += f_y * dt / mass
  end for
  return x, y
end simulate