@Diego

Indeed!

So I imagine that I can start from here:

Projectile3[v_,theta_, windspeed_] := First@NDSolve[{

y''[t] == gravity[y[t]] + yDrag[x'[t], y'[t], y[t]]/mass,

y'[0] == v Sin[theta],(*Initial vertical velocity*)

y[0] == 0, (*Initial altitude*)

x''[t] == xDrag[(x'[t] + windspeed), y'[t], y[t]]/mass,

x'[0] == v Cos[theta], (*Initial horizontal velocity*)

x[0] == 0 (*Initial horizontal position*)},

{x[t], y[t]}, {t, 0, 120}]

and create two windspeed functions like:

windspeedX[x_,y_]:=(*horizontal component of the wind speed at point {x,y}*)

windspeedY[x_,y_]:=(*vertical component of the wind speed at point {x,y}*)

And then specify it in the

**x'' [ t ]** and the

**y'' [ t ]**:

Projectile3[v_,theta_, windspeed_] := First@NDSolve[{

y''[t] == gravity[y[t]] + yDrag[x'[t], (y'[t] + windspeedY[x[t], y[t]]), y[t]]/mass,

y'[0] == v Sin[theta],(*Initial vertical velocity*)

y[0] == 0, (*Initial altitude*)

x''[t] == xDrag[(x'[t] + windspeedX[x[t], y[t]]), y'[t], y[t]]/mass,

x'[0] == v Cos[theta], (*Initial horizontal velocity*)

x[0] == 0 (*Initial horizontal position*)},

{x[t], y[t]}, {t, 0, 120}]

It sounds correct. Any thoughts?

What is then the best approach to retrieve the wind velocity at a given point of the domain? I have always struggled with the Interpolation on the 2D domain... Would a package like "NonGridInterpolation" better fit my needs? Are there other alternatives? Please remember that I have black spots (the "buildings").