most recent

%% ENSC180-Assignment2

% Student Name 1: Jose Lopez

% Student 1 #: 301317479

% Student 1 userid (email): jlopezsa@sfu.ca

% Student Name 2: Cyrus WaChong

% Student 2 #: 301306459

% Student 2 userid (email): cwachong@sfu.ca

% Below, edit to list any people who helped you with the assignment,

% Helpers: Selikem Kwaznadovia, Sterling Smith, Yutaka Yen, Tal Kazakov, Rameshwar
%  Kannan, Santhosh Nandakumar, Sirpreet Kaur Dhillon

%% Instructions:
% * Put your name(s), student number(s), userid(s) in the above section.
% * Edit the "Helpers" line.
% * Your group name should be "A2_<userid1>_<userid2>" (eg. A2_stu1_stu2)
% * Form a group
%   as described at:  https://courses.cs.sfu.ca/docs/students
% * Replace "% <place your work here>" below, or similar, with your own answers and work.
% * You can copy your work from your other functions and (live) scripts and as needed.
% * Nagvigate to the "PUBLISH" tab (located on top of the editor)
%   * Choose pdf as "Output file format" under "Edit Publishing Options..."
%   * Click "Publish" button. Ensure a report is automatically generated
% * You will submit THIS file (assignment2.m),
%   and the PDF report (assignment2.pdf).
% Craig Scratchley, Spring 2017

%% main

function main

clf

% constants -- you can put constants for the program here
%MY_CONST = 123;

% variables -- you can put variables for the program here
%myVar = 456;

% prepare the data
% <place your work here>
X=xlsread('data_clean_more_fixed_simplest.xlsx', 'data_clean_more_fixed_label', 'B1:B15348');
X=X(isfinite(X(:, 1)), :);
V=xlsread('data_clean_more_fixed_simplest.xlsx', 'data_clean_more_fixed_label', 'C1:C15348');
V=V(isfinite(V(:, 1)), :);
T=xlsread('data_clean_more_fixed_simplest.xlsx', 'data_clean_more_fixed_label', 'A1:A15348');
T=T(isfinite(T(:, 1)), :);

% myVector(isnan(myVector))=[];

% <put here any conversions that are necessary>

%% Part 1
% Answer some questions here in these comments...
% How accurate is the model for the first portion of the minute?
% They appear to be almost identical, therefore rendering the accuracy
% quite high.

% How accurate is the model for the last portion of that first minute?
% it diminishes as time continues; accuracy is less than the first portion.

% Comment on the acceleration calculated from the measured data.
% Is there any way to smooth the acceleration calculated from the data?
% We could use the smooth function, which we found through Sirpreet, or we
% could take more info from an increment in time, which would smooth the
% graph.

part = 1;
endtime=60
starttime=0
%[t1,x1] = ode45(@fall, [starttime,endtime], [38969.4, 0]);
%plotter(X,V,T,'Part1 - Freefall')
% <call here your function to create your plots>


%% Part 2
% Answer some questions here in these comments...
% Estimate your uncertainty in the mass that you have chosen (at the
%     beginning of the jump).
% between 1- 3% of the total mass, due to small things, like food in his
% sytems, liquids, equipment on him, and many other small differences.

% How sensitive is the velocity and altitude reached after 60 seconds to
%    changes in the chosen mass?
% so after a bit of testing, small changes such as 2 to 3 kilos, change
% your altitude from the surface by around 5x10^3 m  in the 60 second interval,
% and around an increase in 8-11 m/s at 60 seconds after the jump.

part = 2;
endtime=60
starttime=0
%[t1,x1] = ode45(@fall, [starttime,endtime], [38969.4, 0]);

% <call here your function to create your plots>
%plotter(X,V,T,'Part2 - Simple Air Resistance' )

%% Part 3
% Answer some questions here in these comments...
% Felix was wearing a pressure suit and carrying oxygen. Why?
%     What can we say about the density of air in the stratosphere?
%     How is the density of air different at around 39,000 meters than it
%     is on the ground?
%
% The density of the air in the stratosphere is less than the surface of
% the earth, meaning there is less molecules per specified area, and
% without the tank of oxygen and suit, he would die from lack of oxygen,
% dangerous radiation, and the cold. The density is about 272 times
% stronger on the crust compared to 39.000 feet up, on earth, stdatmo
% returns a value of 1.225, while up at 39.000, it returns 0.0045, and
% looking at standard graphs for air pressure at atmospheres ( I’m aware
% they are not accurate due to all the factors that affect it) is ~250 times
% also.

% What are the factors involved in calculating the density of air?
%     How do those factors change when we end up at the ground but start
%     at the stratosphere?  Please explain how calculating air density up
%     to the stratosphere is more complicated than say just in the troposphere.

% There are quite a few factors, like temperature, air pressure, altitude and
% weather (and humidity). The air pressure will increase significantly, the
% weather would change quite a bit, altitude would obviously drop, and
% temperature would decrease, due to the ozone layer absorbing the UV rays creating heat energy.
% and leading to a higher overall density of air. The troposphere has the
% most density, due to all layers on top pushing on it, squeezing molecules
% closer together.


% What method(s) can we employ to estimate [the ACd] product?
% We can take the drag force and divide by air density (rho), divide by
% velocity squared, and all multiplied by 2. This is the exact way to get
% it. We can also use newton's second law to find the ACd.

% What is your estimated [ACd] product?
% 0.12
%
% [Given what we are told in the textbook about the simple drag constant, b,]
%   does the estimate for ACd seem reasonable?
% <put your answer here in these comments>

part = 3;
endtime=270
starttime=0
% <place your work here>
%[t1,x1] = ode45(@fall, [starttime,endtime], [38969.4, 0]);
%plotter(X,V,T,'Part3 - Simple Air Resistance' )

%% Part 4
% Answer some questions here in these comments...
% What is the actual gravitational field strength around 39,000 meters?
%   (See Tipler Volume 1 6e page 369.)

% Using the gravitational Field Strength Formula, we can calculate that
% (6.67*10^-11)*(5.972*10^24)/(6371000+39000)^2=9.69 m/s^2

% How sensitive is the altitude reached after 4.5 minutes to simpler and
%   more complicated ways of modelling the gravitational field strength?
%
% Since the acceleration at 39,000 meters is 9.69 m/s^2 and at the surface
% of the earth is 9.81, we can conclude that a more complicated way of
% calculating the gravitiational field strength would not affect the value
% that much since the difference between the two is only 0.1m/s^2

% What other changes could we make to our model? Refer to, or at least
%   attempt to explain, the physics behind any changes that you propose.
%

% What is a change that we could make to our model that would result in
%   insignificant changes to the altitude reached after 4.5 minutes?
% <put your answer here in these comments>

% How can we decide what change is significant and what change is
%   insignificant?
% <put your answer here in these comments>

% [What changes did you try out to improve the model?  (Show us your changes
%   even if they didn't make the improvement you hoped for.)]
% <put your answer here in these comments>

part = 4;

% <place your work here>
%plotter(X,V,T,'Part 4')
%% Part 5
% Answer some questions here in these comments...
% At what altitude does Felix pull the ripcord to deploy his parachute?
% <put your answer here in these comments>

% Recalculate the ACd product with the parachute open, and modify your
%   code so that you use one ACd product before and one after this altitude.
%   According to this version of the model, what is the maximum magnitude
%   of acceleration that Felix experiences?
% <put your answer here in these comments>

%   How safe or unsafe would such an acceleration be for Felix?
% <put your answer here in these comments>

part = 5;

%Make a single acceleration-plot figure that includes, for each of the
%model and the acceleration calculated from measurements, the moment when
%the parachute opens and the following 10 or so seconds. If you have
%trouble solving this version of the model, just plot the acceleration
%calculated from measurements.
% <place your work here>

%% Part 6
% Answer some questions here in these comments...
% How long does it take for Felix’s parachute to open?
% 5.6 seconds, measured by looking at the video

part = 6;

%Redraw the acceleration figure from the previous Part but using the new
%   model. Also, using your plotting function from Part 1, plot the
%   measured/calculated data and the model for the entire jump from
%   stratosphere to ground.
% <place your work here>
endtime = 450
starttime =400
[t1,x1] = ode45(@fall, [starttime,endtime], [38969.4, 0]);
plotter(X,V,T,'Part6' )
%% nested functions
% nested functions below are required for the assignment.
% see Downey Section 10.1 for discussion of nested functions

function res = fall(t, X)
    %FALL <Summary of this function goes here>
    %   <Detailed explanation goes here>

    % do not modify this function unless required by you for some reason!

    p = X(1); % the first element is position
    v = X(2); % the second element is velocity

    dpdt = v; % velocity: the derivative of position w.r.t. time
    dvdt = acceleration(t, p, v); % acceleration: the derivative of velocity w.r.t. time

    res = [dpdt; dvdt]; % pack the results in a column vector
end

function res = acceleration(t, p, v)
    % <insert description of function here>
    % input...
    % t: time
    % p: position
    % v: velocity
    % output...
    % res: acceleration

    % do not modify this function unless required by you for some reason!

    a_grav = gravityEst(p);

    if part == 1 % variable part is from workspace of function main.
        res = -a_grav;

    elseif part == 6
            G=6.67*10^(-11);
                if t<263
                A=0.7
                end
                if t>263
                    A=25
                end
                C=1.15
                m=118;
                M=5.98*10^(24);
                r=6370000+p;
                R=6370000;
                g=-9.81;
                res1=g*(R)^2/r^2;
                res2=0.5*C*v^2*A*stdatmo(p);
                res2=res2/m
                res=res1+res2;

    else
        m = mass(t, v);
        b = drag(t, p, v, m);

        f_drag = b * v^2;
        a_drag = f_drag / m;
        res = -a_grav + a_drag;
    end
end

% Please paste in or type in code into the below functions as may be needed.

function a_grav = gravityEst(p)
    % estimate the acceleration due to gravity as a function of altitude, p
    A_GRAV_SEA = 9.807;  % acceleration of gravity at sea level in m/s^2

     if part <= 3 % fill in
        a_grav = A_GRAV_SEA;
     else
        a_grav= 3.98866*10^14/(p+6.77*10^6)^2
    end
end

function res = mass(t, v)
    % mass in kg of Felix and all his equipment
    res = 118
end

function res = drag(t, p, v, m)
% <insert description of function here>

    if part == 2
        res = 0.2;
    else
        %air resistance drag = 1/2*rho*A*Cd
        %%%%%%%%%%% input your code here %%%%%%%%%%%%%%%%
        m=118
        g=9.80665
        rho=stdatmo(p)
        a=0.5
        cd=0.2
        acd=(rho*a*cd)/2
        res=acd


        %ACd =2mg/rho*v*v

    end
end

function res = plotter(X,V,T,tit)
V=abs(V)*(1/3.6);
alt=x1(:,1);
vel=abs(x1(:,2));



dv=diff(vel);
dt=diff(t1);
acc=abs(dv./dt);



DV=diff(V)
DT=diff(T)
A=abs(DV./DT)
A=smooth(A)

hold on
grid on
plot(t1,alt,'bo');
plot(T,X,'go');
xlim([starttime,endtime]);
title(tit)
xlabel({'Time','(in Seconds)'})
ylabel({'Altitude','in meters'})
legend('Calculated','Excel')

figure
hold on
grid on
plot(t1,vel,'bo');
plot(T,V,'go');
xlim([starttime,endtime]);
title(tit)
xlabel({'Time','(in Seconds)'})
ylabel({'Velocity','in m/s'})
legend('Calculated','Excel')

figure
hold on
grid on
finalt1=t1(1:end-1);
finalT=T(1:end-1);
plot(finalt1,acc,'bo');
plot(finalT,A,'go');
xlim([starttime,endtime]);
title(tit)
xlabel({'Time','(in Seconds)'})
ylabel({'Acceleration','in m/s^2'})
legend('Calculated','Excel')
end


%% Additional nested functions
% Nest any other functions below.
%Do not put functions in other files when you submit.

% end of nested functions
end % closes function main.