% 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.1 0.2 0.3; ...
    0.0 0.0 0.0 0.0];
%C=[0.0 0.0 0.0 0.0; ...
%    0.0 0.1 0.2 0.3];
% Constants
g = 9.8;
rho = 1.2;  % air density in kg/m^3
A = 0.0375; % cross sectional area in m^2
m = 0.424;  % mass in kg
alpha = rho*A/2/m;
mkr = {'-','--',':','-.'};

% Initial velocity in m/s
v0 = 20;
theta = 50;
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 3],[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','south');
axis image
hold off