Rzut poziomy w jednorodnym ziemskim polu grawitacyjnym

Ruch punktu materialnego dany jest układem równań
x = (Subscript[v, 0])*t,
y = h - (1/2)*g*(t^2),
z=0,
gdzie (Subscript[v, 0]), h i g to dodatnie stałe.
Znaleźć równanie toru, prędkość, przyspieszenie, szarpnięcie, składową styczną przyspieszenia (Subscript[a, st]), promień krzywizny ρ, składową normalną przyspieszenia  (Subscript[a, nor]) oraz drogę (Subscript[s, 12]) , jaką punkt materialny przebędzie od chwili (Subscript[t, 1]) do chwili (Subscript[t, 2]), gdzie 0 < (Subscript[t, 1]) < (Subscript[t, 2]). 
Przykład:

v0=10; (* m/s *)

h=200; (* m *)

g=10; (* m/s^2  *)

tmax=Sqrt[2*h/g];

N[tmax]

(* W tym miejscu, używając ParametricPlot,  pomijamy składową z *)

p=ParametricPlot[{v0*t,h-(1/2)*g*t^2},{t,0,tmax},PlotRange->{{-3,v0*tmax+5},{-5,200+5}},AspectRatio-> 1/2, AxesLabel->{"x [m]","y [m]"},PlotLabel-> "(Subscript[v, 0])=10 m/s, h=200 m, g=10 m/\!\(TraditionalForm\`\(\*SuperscriptBox[\(s\), \(2\)]\(\\ \)\)\)"]

(*  Export[ NotebookDirectory[] <> "rzut_poziomy.png",p]; *)

Animate[Show[p,ListPlot[{{v0*t,h-(1/2)*g*t^2}}, PlotStyle->{PointSize[0.02],Orange}]],{t,0,tmax}]

ClearAll[v0,h,g];

(* wektor położenia *)
r={v0*t,h-(1/2)*g*t^2,0};

(* wektor prędkości *)
v=D[r,t]

(* wektor przyspieszenia *)
a=D[v,t]

(* wektor szarpnięcia *)
sz=D[a,t]

(* iloczyn wektorowy prędkości i przyspieszenia *)
acrossv = Cross[a,v]

FullSimplify[acrossv,v0>0&&h>0&&g>0]

absacrossv=FullSimplify[Sqrt[acrossv.acrossv],v0>0&&h>0&&g>0]

(* długość wektora prędkości, czyli tzw. szybkość *)
absv=FullSimplify[Sqrt[v.v],v0>0&&h>0&&g>0]

(* składowa styczna przyspieszenia *)
ast=FullSimplify[D[absv,t],v0>0&&h>0&&g>0]

(* promień krzywizny *)
rho=FullSimplify[absv^3/absacrossv,v0>0&&h>0&&g>0]

(* składowa normalna przyspieszenia *)
anor=FullSimplify[absv^2/rho,v0>0&&h>0&&g>0]

(* sprawdzenie, że ast i anor są rzeczywiście składowymi wektora przyspieszenia *)
chkabsa=FullSimplify[ast^2+anor^2==a.a,v0>0&&h>0&&g>0]

Limit[anor,t->Infinity]

Simplify[Limit[ast,t->Infinity],g>0]

(* Droga przebyta od chwili t1 do chwili t2 dana jest wzorem 
(∫_t1^t2)|(OverVector[v])(t)|\[DifferentialD]t . Jest to całka oznaczona z długości wektora prędkości w przedziale (t1,t2), więc w zasadzie możemy użyć komendy 

s12= FullSimplify[Integrate[absv,{t,t1,t2}],v0>0&&h>0&&g>0&&t1≥0&&t2≥0&&t2≥t1]

 UWAGA: znacznie szybciej jest jednak policzyć najpierw całkę nieoznaczoną, a potem z niej całkę oznaczoną jako różnicę wartości całki nieoznaczonej w t2 i t1 *)

absv

(* Liczenie całki nieoznaczonej jest operacją odwrotną do liczenia pochodnej ... *)

intabsv=FullSimplify[Integrate[absv,t],v0>0&&h>0&&g>0]

(* ... i łatwo to sprawdzić: *)

Simplify[D[intabsv,t]==absv]

s12= FullSimplify[(intabsv/.{t->t2}) - (intabsv/.{t->t1}), v0>0&&h>0&&g>0]

(* W szczególności dla naszego przykładu liczbowego mielibyśmy dla t1=0 i t2=tmax: *)

s12p=s12/.{h->200,v0->10,g->10,t1->0,t2->tmax}

N[s12p]

(* To jest długość w metrach części paraboli widocznej na rysunku i animacji *)