% Trajectory of spinning tennis ball
% initialization
global g alpha w etha
% basic constants im MKS units
g = 9.81; d = 0.063; m = 0.05; rho = 1.29; 
alpha=pi*d^2/(8*m)*rho;
etha = 1;
w = 20;
% initial conditions
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; 
% solution in vacuum
[tid, xid] = ode23('tennisip', [0 tmaxid], xin);  
% solution without spin
[t0, x0] = ode23('tennis0p', [0 tmaxid], xin);  
% solution with spin
[t1, x1] = ode23('tennis1p', [0 tmaxid], xin);  
N = max(xid(:, 1)); x = 0:N/100:N;
axis([0,max(xid(:,1)), 0, max(xid(:,2))])
hold on;
plot(x, spline(xid(:,1), xid(:, 2), x), ':r');
plot(x, spline(x0(:,1), x0(:, 2), x), '--b');
plot(x, spline(x1(:,1), x1(:, 2), x), '-g'); 
hold off;

% Determine flight time  by inverse interpolation 
[z, j] = min(x0(:, 2)); [y, i] = max(x0(:, 2));
tmax0 = spline(x0(i + 1:j, 2), t0(i + 1:j), 0)
[z, j] = min(x1(:, 2)); [y, i] = max(x1(:, 2));
tmax1 = spline(x1(i+1:j,2), t1(i+1:j), 0)

