A_66=zeros(6:6);
for i=1:6;
    j=1:6;
    A_66(i,j)=i*j;
end
fprintf('mean of A_66 = %1.2f\nstdev of A_66 = %1.2f\n',mean(A_66(:)),std(A_66(:)))

data=dlmread('US_energy_by_sector.csv',',',2,0);
h=figure();
plot(data(:,1),data(:,3),data(:,1),data(:,9));
%plotting residential energy consumption and transportational energy  
%consumption from 1949 to 2016 in trillions of BTUs
set(0, 'defaultAxesFontSize', 16)
set(0,'defaultTextFontSize',14)
set(0,'defaultLineLineWidth',3)
set(gcf, 'Position', [200, 200, 1000, 800])
xlabel('Time (years)')
ylabel('Total Energy Consumed (trillions of BTUs)')
legend('Residential','Transportation','Location','northwest')
title('US Energy Consumption from 1949 to 2016')
saveas(h,'figure_01.png')

rescum=cumsum(data(:,3))./1000000;%cumulative sum of residential energy
%consumption in quintillions of BTUs
transcum=cumsum(data(:,9))./1000000;%cumulative sum of transportational
%energy consumption in quintillions of BTUs
h=figure();
plot(data(:,1),rescum,data(:,1),transcum)
set(gcf, 'Position', [500, 200, 1000, 800])
xlabel('Time (years)')
ylabel('Total Energy Consumed (quintillions of BTUs)')
legend('Residential','Transportation','Location','northwest')
title('Cumulative US Energy Consumption from 1949 to 2016')
saveas(h,'figure_02.png')

g=figure();
hold on
timespan=30;

for h=[.1,1,5]
    freefall(h,timespan)
end
set(0, 'defaultAxesFontSize', 16)
set(0,'defaultTextFontSize',14)
set(0,'defaultLineLineWidth',3)
set(gcf, 'Position', [200, 200, 1000, 800])
xlabel('Time (seconds)')
ylabel('Velocity (meters/second)')
legend('v analytical(0.1)','v numerical(0.1)','v analytical(1)','v numerical(1)','v analytical(5)','v numerical(5)','Location','southeast')
title('Velocity Comparison by Varying Time Steps')

saveas(g,'figure01.png')

function [vx,vy,vz] = my_velocity(x,y,z,t)
    % Help documentation of "my_velocity"
    % This function computes the velocity in the x- and y-directions given
    % three vectors of position in x- and y-directions as a function of time
    % x = x-position
    % y = y-position
    % z = z-position
    % t = time
    % output
    % vx = velocity in x-direction
    % vy = velocity in y-direction
    % vz = velocity in z-direction

    vx=zeros(length(t),1);
    vy=zeros(length(t),1);
    vz=zeros(length(t),1);

    vx(1:end-1) = diff(x)./diff(t); % calculate vx as delta x/delta t
    vy(1:end-1) = diff(y)./diff(t); % calculate vy as delta y/delta t
    vz(1:end-1) = diff(z)./diff(t); % calculate vy as delta y/delta t

    vx(end) = vx(end-1);
    vy(end) = vy(end-1);
    vz(end) = vz(end-1);

end

function [ax,ay,az]=my_acceleration(x,y,z,t)
    % Help documentation of "my_acceleration"
    % This function computes the acceleration in the x- and y-directions given
    % three vectors of position in x- and y-directions as a function of time
    % x = x-position
    % y = y-position
    % z = z-position
    % t = time
    % output
    % ax = acceleration in x-direction
    % ay = acceleration in y-direction
    % az = acceleration in z-direction

    function v=diff_match_dims(x,t)
      v=zeros(length(t),1);
      v(1:end-1)=diff(x)./diff(t);
      v(end)=v(end-1);
    end

    [vx,vy,vz]=my_velocity(x,y,z,t);

    ax = diff_match_dims(vx',t);
    ay = diff_match_dims(vy',t);
    az = diff_match_dims(vz',t);
    %added apostrophe for derivatives of velocity

end