%Below is an example MATLAB script that reads an ASCII text file of numeric data, 
% plots its histogram (normalized as a probability density function), 
% fits a Gaussian function to the histogram data, and then prints the fit parameters. 
% In this example, the Gaussian model is defined as f(x)=aexp⁡(−(x−b)/c)2) 
% where b is the mean and c is related to the standard deviation by σ=c/SQRT(2) ​.

% Clear workspace and close figures
clear; close all; clc;

% Prompt the user for the input file name (e.g., 'data.txt')
filename = input('Enter the input ASCII text file name (with extension): ', 's');

% Load numeric data from the file
data = load(filename);

% Determine the optimal bin width using the Freedman-Diaconis rule:
%   bin width = 2 * IQR(data) / n^(1/3)
n = length(data);
iqr_val = iqr(data);
binWidth = 2 * iqr_val / n^(1/3);

% Calculate the number of bins needed
numBins = ceil((max(data) - min(data)) / binWidth);

% Plot the histogram with optimized bins, normalized to a probability density function
figure;
h = histogram(data, numBins, 'Normalization', 'pdf');
hold on;
%title('Histogram with Gaussian Fit');
xlabel('Nanoparticle Displacements (\mum)');
ylabel('Count');
xlim([-0.8 0.8]); %CHANGES THIS

% Calculate bin centers from the bin edges
edges = h.BinEdges;
centers = (edges(1:end-1) + edges(2:end)) / 2;

% Define the Gaussian model as a string for use with fit:
% f(x) = a * exp(-((x-b)/c)^2)
gaussEqn = 'a*exp(-((x-b)/c)^2)';

% Set initial guesses for the parameters:
% a: approximate peak height from the histogram
% b: mean of the data
% c: related to the standard deviation (sigma = c/sqrt(2))
a0 = max(h.Values);
b0 = mean(data);
c0 = std(data)*sqrt(2);

% Fit the Gaussian model to the histogram data (using the bin centers and counts)
f = fit(centers', h.Values', gaussEqn, 'Start', [a0, b0, c0]);

% Generate fine x-data for plotting the fitted curve
xFit = linspace(min(data), max(data), 10000);
yFit = f.a * exp(-((xFit - f.b) / f.c).^2);
plot(xFit, yFit, 'r-', 'LineWidth', 2);
hold off;

% Obtain the 95% confidence intervals for the fit parameters.
% confint returns a 2-by-3 matrix: first row is the lower bounds, second row is the upper bounds.
CI = confint(f, 0.95);
error_a = (CI(2,1) - CI(1,1)) / 2;
error_b = (CI(2,2) - CI(1,2)) / 2;
error_c = (CI(2,3) - CI(1,3)) / 2;

% Output the fit parameters and their errors
fprintf('Gaussian Fit Parameters (with 95%% confidence intervals):\n');
fprintf('Amplitude, a = %.4f ± %.4f\n', f.a, error_a);
fprintf('Mean, b = %.4f ± %.4f\n', f.b, error_b);
fprintf('Parameter, c = %.4f ± %.4f (sigma = c/sqrt(2))\n', f.c, error_c);


% Output the fit parameters
fprintf('Gaussian Fit Parameters:\n');
fprintf('Amplitude, a = %.4f\n', f.a);
fprintf('Mean, b = %.4f\n', f.b);
fprintf('Parameter, c = %.4f (sigma = c/sqrt(2))\n', f.c);
