Skip to main content

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):


enter image description here



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)$$


enter image description here


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}]]

Mathematica graphics


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]

Mathematica graphics


Testing:


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

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


Comments

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...