% Projectile motion with damping and lift
% We solve
%    ax = -alpha*v*(CD*vx+CL*vz)
%    az = alpha*v*(-CD*vz+CL*vx)-g
% for x,z,vx,vz given ICs x0, v0, v0, theta
clear
global g CD CL alpha

% [Drag; Lift] coefficients
C=[0.0, 0.42 ; ...
    0.0, 0.0];

% Constants
g = 9.8;
rho = 1.21;  % air density in kg/m^3
D = 0.0427;  % Diameter
A = pi*D^2/4; % cross sectional area in m^2
m = 0.046;  % mass in kg 
nu =1.47e-5; % viscosity
alpha = rho*A/2/m;
mkr = {'-','--',':','-.'};

% Initial velocity in m/s
v0 = 36;
theta = 30; %in degrees
v0x = v0*cos(theta*pi/180);
v0z = v0*sin(theta*pi/180);

% Range and height with no drag or lift
range = v0^2*sin(2*theta*pi/180)/g;
height = v0^2*sin(theta*pi/180)^2/g/2;

% Initial position in m
x0 = 0;
z0 = 0;

figure(1)
hold on
mylegend='';
for i=1:length(C(1,:))
    CD=C(1,i);
    CL=C(2,i);
    mylegend=strvcat(mylegend,['CD = ',num2str(CD), ' CL =', num2str(CL)]);
        % Solve differential equations
    [T,Y] = ode45('ballf',[0 4],[x0,z0,v0x,v0z]);
    
    % Plot solution z vs x
    plot(Y(:,1),Y(:,2),'k','linestyle',mkr{i})
end    
%title('Projectile Motion')
xlabel('x (m)')
ylabel('z (m)')
%scale = max([range,height]);
%axis([0,scale,0,scale])
legend(mylegend,'location','northeast');
axis image
axis([0,range,0,height])
hold off

figure(2)
v=Y(:,3).^2+Y(:,4).^2;
R=D*v/nu;
plot(T,R)
xlabel('t (s)')
ylabel('Re ')
title('Reynolds Number vs Time')

figure(3)
% Plot hieght vs time. 

