function [ output_args ] = tc_stab_parameter_guide_02( sweep_fqn, up_Kbit, down_Kbit )
%TC_STAB_PARAMETER_GUIDE Summary of this function goes here
%   try to read in the result from a ping sweep run
%	sweep_fqn: the log file of the ping sweep against the first hop after
%		the DSL link
%	up_Kbit: the uplink rate in Kilobits per second
%	down_Kbit: the downlink rate in Kilobits per second
%
% TODO:
%	find whether the carrier is ATM quantized
%	if yes:
%		what is the RTT step (try to deduce the combined up and down rates from this)
%		try to figure out the overhead for each packet
%	make sure all sizes are filled (use NaN for emty ones???)
%Thoughts:
%	maybe require to give the nominal up and down rates, to estimate the
%		RTT stepsize
%	ask about IPv4 or IPv6 (what about tunneling?)
%	the sweep should be taken directly connected to the modem to reduce
%		non-ATM routing delays

dbstop if error;
if ~(isoctave)
	timestamps.(mfilename).start = tic;
else
	tic();
end
disp(['Starting: ', mfilename]);

% control options
show_mean = 0;		% the means are noisier than the medians
show_median = 1;	% the mean is the way to go
show_min = 1;		% the min should be the best measure, but in the ATM test sweep is too variable
show_max = 0;		% only useful for debugging
show_sem = 0;		% give some estimate of the variance
show_ci = 1;		% show the confidence interval of the mean
ci_alpha = 0.05;	% alpha for confidence interval calculation
use_measure = 'median';
use_processed_results = 1;


if (nargin == 0)
	sweep_fqn = fullfile(pwd, 'ping_sweep_ATM.txt');	% was Bridged, LLC/SNAP RFC-1483/2684 connection (overhead 32 bytes - 14 = 18)
%  	sweep_fqn = fullfile(pwd, 'ping_sweep_CABLE_20120426_230227.txt');
%  	sweep_fqn = fullfile(pwd, 'ping_sweep_CABLE_20120801_001235.txt');
	up_Kbit = 512;
	down_Kbit = 3008;
end
if (nargin == 1)
	up_Kbit = 512;
	down_Kbit = 3008;
end
if (nargin == 2)
	down_Kbit = 3008;
end

%ATM
quantum.byte = 48;	% ATM packets are always 53 bytes, 48 thereof payload
quantum.bit = quantum.byte * 8;
ATM_cell.byte = 53;
ATM_cell.bit = ATM_cell.byte * 8;

% CONFIGURE THE NEXT TWO VALUES
line.down.Kbit = down_Kbit;
line.up.Kbit = up_Kbit;
% DONE
line.down.bit = line.down.Kbit * 1024;
line.up.bit = line.up.Kbit * 1024;


% known packet size offsets in bytes
offsets.IPv4 = 20;		% assume no IPv4 options are used, IP 6 would be 40bytes?
offsets.IPv6 = 40;		% not used yet...
offsets.ICMP = 8;		% ICMP header
offsets.ethernet = 14;	% ethernet header
% unknown offsets is what we need to figure out to feed tc-stab...


[sweep_dir, sweep_name] = fileparts(sweep_fqn);
cur_parsed_data_mat = [sweep_fqn(1:end-4), '.mat'];

if (use_processed_results && ~isempty(dir(cur_parsed_data_mat)))
	disp(['Loading processed ping data from ', cur_parsed_data_mat]);
	load(cur_parsed_data_mat, 'ping');
else
	% read in the result from a ping sweep
	disp(['Processing ping data from ', sweep_fqn]);
	ping = parse_ping_output(sweep_fqn);
	save(cur_parsed_data_mat, 'ping');
end


% analyze the data
min_ping_size = min(ping.data(:, ping.cols.size)) - offsets.ICMP;
disp(['Minimum size of ping payload used: ', num2str(min_ping_size), ' bytes.']);
known_overhead = offsets.ICMP + min_ping_size + offsets.IPv4;
ping.data(:, ping.cols.size) = ping.data(:, ping.cols.size) + known_overhead;	% we know we used ping so add the 8 bytes already (since these will not account for overhead)
size_list = unique(ping.data(:, ping.cols.size));

per_size.header = {'size', 'mean', 'median', 'min', 'max', 'std', 'n', 'sem', 'ci'};
per_size.cols = get_column_name_indices(per_size.header);
per_size.data = zeros([length(size_list), length(per_size.header)]);
per_size.data(:, per_size.cols.size) = size_list;
for i_size = 1 : length(size_list)
	cur_size = size_list(i_size);
	cur_size_idx = find(ping.data(:, ping.cols.size) == cur_size);
	per_size.data(i_size, per_size.cols.mean) = mean(ping.data(cur_size_idx, ping.cols.time));
	per_size.data(i_size, per_size.cols.median) = median(ping.data(cur_size_idx, ping.cols.time));
	per_size.data(i_size, per_size.cols.min) = min(ping.data(cur_size_idx, ping.cols.time));
	per_size.data(i_size, per_size.cols.max) = max(ping.data(cur_size_idx, ping.cols.time));
	per_size.data(i_size, per_size.cols.std) = std(ping.data(cur_size_idx, ping.cols.time), 0);
	per_size.data(i_size, per_size.cols.n) = length(cur_size_idx);
	per_size.data(i_size, per_size.cols.sem) = per_size.data(i_size, per_size.cols.std) / sqrt(length(cur_size_idx));	
	per_size.data(i_size, per_size.cols.ci) = calc_cihw(per_size.data(i_size, per_size.cols.std), per_size.data(i_size, per_size.cols.n), ci_alpha);	
end

figure('Name', sweep_name);
hold on;
legend_str = {};
if (show_mean)
	% means
	legend_str{end + 1} = 'mean';
	plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.mean), 'Color', [0 1 0 ]);
	if  (show_sem)
		legend_str{end + 1} = '+sem';
		legend_str{end + 1} = '-sem';
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.mean) - per_size.data(:, per_size.cols.sem), 'Color', [0 0.66 0]);
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.mean) + per_size.data(:, per_size.cols.sem), 'Color', [0 0.66 0]);
	end
	if  (show_ci)
		legend_str{end + 1} = '+ci';
		legend_str{end + 1} = '-ci';
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.mean) - per_size.data(:, per_size.cols.ci), 'Color', [0 0.37 0]);
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.mean) + per_size.data(:, per_size.cols.ci), 'Color', [0 0.37 0]);
	end

end
if(show_median)
	% median +- standard error of the mean, confidence interval would be
	% better
	legend_str{end + 1} = 'median';
	plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.median), 'Color', [1 0 0]);
	if (show_sem)
		legend_str{end + 1} = '+sem';
		legend_str{end + 1} = '-sem';
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.median) - per_size.data(:, per_size.cols.sem), 'Color', [0.66 0 0]);
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.median) + per_size.data(:, per_size.cols.sem), 'Color', [0.66 0 0]);
	end
	if (show_ci)
		legend_str{end + 1} = '+ci';
		legend_str{end + 1} = '-ci';
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.median) - per_size.data(:, per_size.cols.ci), 'Color', [0.37 0 0]);
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.median) + per_size.data(:, per_size.cols.ci), 'Color', [0.37 0 0]);
	end
end
if(show_min)
	% minimum, should be cleanest, but for the test data set looks quite sad...
	legend_str{end + 1} = 'min';
	plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.min), 'Color', [0 0 1]);
end
if(show_max)
	% minimum, should be cleanest, but for the test data set looks quite sad...
	legend_str{end + 1} = 'max';
	plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.max), 'Color', [0 0 0.66]);
end

title(['If this plot shows a (noisy) step function with a stepping of ', num2str(quantum.byte), ' bytes then the data carrier is quantised, make sure to use tc-stab']);
xlabel('Approximate packet size [bytes]');
ylabel('ICMP round trip times (ping RTT) [ms]');
legend(legend_str);
hold off;

% looking at the different plot we see many side peaks
%figure;
%plot(per_size.data(1:end - 1, strmatch('size', per_size.header , 'exact')), diff(per_size.data(:, strmatch('median', per_size.header , 'exact'))), 'Color', [1 0 0]);


% potentially clean up the data, by interpolating values with large sem
% from the neighbours?


% if this is ATM based the expectancy is there are RTT time quantums of 48
% bytes (the payload of an ATM package) find this
%x = fminsearch(@(x)sin(x^2), x0);


% estimate the RTT step size
% at ADSL down 3008kbit/sec up 512kbit/sec we expect, this oes not include
% processing time
expected_RTT_quantum_ms = (ATM_cell.bit / line.down.bit + ATM_cell.bit / line.up.bit ) * 1000;	% this estimate is rather a lower bound for fastpath , so search for best fits
disp(['lower bound estimate of one ATM cell RTT is ', num2str(expected_RTT_quantum_ms), ' ms.']);
% lets search from expected_RTT_quantum_ms to 1.5 * expected_RTT_quantum_ms
% in steps of expected_RTT_quantum_ms / 100

% to allow for interleaved ATM set ups increase the search space up to 32
% times best fastpath RTT estimate
RTT_quantum_list = (expected_RTT_quantum_ms : expected_RTT_quantum_ms / 100 : 32 * expected_RTT_quantum_ms);
quantum_list = (1 : 1 : quantum.byte);

differences = zeros([length(RTT_quantum_list) length(quantum_list)]);
cumulative_differences = differences;
all_stairs = zeros([length(RTT_quantum_list) length(quantum_list) length(per_size.data(:, per_size.cols.(use_measure)))]);
for i_RTT_quant = 1 : length(RTT_quantum_list)
	cur_RTT_quant = RTT_quantum_list(i_RTT_quant);
	for i_quant = 1 : quantum.byte
		[differences(i_RTT_quant, i_quant), cumulative_differences(i_RTT_quant, i_quant), all_stairs(i_RTT_quant, i_quant, :)] = get_difference_between_data_and_stair( per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.(use_measure)), ...
																	quantum_list(i_quant), quantum.byte, 0, cur_RTT_quant );
	end
end

% for the test DSL set the best x_offset is 21.
[min_cum_diff, min_cum_diff_idx] = min(cumulative_differences(:));
[min_cum_diff_row_idx, min_cum_diff_col_idx] = ind2sub(size(cumulative_differences),min_cum_diff_idx);
best_difference = differences(min_cum_diff_row_idx, min_cum_diff_col_idx);
disp(['remaining ATM cell length after ICMP header is ', num2str(quantum_list(min_cum_diff_col_idx)), ' bytes.']);
disp(['ICMP RTT of a single ATM packet is ', num2str(RTT_quantum_list(min_cum_diff_col_idx)), ' ms.']);

figure('Name', 'Comparing ping data with');
hold on
legend_str = {'ping_data', 'fitted_stair'};
plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.(use_measure)), 'Color', [1 0 0]);
plot(per_size.data(:, per_size.cols.size), squeeze(all_stairs(min_cum_diff_row_idx, min_cum_diff_col_idx, :)) + best_difference, 'Color', [0 1 0]);
title(['Estimated RTT per quantum: ', num2str(RTT_quantum_list(min_cum_diff_col_idx)), ' ms; ICMP header offset in quantum ', num2str(quantum_list(min_cum_diff_col_idx)), ' bytes']);
xlabel('Approximate packet size [bytes]');
ylabel('ICMP round trip times (ping RTT) [ms]');
if (isoctave)
	legend(legend_str);
else
	legend(legend_str, 'Interpreter', 'none');
end
hold off


% find cumulative_differences minimum, get the differences for that matrix
% and plot the corresponding stair with the real data
% next extrapolate to get the y-intercept


% as first approximation use the ATM cell offset and known offsets (ICMP
% IPv4 min_ping_size) to estimate the number of cells used for per paket
% overhead
% this assumes that no ATM related overhead is >= ATM cell size
% -1 to account for matlab 1 based indices
n_bytes_overhead_2nd_cell = quantum.byte - (quantum_list(min_cum_diff_col_idx) - 1);	% just assume we can not fit all overhead into one cell...
pre_IP_overhead = quantum.byte + (n_bytes_overhead_2nd_cell - known_overhead);	% ths is the one we are after in the end
disp(' ');
disp(['Estimated overhead preceeding the IP header: ', num2str(pre_IP_overhead), ' bytes']);



% use http://ace-host.stuart.id.au/russell/files/tc/tc-atm/ to present the
% most likely ATM setup for a given overhead and present a recommendation
% for the stab invocation
display_protocol_stack_information(pre_IP_overhead);



% now turn this into tc-stab recommendations:
disp(['Add the following to both the ingress and egress root qdiscs:']);
disp(' ');
disp(['A) Assuming the router connects over ethernet to the DSL-modem:']);
disp(['stab mtu 2048 tsize 128 overhead ', num2str(pre_IP_overhead - offsets.ethernet), ' linklayer atm']);
disp(' ');
disp(['B) Assuming the router connects via PPP and non-ethernet to the modem:']);
disp(['stab mtu 2048 tsize 128 overhead ', num2str(pre_IP_overhead), ' linklayer atm']);

disp(' ');
if ~(isoctave)
	timestamps.(mfilename).end = toc(timestamps.(mfilename).start);
	disp([mfilename, ' took: ', num2str(timestamps.(mfilename).end), ' seconds.']);
else
	toc
end

disp('Done...');

return
end


function [ ping_data ] = parse_ping_output( ping_log_fqn )
%PARSE_PING_OUTPUT read the putput of a ping run/sweep
% for further processng

cur_sweep_fd = fopen(ping_log_fqn, 'r');
ping_data.header = {'size', 'icmp_seq', 'ttl', 'time'};
ping_data.cols = get_column_name_indices(ping_data.header);

ping_data.data = zeros([1, length(ping_data.header)]);
cur_data_lines = 0;

% skip the first line
% PING netblock-75-79-143-1.dslextreme.com (75.79.143.1): (16 ... 1000)
% data bytes
header_line = fgetl(cur_sweep_fd);

while ~feof(cur_sweep_fd)
	cur_line = fgetl(cur_sweep_fd);
	[first_element, remainder] = strtok(cur_line);
	first_element_as_number = str2num(first_element);
	if isempty(first_element);
		% skip empty lines
		continue;
	end
	if strmatch('---', first_element)
		%we reached the end of sweeps
		break;
	end
	% now read in the data
	% 30 bytes from 75.79.143.1: icmp_seq=339 ttl=63 time=14.771 ms
	if ~isempty(first_element_as_number)
		cur_data_lines = cur_data_lines + 1;
		% size of the ICMP package
		ping_data.data(cur_data_lines, ping_data.cols.size) = first_element_as_number;
		% now process the remainder
		while ~isempty(remainder)
			[next_item, remainder] = strtok(remainder);
			equality_pos = strfind(next_item, '=');
			% data items are name+value pairs
			if ~isempty(equality_pos);
				cur_key = next_item(1: equality_pos - 1);
				cur_value = str2num(next_item(equality_pos + 1: end));
				switch cur_key
					% busybox ping and macosx ping return different key names
					case {'seq', 'icmp_seq'}
						ping_data.data(cur_data_lines, ping_data.cols.icmp_seq) = cur_value;
					case 'ttl'
						ping_data.data(cur_data_lines, ping_data.cols.ttl) = cur_value;
					case 'time'
						ping_data.data(cur_data_lines, ping_data.cols.time) = cur_value;
				end
			end
		end
		
	end
	
end

% clean up
fclose(cur_sweep_fd);
return
end



function [ difference , cumulative_difference, stair_y ] = get_difference_between_data_and_stair( data_x, data_y, x_size, stair_x_step_size, y_offset, stair_y_step_size )
% x_size is the flat part of the first stair, that is quantum minus the
% offset
debug = 0;
difference = [];
x_start_val = min(data_x);
x_end_val = max(data_x);

% construct stair
stair_x = data_x;
proto_stair_y = zeros([data_x(end) 1]);
% make sure the x_size values do not exceed the step size...
if (x_size > stair_x_step_size)
	if mod(x_size, stair_x_step_size) == 0
		x_size = stair_x_step_size;
	else
		x_size = mod(x_size, stair_x_step_size);
	end
end

stair_y_step_idx = (x_start_val + x_size : stair_x_step_size : x_end_val);
proto_stair_y(stair_y_step_idx) = stair_y_step_size;
stair_y = cumsum(proto_stair_y);
stair_y = stair_y(x_start_val:x_end_val) + y_offset;
if (debug)
	figure
	hold on;
	title(['x offset used: ', num2str(x_size), ' with quantum ', num2str(stair_x_step_size)]);
	plot(data_x, data_y, 'Color', [0 1 0]);
	plot(stair_x, stair_y, 'Color', [1 0 0]);
	hold off;
end
difference = sum(abs(data_y - stair_y)) / length(data_y);
cumulative_difference = sum(abs(data_y - (stair_y + difference)));

return
end

% function [ stair ] = build_stair(x_vector, x_size, stair_x_step_size, y_offset, stair_y_step_size )
% stair = [];
% 
% return
% end


function [columnnames_struct, n_fields] = get_column_name_indices(name_list)
% return a structure with each field for each member if the name_list cell
% array, giving the position in the name_list, then the columnnames_struct
% can serve as to address the columns, so the functions assitgning values
% to the columns do not have to care too much about the positions, and it
% becomes easy to add fields.
n_fields = length(name_list);
for i_col = 1 : length(name_list)
	cur_name = name_list{i_col};
	columnnames_struct.(cur_name) = i_col;
end
return
end




function [ci_halfwidth_vector] = calc_cihw(std_vector, n, alpha)
%calc_ci : calculate the half width of the confidence interval (for 1 - alpha)
%	the t_value lookup depends on alpha and the samplesize n; the relevant
%	calculation of the degree of freedom is performed inside calc_t_val.
%	ci_halfwidth = t_val(alpha, n-1) * std / sqrt(n)
%	Each groups CI ranges from mean - ci_halfwidth to mean - ci_halfwidth, so
%	the calling function has to perform this calculation...
%
% INPUTS:
%	std_vector: vector containing the standard deviations of all requested
%		groups
%	n: number of samples in each group, if the groups have different
%		samplesizes, specify each groups samplesize in a vector
%	alpha: the desired maximal uncertainty/error in the range of [0, 1]
% OUTPUT:
%	ci_halfwidth_vector: vector containing the confidence intervals half width
%		for each group

% calc_t_val return one sided t-values, for the desired two sidedness one has
% to half the alpha for the table lookup
cur_alpha = alpha / 2;

% if n is scalar use same n for all elements of std_vec
if isscalar(n)
	t_ci = calc_t_val(cur_alpha, n);
	ci_halfwidth_vector = std_vector * t_ci / sqrt(n);
	% if n is a vector, prepare a matching vector of t_ci values
elseif isvector(n)
	t_ci_vector = n;
	% this is probably ugly, but calc_t_val only accepts scalars.
	for i_pos = 1 : length(n)
		t_ci_vector(i_pos) = calc_t_val(cur_alpha, n(i_pos));
	end
	ci_halfwidth_vector = std_vector .* t_ci_vector ./ sqrt(n);
end

return
end

%-----------------------------------------------------------------------------
function [t_val] = calc_t_val(alpha, n)
% the t value for the given alpha and n
% so call with the n of the sample, not with degres of freedom
% see http://mathworld.wolfram.com/Studentst-Distribution.html for formulas
% return values follow Bortz, Statistik fuer Sozialwissenschaftler, Springer 
% 1999, table D page 775. That is it returns one sided t-values.
% primary author S. Moeller

% TODO:
%	sidedness of t-value???

% basic error checking
if nargin < 2
	error('alpha and n have to be specified...');
end

% probabilty of error
tmp_alpha = alpha ;%/ 2;
if (tmp_alpha < 0) || (tmp_alpha > 1)
	msgbox('alpha has to be taken from [0, 1]...');
	t_val = NaN;
	return
end
if tmp_alpha == 0
	t_val = -Inf;
	return
elseif tmp_alpha ==1
	t_val = Inf;
	return
end
% degree of freedom
df = n - 1;
if df < 1
	%msgbox('The n has to be >= 2 (=> df >= 1)...');
% 	disp('The n has to be >= 2 (=> df >= 1)...');
	t_val = NaN;
	return
end


% only calculate each (alpha, df) combination once, store the results
persistent t_val_array;
% create the t_val_array
if ~iscell(t_val_array)
	t_val_array = {[NaN;NaN]};
end
% search for the (alpha, df) tupel, avoid calculation if already stored
if iscell(t_val_array)
	% cell array of 2d arrays containing alpha / t_val pairs
	if df <= length(t_val_array)
		% test whether the required alpha, t_val tupel exists
		if ~isempty(t_val_array{df})
			% search for alpha
			tmp_array = t_val_array{df};
			alpha_index = find(tmp_array(1,:) == tmp_alpha);
			if any(alpha_index)
				t_val = tmp_array(2, alpha_index);
				return
			end
		end
	else
		% grow t_val_array to length of n
		missing_cols = df - length(t_val_array);
		for i_missing_cols = 1: missing_cols
			t_val_array{end + 1} = [NaN;NaN];
		end
	end
end

% check the sign
cdf_sign = 1;
if (1 - tmp_alpha) == 0.5
	t_val = t_cdf;
elseif (1 - tmp_alpha) < 0.5 % the t-cdf is point symmetric around (0, 0.5)
	cdf_sign = -1;
	tmp_alpha = 1 - tmp_alpha; % this will be undone later
end

% init some variables
n_iterations = 0;
delta_t = 1;
last_alpha = 1;
higher_t = 50;
lower_t = 0;
% find a t-value pair around the desired alpha value 
while norm_students_cdf(higher_t, df) < (1 - tmp_alpha);
	lower_t = higher_t;
	higher_t = higher_t * 2;
end

% search the t value for the given alpha...
while (n_iterations < 1000) && (abs(delta_t) >= 0.0001)
	n_iterations = n_iterations + 1;
	% get the test_t (TODO linear interpolation) 
 	% higher_alpha = norm_students_cdf(higher_t, df);
 	% lower_alpha = norm_students_cdf(lower_t, df);
	test_t = lower_t + ((higher_t - lower_t) / 2);
	cur_alpha = norm_students_cdf(test_t, df);
	% just in case we hit the right t spot on...
	if cur_alpha == (1 - tmp_alpha)
		t_crit = test_t;
		break;
	% probably we have to search for the right t
	elseif cur_alpha < (1 - tmp_alpha)
		% test_t is the new lower_t
		lower_t = test_t;
		%higher_t = higher_t;	% this stays as is...
	elseif cur_alpha > (1 - tmp_alpha)
		% 
		%lower_t = lower_t;	% this stays as is...
		higher_t = test_t;
	end
	delta_t = higher_t - lower_t;
	last_alpha = cur_alpha;
end
t_crit = test_t;

% set the return value, correct for negative t values
t_val = t_crit * cdf_sign;
if cdf_sign < 0
	tmp_alpha = 1 - tmp_alpha;
end

% store the alpha, n, t_val tupel in t_val_array
pos = size(t_val_array{df}, 2);
t_val_array{df}(1, (pos + 1)) = tmp_alpha;
t_val_array{df}(2, (pos + 1)) = t_val;

return
end

%-----------------------------------------------------------------------------
function [scaled_cdf] = norm_students_cdf(t, df)
% calculate the cdf of students distribution for a given degree of freedom df,
% and all given values of t, then normalize the result
% the extreme values depend on the values of df!!!

% get min and max by calculating values for extrem t-values (e.g. -10000000,
% 10000000)
extreme_cdf_vals = students_cdf([-10000000, 10000000], df);

tmp_cdf = students_cdf(t, df);

scaled_cdf =	(tmp_cdf - extreme_cdf_vals(1)) /...
				(extreme_cdf_vals(2) - extreme_cdf_vals(1));
return
end

%-----------------------------------------------------------------------------
function [cdf_value_array] = students_cdf(t_value_array, df)
%students_cdf: calc the cumulative density function for a t-distribution
% Calculate the CDF value for each value t of the input array
% see http://mathworld.wolfram.com/Studentst-Distribution.html for formulas
% INPUTS:	t_value_array:	array containing the t values for which to
%							calculate the cdf
%			df:	degree of freedom; equals n - 1 for the t-distribution

cdf_value_array = 0.5 +...
		((betainc(1, 0.5 * df, 0.5) / beta(0.5 * df, 0.5)) - ...
		(betainc((df ./ (df + t_value_array.^2)), 0.5 * df, 0.5) /...
		beta(0.5 * df, 0.5))) .*...
		sign(t_value_array);
	
return
end

%-----------------------------------------------------------------------------
function [t_prob_dist] = students_pf(df, t_arr)
%  calculate the probability function for students t-distribution

t_prob_dist =	(df ./ (df + t_arr.^2)).^((1 + df) / 2) /...
				(sqrt(df) * beta(0.5 * df, 0.5));

% % calculate and scale the cdf by hand...
% cdf = cumsum(t_prob_dist);
% discrete_t_cdf = (cdf - min(cdf)) / (max(cdf) - min(cdf));
% % numericaly get the t-value for the given alpha
% tmp_index = find(discrete_t_cdf > (1 - tmp_alpha));
% t_crit = t(tmp_index(1));
				
return
end

function in = isoctave ()
 persistent inout;

 if isempty(inout),
   inout = exist('OCTAVE_VERSION','builtin') ~= 0;
 end;
 in = inout;

return;
end


function [] = display_protocol_stack_information(pre_IP_overhead)
% use http://ace-host.stuart.id.au/russell/files/tc/tc-atm/ to present the
% most likely ATM prtocol stack setup for a given overhead so the user can
% compare with his prior knowledge

disp(' ');
disp('According to http://ace-host.stuart.id.au/russell/files/tc/tc-atm/');
disp(['', num2str(pre_IP_overhead), ' bytes overhead indicate']);

switch pre_IP_overhead
	case 8
		disp('Connection: IPoA, VC/Mux RFC-2684');
		disp('Protocol (bytes): ATM AAL5 SAR (8) : 8');
		
	case 16
		disp('Connection: IPoA, LLC/SNAP RFC-2684');
		disp('Protocol (bytes): ATM LLC (3), ATM SNAP (5), ATM AAL5 SAR (8) : 16');
		
	case 24
		disp('Connection: Bridged, VC/Mux RFC-1483/2684');
		disp('Protocol (bytes): Ethernet Header (14), ATM pad (2), ATM AAL5 SAR (8) : 24');

	case 28
		disp('Connection: Bridged, VC/Mux+FCS RFC-1483/2684');
		disp('Protocol (bytes): Ethernet Header (14), Ethernet PAD [8] (0), Ethernet Checksum (4), ATM pad (2), ATM AAL5 SAR (8) : 28');
	
	case 32
		disp('Connection: Bridged, LLC/SNAP RFC-1483/2684');
		disp('Protocol (bytes): Ethernet Header (14), ATM LLC (3), ATM SNAP (5), ATM pad (2), ATM AAL5 SAR (8) : 32');
		disp('OR');
		disp('Connection: PPPoE, VC/Mux RFC-2684');
		disp('Protocol (bytes): PPP (2), PPPoE (6), Ethernet Header (14), ATM pad (2), ATM AAL5 SAR (8) : 32');
	
	case 36
		disp('Connection: Bridged, LLC/SNAP+FCS RFC-1483/2684');
		disp('Protocol (bytes): Ethernet Header (14), Ethernet PAD [8] (0), Ethernet Checksum (4), ATM LLC (3), ATM SNAP (5), ATM pad (2), ATM AAL5 SAR (8) : 36');
		disp('OR');
		disp('Connection: PPPoE, VC/Mux+FCS RFC-2684');
		disp('Protocol (bytes): PPP (2), PPPoE (6), Ethernet Header (14), Ethernet PAD [8] (0), Ethernet Checksum (4), ATM pad (2), ATM AAL5 SAR (8) : 36');
	
	case 10
		disp('Connection: PPPoA, VC/Mux RFC-2364');
		disp('Protocol (bytes): PPP (2), ATM AAL5 SAR (8) : 10');
	
	case 14  
		disp('Connection: PPPoA, LLC RFC-2364');
		disp('Protocol (bytes): PPP (2), ATM LLC (3), ATM LLC-NLPID (1), ATM AAL5 SAR (8) : 14');
	
	case 40
		disp('Connection: PPPoE, LLC/SNAP RFC-2684');
		disp('Protocol (bytes): PPP (2), PPPoE (6), Ethernet Header (14), ATM LLC (3), ATM SNAP (5), ATM pad (2), ATM AAL5 SAR (8) : 40');
	
	case 44
		disp('Connection: PPPoE, LLC/SNAP+FCS RFC-2684');
		disp('Protocol (bytes): PPP (2), PPPoE (6), Ethernet Header (14), Ethernet PAD [8] (0), Ethernet Checksum (4), ATM LLC (3), ATM SNAP (5), ATM pad (2), ATM AAL5 SAR (8) : 44');  
	
	otherwise
		disp('a protocol stack this program does know (yet)...');
end

disp(' ');
return;
end

