I think the only real problem keeping you from using the more modern `Quantities`

framework is the fact that `Integrate(a0, t, t)`

and `a0 Integrate(1, t, t)`

are not equivalent if one or more accelerations in `a0`

are zero. In the first case we get `{0 m/s^2, -4.9 t^2 m/s^2}`

and in the second case we get `{t^2 (0 m/s^2), t^2 (-4.9 m/s^2)}`

. Since we’re integrating constants, we can pull the constants out of the integral and put them in front. I’ve done it for both `v0`

and `a0`

, though it was only absolutely necessary for `a0`

. Unfortunately, because of the way the replacement works, Mathematica will print a warning that m/s and m are incompatible units, *then* the substitution happens and it realizes the units are fine.

```
vectorDirection(vec_) := ArcTan(vec((2))/vec((1)))*180/(Pi)
r0 = Quantity({0, 0}, "Meters");
s0 = Quantity(37, "Meters"/"Seconds");
v0 = AngleVector(s0, 53.1 Degree)
a0 = Quantity({0, -9.8}, "Meters"/"Seconds"^2);
r0 + v0 Integrate(1, t) + a0 Integrate(1, t, t) /.
t -> Quantity(2, "Seconds")
{Norm(
v0 + a0 Integrate(1, t) /. t -> Quantity(2, "Seconds")),
vectorDirection(
v0 + a0 Integrate(1, t) /. t -> Quantity(2, "Seconds"))}
r1 = {x, y};
v1 = {v1x, Quantity(0, ("Meters")/("Seconds"))};
Solve({r1 == r0 + v0 Integrate(1, t) + a0 Integrate(1, t, t),
v1 == v0 + a0 Integrate(1, t)}, {x, y, t, v1x})
r1 = {x, Quantity(0, "Meters")};
v1 = {v1x, v1y};
Last@Solve({r1 == r0 + v0 Integrate(1, t) + a0 Integrate(1, t, t),
v1 == v0 + a0 Integrate(1, t)}, {x, t, v1x, v1y})
```

Of course all the quantities and integrals can be entered in the more traditional format when using Mathematica, but it doesn’t translate well to StackExchange. I do think it’s a bit strange to carry around all these integrals and repeatedly solve the same equations over and over. My usual process would be to solve the integral once, if possible, and then be done with it. Perhaps you have a reason you want to keep the integrals, but I think it’s cleaner without. With the initial values defined as above:

```
pos = Integrate(a, t, t, GeneratedParameters -> C) /. {C(1) -> r0, C(2) -> v0, a -> a0}
pos /. t -> Quantity(2, "Seconds")
{Norm(#), vectorDirection(#)} &@(D(pos, t) /.
t -> Quantity(2, "Seconds"))
Solve({{x, y} == pos, {v1x, Quantity(0, ("Meters")/("Seconds"))} ==
D(pos, t)}, {x, y, t, v1x})
Last@Solve({{x, Quantity(0, "Meters")} == pos, {v1x, v1y} ==
D(pos, t)}, {x, t, v1x, v1y})
```

If we were using the derivative of the position a lot, we could also just define a new equation like `vel = D(pos, t)`

. I also switched the definition of position here. From a physics perspective, I don’t think the original definition is quite right. It gets the same answer, but a double indefinite integral should have two unknowns (initial velocity and initial position, here). When using the double integral, we’re really solving the differential equation `r''(t) == a`

and integrating both sides twice to get `r(t) = 1/2 a t^2 + v0 t + r0`

, so we don’t need the definition to contain 3 separate parts.