% Trajectory of rotating tennis ball
% initialization
global g alpha w etha
g= 9.81; d= 0.063; m= 0.05;  rho= 1.29; 
alpha= pi*d^2/(8*m)*rho;
etha= 1;  
w=20;
%initial condition
h=1; v0=25; theta=pi/180*15;
xin=[0,h,v0*cos(theta),v0*sin(theta)]
% flight time for vacuum
tmaxid=(xin(4)+sqrt(xin(4)^2+2*g*xin(2)))/g
% setting options for the integrator
options= odeset('Events','on','Maxstep',tmaxid/100);
% solution in vacuum
[tid,xid]= ode23('tennisip2',[0 tmaxid],xin,options); 
% solution without spin
[t0,x0]= ode23('tennis0p2',[0 tmaxid],xin,options);  
% solution with spin
[t1,x1]= ode23('tennis1p2',[0 tmaxid],xin,options); 
clf;
axis([0,max(xid(:,1))*1.13,0,max(xid(:,2))*1.13]);
hold on;
plot(xid(:,1),xid(:,2),':r');
plot(x0(:,1),x0(:,2),'--b');
plot(x1(:,1),x1(:,2),'-g');
hold off;
% flight time for the trajectory without spin
tmax0= t0(length(t0))
% flight time for the trajectory with spin
tmax1= t1(length(t1))

