%% Galileo's Projectile Motion
%
% For the tracker parachute worksheet, you used the software to analyze the trajectory of a falling "grocery store bag" parachute. Now, let's simulate the same scenario and make sure that it matches the Tracker results.
%
% Then, we will modify the design parameters to simulate an actual parachute.
%%
%% Activity: Simulate the parachute's falling motion
%
% Establish relationship between variables & determine equations and functions.
% We developed these in class, and they are repeated below:
% F_g =\textrm{mg}\;\;\;\;\left(1\right)
% F_D =\frac{1}{2}\textrm{𝜌𝐴}C_D {v_y \left(t\right)}^2 \;\;\left(2\right)
% V_y \left(t+\textrm{∆t}\right)=V_y \left(t\right)+\left(F_D +F_g \right)\frac{\textrm{∆t}}{m}\;\;\;\;\left(3\right)
% y\left(t+\textrm{∆t}\right)=y\left(t\right)+V_y \left(t\right)\bullet \textrm{∆t}\;\;\;\;\;\left(4\right)
% In class, you converted these to Matlab code, which you will use below.
% You will use these discretized equations (eqs 1,2 above) to determine the velocity and position in the y direction using the timestep ∆t=dt=0.01 sec.
% The initial conditions that we will use are:
% y(1) = <the measured initial height of your grocery bag parachute>
% Vy(1) = 0 m/sec
%
%% Task 1:
% Similar to the bike share simulation, you will initialize arrays, setup initial values, and then within a for loop, compute new values for the arrays, based on the previous values. For this task, you will compute these values across 40 time steps.
% Notes:
% To help you get started, there is some initial code below.
% Use Nsteps=40 to set the number of time steps for the computation. This variable is also used to initialize the arrays (see below). This makes it easy to change this parameter and it is good programming practice!
% Remember to include a ";" at the end of EVERY assignment statement. Otherwise, Matlab will display assignment results for all 40 times through the loop, which will make the output cluttered and slow down your program execution!
% Also, think about the parameters to use in the for statement. If you use for i=1:Nsteps and Nsteps=40, then you will execute the loop 40 times. Including the initial conditions, this is 41 time steps, which is one extra step than I've asked you do to. So the for statement should be for i=1:Nsteps-1.
%

clc; clear; %Clean the workspace
close all; %Close all the existing figures 

%initialize arrays
Nsteps = 40;
y=zeros(Nsteps,1);
Vy=zeros(Nsteps,1);

% Next, you will setup all initial conditions. This should include dt, as well as all of the parameters needed for your calculations
% density = 1.225 kg/m^3 (note that this is in kg -- make sure that all of your weight mesurements are consistently in kg or g). This is assuming 1 atm pressure and 15˚C temp, as shown here: https://www.engineeringtoolbox.com/air-density-specific-weight-d_600.html. Note that in reality the density will change as the parachute falls and pressure and temperature change, so we are making an approximation by using a fixed value and this is a potential source of error.
% C=1.5. This is slightly larger than the drag coefficient for an "hollow semi-sphere", and we're assuming that our parachute will have slightly higher drag. One good source for drag coefficients of different shapes is below, but note that drag coefficients can vary widely depending on the parachute shape, so this is a potential source of error. https://www.engineeringtoolbox.com/drag-coefficient-d_627.html
% A = area of parachute computed from your diameter measurement
% m = mass of parachute that you measured in class
% g = -9.8 m/sec^2 (note that this is negative because parachute is falling in the -y direction)
% y(1) is the initial height of your parachute that you measured
% Vy(1) is the initial velocity of your parachute just before you dropped it (I'll leave it to you to determine what is the appropriate value of this!)
% dt=0.01 sec

%setup initial conditions for all variables


% Here is your for loop where you will compute y and Vy at each time step

%for loop
for i=1:Nsteps-1    %see note above for why this is Nsteps-1
    %compute values for Fd, y(i+1), Vy(i+1) using the
    %equations that were given for these parameters

end

% Display the simulation values for y and Vy. What do you expect for these parameters?  (Hint: look at your Tracker plots). Are you getting what you would expect?

disp(y)
disp(Vy)

%
%% Task 2:
% For this task, you will modify the code from task 1 to display the simulation results as a plot instead of a table.
% To make a plot, you need data for the horizontal and vertical axes. The horizontal axis is time, so you need to generate a time array. The time array is 0, dt, dt*2, dt*3, etc.
% There are several ways to generate this: using the colon operator, using the linspace function, or by initializing t using the zeros function, and computing t(i+1)= t(i)+dt within the for loop. Any of these methods are fine! If you are not sure how to do this, then look it up!
% Here are the steps for this task:
% Copy your code from task 1 and paste it below, but delete the two lines of code at the end that display the results of y and Vy
% Add code to generate values for the t array using one of the methods described above
% Plot y(t) vs t, Vy(t) vs t
% Put relevant titles for the plots and add x-axis and y-axis labels.

% Your code begins here:


%
%% Task 3:
% Note that after 40 time steps, the parachute never reaches the ground, but we want to keep running the simulation until the parachute hits the ground. We don't know how many time steps this will take, so a for loop is not really appropriate here. Instead let's use a while loop and keep repeating it as long as y(i) >= 0
% Copy the code from task 2 and paste it below
% Replace the for loop with a while loop. The condition of the while loop should repeat the loop as long as y(i) >= 0.
% Remember that with a while loop, you'll need to add i=1;  before the loop, and i=i+1; right before the end statement of the loop
% When initializing the arrays with the zeros function, use a length of Nsteps=1000; which should be large enough so that you capture the entire trajectory of the parachute reaching the ground
% Plot y and Vy as above, and your figures should have titles and you should label the axes. Since you only executed the while loop i times, your plot statement should only plot the first i elements. For example: plot(t(1:i), y(1:i))

% Your code begins here: 




%
%% Task 4: Validate the simulation by comparing the simulation results to experimental results from Tracker.
% Note that you should NOT move on to task 5 until you have validated your data here. You will receive a lower grade if you do not properly validate your data.
% Step 1: Copy your code from Task 3, but remove the plotting commands at the end
% Step 2: Export your tracker data to a .csv file. Look at the video linked from the W4 instructions for how to do this. Move the .csv file to the same folder as this .mlx file, and set your Matlab current folder to thsi folder.
% Step 3: Now, you can load the .csv file into Matlab as a single variable. The Matlab commands to accomplish this are below.
% Step 4: Compare the experimental data (from Tracker) to the simulated data generated by Matlab. Plot y vs t and y_measured vs t_measured on one graph.  Do the same for Vy. Your figures should have a legend that identifies the two curves on each plot as simulation and experimental results. The measured data should be plotted with 'o' markers to show each measured data point, and the experimental data should be plotted as a line. And as in the above tasks, your figures should have titles and you should label the x and y axes.
% Note: For each figure, you should plot the experimental data, followed by "hold on", and then plot the simulation data followed by "hold off". This will plot both graphs on the same axes.
% Step 5: Does the experimental data match the simulated data? It does not need to be a perfect match between simulation and measurement results, but it should be reasonably close. If the answer is no, then try to figure out the source of the error before you move on.

% Copy the code from task 3, but leave off the 
%  plotting commands


% Import the Tracker data
data=readtable('trackerdata.csv');
t_measured=data.Var1;
y_measured=data.Var2;
Vy_measured=data.Var3;

%Plot the simulation data for y vs t and label axes

%hold on preserves this plot so that the next plot 
% goes on the same set of axes
hold on

%Plot the measured Tracker data for y_measured vs t_measured. 
% Be sure to include a legend to distinguish the 
% simulation from the experiment


%hold off ends this plot
hold off

%Repeat the above for Vy vs t and Vy_measured vs t_measured


%
%
%% Task 5: Now that you have verified your simulation, let's try this for a real parachute
%
% Step 1: Copy your code from task 4, but eliminate the last section where you plot measured vs simulation data
% Now that you have validated your simulation using the shopping bag parachute, here are the new parameters to use for your skydiver:
% d (m): diameter of an actual parachute (do some research online to find out a good starting point for this value -- you will want to adjust this to get an appropriate terminal velocity)
% m (kg): use the mass of your body (in kg) and add 25 lb (11.3 kg) to account for the weight of the parachute and other gear; note that you do not need to reveal your own weight, and you can choose another person's weight here if you wish!
% y(1) (m): assume that you are opening the parachute at 1000 ft, or 305 m
% Vy(1) (m/sec): initial y velocity is zero
% Next steps
% Plot y vs t and Vy vs t. What is the terminal velocity?
% Do research to identify a "reasonable" velocity for a safe landing (save this citation for your report). Did your simulation result in a reasonable value for a safe landing?
% Keep adjusting your value for d (parachute diameter) until you get a "reasonable" value without making the parachute so large in diameter so that it would be too heavy or cumbersome.
%

% Copy your task 4 code here and modify as described above
% You only need to plot y vs t and Vy vs t
% 