%% 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;