How to calculate position from 3 dimensional acceleration data?
Assuming that all the initial positions (x[0]
, y[0]
and z[0]
) and the initial velocities (x'[0]
, y'[0]
and z'[0]
) are equal to 0
you can do:
adat = Rest@Import["http://pastebin.com/raw.php?i=jZ57mqZT"];
{ax, ay, az} = Interpolation /@ (adat[[All, {1, #}]] & /@ {2, 3, 4});
{xt, yt, zt} = (x /. Quiet@First@NDSolve[{
x[0] == 0, x'[0] == 0,
x''[t] == #[t]
}, x, {t, First@adat[[All, 1]], Last@adat[[All, 1]]}] & /@ {ax, ay, az});
Plot[{ax[t], ay[t], az[t]}, {t, First@adat[[All, 1]],Last@adat[[All, 1]]},
PlotLabel -> "Acceleration", Frame -> True,
FrameLabel -> {"Time", "Aceleration"}]
Plot[{xt[t], yt[t], zt[t]}, {t, First@adat[[All, 1]], Last@adat[[All, 1]]},
PlotLabel -> "Positions", Frame -> True,
FrameLabel -> {"Time", "Position"}]
ParametricPlot3D[{xt[t], yt[t], zt[t]}, {t, First@adat[[All, 1]], Last@adat[[All, 1]]},
PlotRange -> {{0, -2}, {-1, 1}, {0, -2}},
AxesLabel -> {"X", "Y", "Z"}]
A different interpretation would be that the acceleration is always in the same direction, but that the sensor actually rotates. In that case:
{pt, px, py, pz} = Transpose@Rest@Import["http://pastebin.com/raw.php?i=jZ57mqZT"];
{pa, pθ, pϕ} = Transpose@CoordinateTransform["Cartesian" -> "Spherical", Transpose@{px, py, pz}];
θi = Interpolation[Transpose@{pt, pθ}];
ϕi = Interpolation[Transpose@{pt, pϕ}];
ai = Interpolation[Transpose@{pt, pa}];
Using NDSolve
and assuming that the acceleration data is the acceleration observed in an object (as opposed to measured within the object)
hi = (h /.
First@Quiet@
NDSolve[{h[0] == 0, h'[0] == 0, h''[t] == -ai[t]},
h, {t, Min@pt, Max@pt}]);
So we can plot the elevation and azimuthal angle
ListLinePlot[{pθ, pϕ}, PlotRange -> {-2 π, 2 π}, FrameTicks ->{Automatic, π Range[-2, 2]}, Frame -> True, PlotLabel -> "Orientation", FrameLabel -> {"Time", "Angles"}]
The acceleration
Plot[hi[t], {t, Min@pt, Max@pt}, Frame -> True, PlotLabel -> "Trajectory", FrameLabel -> {"Time", "Hight"}]
And the trajectory of an object accelerating as in the data
Plot[hi[t], {t, Min@pt, Max@pt}, Frame -> True, PlotLabel -> "Trajectory", FrameLabel -> {"Time", "Hight"}]
To visualize the orientation over time
Animate[
Graphics3D[
{
Opacity[1]
, Red
, Arrow[{{0, 0, 0}, {1, 0, 0}}]
, Green
, Arrow[{{0, 0, 0}, {0, 1, 0}}]
, Blue
, Arrow[{{0, 0, 0}, {0, 0, 1}}]
, Opacity[0.2]
, GeometricTransformation[
Cuboid[{-(1/2), -(1/2), -(1/2)}, {1/2, 1/2, 1/2}]
, Composition @@ {RotationTransform[ϕi[t], {0, 0, 1}],
RotationTransform[θi[t], {0, 1, 0}]}
]
}, PlotRange -> {{-2, 2}, {-2, 2}, {-2, 2}}], {t, Min@pt, Max@pt}]
Now that I figured out how to do Composition
transformations just for eye candy we could mix the trajectory with the rotation angle.
WARNING I'm mixing pears and oranges here, this is misleading and the reality here depends on the details on how and where the acceleration was measured. The data can NOT come from an accelerometer on a falling cube with initial velocity zero (as in the animation) as the acceleration measured in free fall would be minimal until terminal velocity is reached. If the acceleration if measured from 'outside' then we could not measure the angles.
The angles were estimated assuming that gravity acceleration always points down, and the observed magnitude (Norm
) was fairly constant and close to $\vec{g}$ assuming units of $m/s^2$. If this comes from an accelerometer, like the ones in a phone, implies that it was not falling, but just hanging.
Therefore, this animation is not a physically accurate picture, just fun with Mathematica.
Animate[
Grid[{{
Graphics3D[
{
Opacity[1]
, Red
, Arrow[{{0, 0, 0}, {1, 0, 0}}]
, Green
, Arrow[{{0, 0, 0}, {0, 1, 0}}]
, Blue
, Arrow[{{0, 0, 0}, {0, 0, 1}}]
, Opacity[0.2]
, GeometricTransformation[
Cuboid[{-(1/2), -(1/2), -(1/2)}, {1/2, 1/2, +(1/2)}]
, Composition @@ {
TranslationTransform[{0, 0, hi[t]}]
, RotationTransform[ϕi[t], {0, 0, 1}]
, RotationTransform[θi[t], {0, 1, 0}]
}]
}
, PlotRange -> {{-2, 2}, {-2, 2}, {-5, 2}}
, ImageSize -> 300]
,
, Plot[hi[tt], {tt, Min@pt, Max@pt}, Frame -> True,
PlotLabel -> "Trajectory", FrameLabel -> {"Time", "Hight"},
ImageSize -> 600,
Epilog -> {PointSize[Medium], Point[{t, hi[t]}]}]
}}], {t, Table[ti[k], {k, 1, 24, 0.1}]}]