% Projectile motion with damping and lift problem 27
% 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

theta = -.2692

% [Drag; Lift] coefficients
% Parts a b c
C=[ 0.0,  0.0   0.0; ...
    0.0, -0.031 0.031];
C=[ 0.55,  0.55   0.55; ...
    0.0, -0.031 0.031];




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

% Initial velocity in m/s
v0 = 50;
%theta = -.2740765786; %in radians
v0x = v0*cos(theta);
v0z = v0*sin(theta);

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

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

figure(1)
clf
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 .18],[x0,z0,v0x,v0z]);
    
    % Plot solution z vs x
    plot(Y(:,1),Y(:,2),'k','linestyle',mkr{i})
end    

x=[6.4,6.4];
y=[0,0.914];
plot(x,y,'linewidth',2)
xlabel('x (m)')
ylabel('z (m)')

legend(mylegend,'location','northeast');
axis image
axis([6.3,6.65,.85,.95])
hold off
