function dr = motion(t,r )
global g c;
vv=sqrt(r(3)*r(3)+r(4)*r(4));
dr=[r(3);r(4);0.-c*vv*r(3);-g-c*vv*r(4)];



