plotting - 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] == 0, z'[0] == v Cos[theta], 
    m x''[t] == 0, x[0] == 0, x'[0] == v Sin[theta] Cos[phi], 
    m y''[t] == 0, y[0] == 0, y'[0] == 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] == 0, 
    z'[0] == v Cos[theta], 
    m x''[t] == -Tanh[x'[t]]  1/2 rho Cd A (x'[t] - Ux)^2, x[0] == 0, 

    x'[0] == v Sin[theta] Cos[phi], 
    m y''[t] == -Tanh[y'[t]]  1/2 rho Cd A (y'[t] - Uy)^2, y[0] == 0, 
    y'[0] == 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}], 
 Graphics3D[{Red, Arrowheads[0.04], 
   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.

Answer

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] == 0, 
     z'[0] == v Cos[theta], 
     m x''[t] == -Tanh[x'[t]] 1/2 rcd (x'[t] - Ux)^2,        x[0] == 0, 
     x'[0] == v Sin[theta] Cos[phi], 
     m y''[t] == -Tanh[y'[t]] 1/2 rcd (y'[t] - Uy)^2,        y[0] == 0, 
     y'[0] == 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][[2]]]], 
  {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.