# Computing launch parameters for hitting a point in 3D with projectile under influence of wind

The end goal of this problem is to compute functions which describe the launch parameters which are needed to hit a specified goal in 3D in the presence of wind disturbances. This is as far I have come ($\theta$ and $\varphi$ are variables from the spherical coordinate system):

g = 9.81;
m = 10;
rho = 1.225;
Cd = 0.5;
A = 0.1;
theta = 45 Degree;(*0-180 degrees*)
phi = 45 Degree;(*0-360 degrees*)
v = 35;(*Exit velocity*)

(*wind speeds*)
Ux = -7;
Uy = -3;
Uz = 0;

(*Calculations from classic projectile motion*)
EndTime = (2 v Sin[theta])/g + 5;
h = (v^2 Sin[theta]^2)/(2 g) + 5;
d = v^2/g Sin[2 theta] + 5;

(***SOLVE***)

(*No drag*)
EquationsNoDrag =
NDSolve[{m z''[t] == -m g, z == 0, z' == v Cos[theta],
m x''[t] == 0, x == 0, x' == v Sin[theta] Cos[phi],
m y''[t] == 0, y == 0, y' == v Sin[phi] Sin[theta],
WhenEvent[
z[t] == 0, {tMaxNoDrag = t, "StopIntegration", z'[t] -> 0,
y'[t] -> 0, x'[t] -> 0}]}, {x[t], y[t], z[t]}, {t, 0, EndTime}];

(*With drag*)
EquationsDrag =
NDSolve[{m z''[t] == -m g -
Tanh[z'[t]]  1/2 rho Cd A (z'[t] - Uz)^2, z == 0,
z' == v Cos[theta],
m x''[t] == -Tanh[x'[t]]  1/2 rho Cd A (x'[t] - Ux)^2, x == 0,
x' == v Sin[theta] Cos[phi],
m y''[t] == -Tanh[y'[t]]  1/2 rho Cd A (y'[t] - Uy)^2, y == 0,
y' == v Sin[phi] Sin[theta],
WhenEvent[z[t] == 0, {tMaxDrag = t, "StopIntegration"}]}, {x[t],
y[t], z[t]}, {t, 0, EndTime}];

(***PLOTTING***)

plots = {x[t], y[t], z[t]} /. {EquationsNoDrag, EquationsDrag};
LandingPointNoDrag = {x[t], y[t], z[t]} /. EquationsNoDrag /.
t -> tMaxNoDrag // Chop // Flatten;
LandingPointDrag = {x[t], y[t], z[t]} /. EquationsDrag /.
t -> tMaxDrag // Chop // Flatten;
SetPoint = {75, 75, 0};
Error = SetPoint - LandingPointDrag;
Show[ParametricPlot3D[
plots, {t, 0, If[tMaxDrag > tMaxNoDrag, tMaxDrag, tMaxNoDrag]},
PlotRange -> All, PlotTheme -> "Classic",
PlotLegends -> {"No drag", "Drag"}, Boxed -> False,
AxesOrigin -> {0, 0, 0}, AxesLabel -> {"x", "y", "z"},
PlotLabel -> v "m/s exit velocity"],
ListPointPlot3D[{LandingPointDrag, LandingPointNoDrag, SetPoint},
PlotStyle -> {Black, PointSize -> Large}],
Arrow[Tube[{{0, 0, 0}, {Ux, Uy, Uz} d/40}, 0.4]], Red, Thick,
Dashed, Line[{LandingPointDrag, SetPoint}]}], ImageSize -> Large]


And I get a plot like this, where the dashed line is the distance between the current landing point and an arbitrary setpoint which I want to hit, and the arrow is the wind vector (the ideal curve is included for reference): This is where things get tricky. I was thinking maybe to map a great number of landing points and make a function where the landing point in the $xy$ plane is a function of the wind speeds $U_x$, $U_y$, exit velocity $v$, and the angles $\theta$ and $\varphi$ using FindFit. I am not nearly sure how to attack this, however. The final goal is to have some estimated functions where the angles and initial speed are given like functions of the wind speeds, for example polynomial approximations, such that the setpoint is hit as accurately as possible (see equations below). This makes it possible to estimate the necessary configuration easily on other platforms without the need to solve DE and so forth. I would like to be able to calculate the parameters $v$, $\theta$ and $\varphi$ on other platforms which only have access to basic math operations, e.g. an advanced calculator. The functions could be 100th degree polynomials, it doesn't matter, but they will have to be stand alone without the need for solving DEs. I've attempted to clarify this further with the figure below (SP=setpoint).

$$\theta=f(U_x,U_y)$$ $$\varphi=g(U_x,U_y)$$ $$v=h(U_x,U_y)$$ For convenience, one may put on some constraints, such as $(U_x,U_y)\in [-10,10]m/s$, $v\in [0,40]m/s]$ and $\theta\in [20,70]^\circ$ or something similar.

Any help or alternative solution methods are greatly appreciated.

• (at) Julian: Interesting Project. Maybe it's no a bad idea to simplify the Problem drastically, e.g. to y==0, Uy = 0, and study the evaluation and presentation of the results. Then you can generalize it. – Dr. Wolfgang Hintze Feb 26 '15 at 11:03

You may use NMinimize[] on the results of ParametricNDSolve[] like this:

g = 9.81; m = 10; rho = 1.225; Cd = 0.5; A = 0.1; rcd = rho Cd A;
vMax = 40;
EndTime[theta_] := (2 vMax Sin[theta])/g + 5;

sol[Ux_, Uy_, Uz_] :=
Quiet@ParametricNDSolve[{
m z''[t] == -m g - Tanh[z'[t]] 1/2 rcd (z'[t] - Uz)^2,  z == 0,
z' == v Cos[theta],
m x''[t] == -Tanh[x'[t]] 1/2 rcd (x'[t] - Ux)^2,        x == 0,
x' == v Sin[theta] Cos[phi],
m y''[t] == -Tanh[y'[t]] 1/2 rcd (y'[t] - Uy)^2,        y == 0,
y' == v Sin[phi] Sin[theta],
WhenEvent[z[t] == 0, {"StopIntegration"}]}, {x, y, z},
{t, EndTime[theta]}, {v, theta, phi}];

params[target_, {Ux_, Uy_, Uz_}] :=
Quiet@NMinimize[ EuclideanDistance[target, Through[#[#[[1, 1, -1, -1]]]] &@
Through[({x, y, z} /. sol[Ux, Uy, Uz])[v, th, p]]], {v, th, p}]

target = {67, 72, 0};
U = {-7, -3, 0};
params[target, U]
(* {1.29988*10^-7, {v -> 36.2713, th -> 0.642614, p -> 0.779455}} *)


Testing:

Show[
ParametricPlot3D[
Evaluate[Through[#[t]] &@ Through[({x, y, z} /. (sol @@ U))[v, th, p] /. params[target, U][]]],
{t, 0, 10}, PlotRange -> {{0, 80}, {0, 80}, {0, 40}}],
Graphics3D[{PointSize[Large], Point@target}]] Direct hit!

Edit

This is a guide on how to build a model to be able to transfer these calcs to another platform. Please be aware that the full fledge thing requires a lot of computing power because you'll need to solve the ODE a lot of times.

We will solve here for a small bounded rectangle in the target and wind spaces and we will fit (only the velocity, as an example) it with a complete quadric, although fitting larger rectangles may need more sophisticated functions (impossible to know without doing either the calcs or thinking about the physics). Like this:

tab = Table[{{a, b, c, d} = {xt, yt, ux, uy}, params[{xt, yt, 0}, {ux, uy, 0}]},
{xt, 65, 67}, {yt, 70, 72}, {ux, 5, 7}, {uy, 2, 4}];
grid = Partition[Flatten[tab, 4], 2] /. {x_List, y_List} :> {x, Last@y};
vgrid = grid /. {x_List, y_List} :> {Sequence @@ x, v /. y};
cc = Flatten@ Table[coef[i, j, k, l], {i, 0, 2}, {j, 0, 2}, {k, 0, 2}, {l, 0, 2}];
pows = Flatten@ Table[x^i y^j ux^k uy^l, {i, 0, 2}, {j, 0, 2}, {k, 0, 2}, {l, 0,  2}];

nlm = NonlinearModelFit[vgrid, { cc.pows}, cc, {x, y, ux, uy}];
Normal[nlm] Testing:

nlm["FitResiduals"] // Max
(* 0.0590728 *)


So everything seems good. You'll need to calculate some intermediate values to be sure, though.

• @belisarious this is good, but my aim is to generalize and create some estimation functions to predict the landing point and hit a given point at given wind speeds. I've tried to clarify the question a bit. – Julian Feb 26 '15 at 18:03
• @Julian You only need to make a function from my answer. – Dr. belisarius Feb 26 '15 at 18:05
• @Julian Take a look now – Dr. belisarius Feb 26 '15 at 18:42
• better still, but what I am trying to get are some simpler functions, for example polynomial approximations. Would it be possible to map an area with the function you have made and take a power series expansion (in terms of wind speed/SP) in order to approximate it? – Julian Feb 26 '15 at 18:55
• @Julian Pretty simple. For each possible target point you can sample the {Ux,Uy,Uz} space with params[target, U] and then build up an Interpolation[] for each one of {v,theta,phi}. Or use NonLinearModelFit[] if you know or presume some particular functional dependence. – Dr. belisarius Feb 26 '15 at 18:58