Spadek swobodny z wysokości h z uwzględnieniem siły Coriolisa

(* Wektor położenia (OverVector[r]) punktu materialnego w układzie nieinercjalnym *)

r={x[t],y[t],z[t]}

(* Prędkość (OverVector[v]) punktu materialnego *)

v=D[r,t]

(* Przyspieszenie (OverVector[a]) punktu materialnego *)

acc=D[v,t]

(* Wektor prędkości kątowej (OverVector[ω])  *)

omega={0,w*Cos[ϕ],w*Sin[ϕ]}

(* Wektor przyspieszenia ziemskiego (OverVector[g]) *)

gvec={0,0,-g}; (* g > 0 *)

(* II zasada dynamiki w układzie nieinercjalnym *)

eq1=Simplify[acc==gvec-2*omega\[Cross]v]

(* Warunki początkowe *)

WP={x[0]==0,y[0]==0,z[0]==h,x'[0]==0,y'[0]==z'[0]==0}

(* Dla tego problemu szukam rozwiązań, które uwzględniają już warunki początkowe *)

solWP=FullSimplify[DSolve[{eq1,WP},{x[t],y[t],z[t]},t],w>0&&h>0&&g>0&&Pi/2>ϕ>0]

xs[t_]=Simplify[x[t]/.solWP[[1]]]

ys[t_]=Simplify[y[t]/.solWP[[1]]]

zs[t_]=Simplify[z[t]/.solWP[[1]]]

Simplify[xs[0]]

Simplify[xs'[0]]

Simplify[ys[0]]

Simplify[ys'[0]]

Simplify[zs[0]]

Simplify[zs'[0]]

(* Dane dla swobodnego spadku z wysokości h w Krakowie *)

g=981/100; (* m/s^2 *)

h=50; (* m *)

ϕ=50 Degree;

w=2*Pi/(24*3600); (* 1/s *)

N[w]

t0=Sqrt[2*h/g] (* czas spadku dla przypadku bez siły Coriolisa *)

N[t0]

Plot[xs[t],{t,0,t0},PlotRange->All,AxesLabel->{"t [s]","x [m]"}]

Mamy (bardzo małe) odchylenie w kierunku wschodnim !

Plot[ys[t],{t,0,t0},PlotRange->All,AxesLabel->{"t [s]","y [m]"},PlotPoints->10000]

Plot[{h-(1/2)*g*t^2,zs[t]},{t,0,t0},PlotRange->All,AxesLabel->{"t [s]","z [m]"}]