%%% Filename:   flight_data_analysis_060110.m
%%% Author:     Brock J. LaMeres
%%% Date:       6/25/10
%%%
%%% Description:	This script analyzes data taken on a mock
%%%                 high altitude balloon flight which recorded
%%%                 radiation levels vs. altitude
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% First, let's clean up the workspace to avoid variable issues
clc             % clear command window
clear           % clear variables in the workspace
warning off     % suppress warnings, the polyfit/polyvar commands 
                %    can generate warnings



%%% Read in the raw data file into an m x n matrix where 
%%% column 1 = altitude (ft) and column 2 = radiation (MeV)
balloon_data = load('alt_v_rad.dat');
 
%%% Now let's parse the data into column vectors with descriptive names
altitude  = balloon_data(:,1);
radiation = balloon_data(:,2);

%%% Let's plot the altitude vs. radiation
subplot (2,3,1)
plot(altitude, radiation, '.')
title('Radiation Sensor Flight Data - 060110')
xlabel('Altitude (ft)')
ylabel('Radiation Level (MeV)')
hold on

%%% Let's now perform some basic statistical analysis of the data and
%%% print it to the screen
mean_radiation   = mean(radiation);                    % Find mean
median_radiation = median(radiation);                  % Find median
[max_radiation, max_radiation_index] = max(radiation); % Find max and indx

fprintf('\n\n')
fprintf('----------------------------------------------------\n')
fprintf('---            Flight Data Analysis              ---\n')
fprintf('----------------------------------------------------\n')
fprintf(' Mean radiation level:\t\t %2.2f MeV\n', mean_radiation);
fprintf(' Median radiation level:\t %2.2f MeV\n', median_radiation);
fprintf(' Maximum radiation level:\t %2.2f MeV (@ %2.0f ft)\n', max_radiation, altitude(max_radiation_index) );
fprintf('\n')
fprintf('----------------------------------------------------\n')

%%% Let's perform a regression analysis on the data using the function
%%% "polyfit".  This will allow us to construct a function which matches
%%% the data.  In this analysis, we assume that altitude is independant.

%%% First, let's calculate the polynomial coefficents using "polyfit"
poly_coef_1st_order = polyfit(altitude, radiation, 1);
poly_coef_2nd_order = polyfit(altitude, radiation, 2);
poly_coef_3rd_order = polyfit(altitude, radiation, 3);
poly_coef_4th_order = polyfit(altitude, radiation, 4);
poly_coef_5th_order = polyfit(altitude, radiation, 5);

%%% Now let's generate the polynomial values against a tightly spaced
%%% independant variable (alt_ideal) using the function "polyval"

alt_ideal = [0:10:100000];

poly_val_1st_order = polyval(poly_coef_1st_order, alt_ideal);
poly_val_2nd_order = polyval(poly_coef_2nd_order, alt_ideal);
poly_val_3rd_order = polyval(poly_coef_3rd_order, alt_ideal);
poly_val_4th_order = polyval(poly_coef_4th_order, alt_ideal);
poly_val_5th_order = polyval(poly_coef_5th_order, alt_ideal);

%%% Now let's plot it over the original data to see how it matches 
plot(alt_ideal, poly_val_1st_order, 'k')    
plot(alt_ideal, poly_val_2nd_order, 'c')    
plot(alt_ideal, poly_val_3rd_order, 'g')    
plot(alt_ideal, poly_val_4th_order, 'm')    
plot(alt_ideal, poly_val_5th_order, 'r')    
% Let's put a legend on the plot, notice that I can specify the location
legend ('original data', '1st order', '2nd order', '3rd order', '4th order', '5th order', 'Location', 'NorthWest')


%%% Let's plot some histograms to analyze the number of radiation strikes

%%%  First let's use the built in histogram function
%%%  With no options provided, this function will sort the data into 
%%%  10 evenly spaced bins between the minimum and maximum.

subplot (2,3,2)
hist(radiation)
title('Histogram of Flight Data -  Default Plotting')
xlabel('Radiation Level (MeV)')
ylabel('# of Occurances')


%%% Now, let's create a new histogram that has more customized bins.
%%% We can use hist(y,n) which sort the data into 'n' bins between the
%%% minimum and maximum values in 'y'

subplot (2,3,3)
hist(radiation, 6);
title('Histogram of Flight Data - Custom # of bins')
xlabel('Radiation Level (MeV)')
ylabel('# of Occurances')

%%% To get a different look, we can plot the data using 'bar'.
%%% The 'hist' command can still be used but instead of plotting, we will
%%% put the results into new vectors [z,x] = hist(y,n)
%%%    where: n provides the number of bins to sort the y into
%%%           y is the original data
%%%           z is a vector of size n that contains the counts per bin
%%%           x is are the bin locations

subplot(2,3,4)
[radiation_count, radiation_bins] = hist(radiation,6);
bar(radiation_bins,radiation_count)
title('Histogram of Flight Data - Using Bar Plot')
xlabel('Radiation Level (MeV)')
ylabel('# of Occurances')

%%% We can use different sizes of bins in the histogram by providing a
%%% vector that contains the desired bin centers.  The bin widths are
%%% the differences between centers.  

subplot(2,3,5)
hist(radiation, [1, 1.5, 2.25, 2.5, 2.75, 3 ]);  % this creates a large bin at 
title('Histogram of Flight Data - Different bin sizes')
xlabel('Radiation Level (MeV)')
ylabel('# of Occurances')


%%% "Scaled Histogram" it is useful from a probability standpoint to
%%% scale the histogram so that the total area of each histogram = 1.
%%% Using this scaling, instead of the y-axis containing the absolute
%%% frequency, it contains the relative frequency, or the percentage
%%% of the time that the bin occurs.  This can be extended to probability.

radiation_count_sum = sum(radiation_count);
radiation_count_scaled = radiation_count/radiation_count_sum;
subplot(2,3,6)
bar(radiation_bins,radiation_count_scaled)
title('Histogram of Flight Data - Scaled for Relative Frequency')
xlabel('Radiation Level (MeV)')
ylabel('% of Occurances')

%%% cumsub(y) is a command that calculates the area under the scaled
%%% frequency data and hence, probability. The returned value is the
%%% successive sum of everything before it.

rel_freq_prob = cumsum(radiation_count_scaled)*100;

fprintf('\n')
fprintf('----------------------------------------------------\n')
fprintf('---            Flight Data Probability           ---\n')
fprintf('----------------------------------------------------\n')
fprintf(' Probability of a strike below 0.5 MeV is:\t %3.2f %% \n', rel_freq_prob(1));
fprintf(' Probability of a strike below 1.0 MeV is:\t %3.2f %% \n', rel_freq_prob(2));
fprintf(' Probability of a strike below 1.5 MeV is:\t %3.2f %% \n', rel_freq_prob(3));
fprintf(' Probability of a strike below 2.0 MeV is:\t %3.2f %% \n', rel_freq_prob(4));
fprintf(' Probability of a strike below 2.5 MeV is:\t %3.2f %% \n', rel_freq_prob(5));
fprintf(' Probability of a strike below 3.0 MeV is:\t %3.2f %% \n', rel_freq_prob(6));
fprintf('\n')
fprintf('----------------------------------------------------\n')

hold off



