%% This script calculates the fall of a baseball with and without drag.
% It uses Euler's method: ValueNext = ValueNow + RateOfChangeNow*deltat


%% 1 Housekeeping
close all;clc;  %Clean up previous work, start "fresh"

%% 2 Define Constants.  

%End the lines  with ";" to supress echo

Cd = 0.35;      % Drag coefficient (dimensionless)
A = 4.3e-3;     % Cross-sectional area of projectile (m^2)
g = 9.8;        % Gravitational acceleration (m/s^2)
m = 0.145;      % Mass of baseball (kg)
density = 1.2;  %Air density in kg/m^3

C =0.5*Cd*density*A; %Drag Coefficient C in Fdrag = C v^2

% To input the data at run time,  we use name = input('Message To User'), 
% for example:  v0_y = input('Enter initial y velocity (m/s): '); 

% Here we just give all values in the script.
v0 = 100;
alpha0 = 35; %In degrees!
v0_x = v0*cosd(alpha0);
v0_y = v0*sind(alpha0);
x0 = 0;
y0 = 0.7;

%% 3 Define Calculation parameters (number of steps, time step, etc.)
N = 400; 
dt = 0.05; 


%% 4 Preallocate memory for the quantities we want to track at every step.

%Usually they are Nx1 matrices (N rows, 1 column).
t   = zeros(N,1); 
x   = zeros(N,1); 
v_x = zeros(N,1);
Fnet_x = zeros(N,1);
xideal = zeros(N,1);

y   = zeros(N,1); 
v_y = zeros(N,1);
Fnet_y = zeros(N,1);
yideal = zeros(N,1);
 

%% 5 Define Update Loop Starting Conditions (Values on Row 1)
t(1) = 0.0;
x(1) = x0;
y(1) = y0;
v_x(1) = v0_x;
v_y(1) = v0_y;
Fnet_x(1) = 0;
Fnet_y(1) = -m*g;
xideal(1)= x0;
yideal(1)= y0;

%% 6 Run the Update Loop for N steps, for a maximum simulation time of N*dt
for step=1:N-1 % loop over the timesteps

%Careful how you write Fy: y is positive up (-mg) and Fdrag is in the -v
%direction  : Fx = -C abs(V)*Vx, Fy = -C abs(V)*Vy    

Fnet_x(step) = -C*sqrt((v_x(step))^2 + (v_y(step))^2)*v_x(step) ;
Fnet_y(step) = -m*g -C*sqrt((v_x(step))^2 + (v_y(step))^2)*v_y(step) ;
% MP: Momentum principle: pnext = pnow + F*dt, divide by m. 
v_x(step+1) = v_x(step) + (Fnet_x(step)/m)*dt  ;
v_y(step+1) = v_y(step) + (Fnet_y(step)/m)*dt  ;

%PU: Position Update
x(step+1) = x(step) + v_x(step)*dt;
y(step+1) = y(step) + v_y(step)*dt;

%Time Update
t(step+1) = t(step) + dt;

%The ideal case has a known solution
xideal(step+1) = x0+ v0_x*t(step);
yideal(step+1) = y0+ v0_y*t(step) -0.5*g*(t(step))^2;

%Do not let the ball go "underground"
if yideal(step+1)<0
    yideal(step+1)=0;
end

if y(step+1)<0
    y(step+1)=0;
end

end

%% 7 Plot the solutions
plot(t,y,'.r'); 
%keep the previously plotted lines when plotting the next one
hold on;
%plots the ideal (no drag) solution 
plot(t,yideal,'.b');

legend('With Drag','Without Drag');
xlabel('Time (s)');  ylabel('Height (m)');
title('Fall with and without Drag');
hold on;
%release the hold
hold off;