restart; # motion of ball in vacuum eqnid := diff(x(t), t$2) = 0, diff(z(t),t$2) = -g: varid := {x(t), z(t)}: initcid := x(0) = 0, z(0) = h, D(x)(0) = v0*cos(theta), D(z)(0) = v0*sin(theta): #solution for ball in vacuum resid := dsolve({eqnid, initcid}, varid); # numeric solution of ball in air # Basic constants in S g := 9.81: # gravitational acceleration: d := 0.063: m := 0.05: # diameter and mass of ball rov := 1.29: # density of air alpha := evalf(Pi*d^2/(8*m)*rov): v := (dx^2 + dz^2)^(1/2): # definition of velocity Cd := .508+1/(22.503 + 4.196/(w/v(t))^(5/2))^(2/5): Cm := eta/(2.202 + .981/(w/v(t))): # eta =+-1 defines direction of rotation, # for top spinning eta = 1 var := {x(t), z(t), dx(t), dz(t)}: # initial conditions initc := x(0) = 0, z(0) = h, dx(0) = v0*cos(theta), dz(0) = v0*sin(theta): # equations of motion of ball in air # rotation (drag force only) eqnt0:= diff(x(t),t) = dx(t), diff(dx(t),t)= -0.508*alpha*dx(t)*v(t), diff(z(t),t) = dz(t), diff(dz(t),t)= -g-0.508*alpha*dz(t)*v(t): # equations of motion of rotating ball in air # (influence of Magnus effect) eqnt1:= diff(x(t),t) = dx(t), diff(dx(t),t)= (-Cd*dx(t) + Cm*dz(t))*alpha*v(t), diff(z(t),t) = dz(t), diff(dz(t),t)= -g-(Cd*dz(t) + Cm*dx(t))*alpha*v(t): #numeric values of initial parameters h := 1: v0 := 25: theta := Pi/180*15: # theta= 15 degrees w := 20: eta := 1: # solution for non rotating ball res0 := dsolve({eqnt0,initc},var,numeric); # result is a set of equations res0(0.5); # solution for rotating ball res1 := dsolve({eqnt1, initc}, var, numeric); # calculation of tmax - time of falling to earth # for ball in vacuum tmaxid := fsolve(subs(resid, z(t)) = 0, t, t = 0..5); read `zzero.map`: #!!!S.B. Not present in the book # calculation of the flight time for the other models tmax0:= zzero(res0, tmaxid, z(t), dz(t)); tmax1:= zzero(res1, tmaxid, z(t), dz(t)); # making graphs: Gid := plot([subs(resid, x(t)), subs(resid, z(t)), t = 0..tmaxid], linestyle = 2): # for models with numeric solution # calculation of tables [x(t),z(t)] for plotting tablepar := proc(u, x, y, xmin, xmax, npoints) local i,Step; Step := (xmax - xmin)/npoints; [seq([subs(u(xmin + i*Step), x), subs(u(xmin + i*Step) ,y)], i = 0 .. npoints)] end: G0 := plot(tablepar(res0, x(t), z(t), 0, tmax0, 15), linestyle = 3): G1 := plot(tablepar(res1, x(t), z(t), 0, tmax1, 15), linestyle = 1): # plotting of all graphs plots[display]({Gid, G0,G1});