%% load all the .csv files in path folder and put them into a cell array

files = dir('*.csv');   % dir struct of all pertinent .csv files
n = length(files);      % number of .csv files

% matlab may not sort the files correctly by their names (check 'files' in workspace)
% use this to sort by date if needed
[~,index] = sortrows([files.datenum].');
files = files(index);
clear index
%% continue
data = cell(n,1);  % preallocate cell array for all the data

for i = 1:n
    signal = readtable(files(i).name);  % read each file
    signal = table2array(signal);  % convert to array
    data{i} = signal;  % save
end

% individual file can be accessed using
%signal = readtable('filename.csv');

%% gain computation

gains = zeros(n,1);  % preallocate array

% find extremes for each signal and compute the gain of the PM tube
for i = 1:n
    signal = data{i};
    CH1 = signal(:,2);
    CH2 = signal(:,3);

    gain = (abs(min(CH2))/max(CH1));
    gains(i) = gain;
end

%% use this to account for signals where attenuator was used (bounds set manually)
% gains(1:?) = gains(1:?)/10;  % careful not to repeatedly run this part

%% gain results processing
ns = 2;  % number of samples per voltage meassured

results = zeros(2,n/ns);  % preallocate space for the result: 2 rows for mean and SD and n/ns columns for each voltage

temp = 1;
% compute mean and standard deviation for number of samples of each voltage
for i = 1:ns:n
    samples = gains(i:i+ns-1);  % extract corresponding gains
    
    results(1,temp) = mean(samples);
    results(2,temp) = std(samples);
    temp = temp + 1;
end

% save matrix as excel file to further process the results
writematrix(results, 'gains.xls')  % .csv puts the resutls array to one cell

%% compute time interval between both extremes

times = zeros(n,1);  % preallocate array

% find extremes for each signal and compute the gain of the PM tube
for i = 1:n
    signal = data{i};
    time = signal(:,1);
    CH1 = signal(:,2);
    CH2 = signal(:,3);

    [~, I] = max(CH1);  % find index of the extreme
    t1 = time(I);       % find corresponding time
    
    [~, I] = min(CH2);  % find index of the extreme
    t2 = time(I);       % find corresponding time
    
    t = t2 - t1;  % time difference between the extreems
    
    times(i) = t;
end

%% transition time computation

% set accoring to your setup
PMT_cable_time = 1.105*5 + 0.058*3 + 1.5*5;
PD_cable_time = 0.08*3.3 + 1.5*5;

% convert to ns
times = times*10^9;

% PMT transition time = PMT_cable_time + extremes_time_diff - PD_cable_time
tr_time = zeros(n,1);  % prealloc

for i = 1:n
    tr_time(i) = times(i) + PMT_cable_time - PD_cable_time;
end

%% transition time result processing
ns = 2;  % number of samples per voltage meassured

t_results = zeros(2,n/ns);  % preallocate space for the result: 2 rows for mean and SD and n/ns columns for each transition time

temp = 1;
% compute mean and standard deviation for number of samples of each tr_time
for i = 1:ns:n
    samples = tr_time(i:i+ns-1);  % extract corresponding transition times
    
    t_results(1,temp) = mean(samples);
    t_results(2,temp) = std(samples);
    temp = temp + 1;
end

% save matrix as excel file to further process the results
writematrix(t_results, 'times.xls')  % .csv puts the resutls array to one cell

%% compute area under the curve
% access cell data
signal = data{1,1};

% define channels
x = signal(:,1);  % time
y1 = signal(:,2);  % CH1
y2 = signal(:,3);  % CH2
y3 = signal(:,4);  % CH3
y4 = signal(:,5);  % CH4

% integrate (time set manually)
format long
lims = (x >= -10*10^(-9)) & (x <= 90*10^(-9));  % 100 ns

integral = [trapz(x(lims),y1(lims)); trapz(x(lims),y2(lims)); trapz(x(lims),y3(lims)); trapz(x(lims),y4(lims))];
format
%% convert signal to graph

% choose data
signal = data{1,1};

time = signal(:,1);

% choose channels
CH1 = signal(:,2);
CH2 = signal(:,3);
% CH3 = signal(:,3);
% CH4 = signal(:,4);

% write it to excel
graph = [time CH1 CH2];
writematrix(graph, 'graph.xls')

% or plot it here
plot(time,CH1)

grid on

%daspect([1 10^(8.2) 1])  % increase power of 10 to widen the graph

% graph visualisation
axis padded

% change units
xtickformat('%d ns')
ytickformat('%d mV')
ax = gca;
ax.XAxis.Exponent = -9;
ax.YAxis.Exponent = -3;

figure
set(gca, 'FontSize', 15)
% settings for the thesis latex version
% set(groot, 'defaultAxesTickLabelInterpreter', 'latex')
% set(groot, 'defaultTextInterpreter', 'latex')
% set(groot, 'defaultLegendInterpreter', 'latex')

% set specific numbers on xy axes
%xticks([-300 -200 -100 0 100 200 300 400 500 600 700]*10^(-9))
%yticks([-12 -10 -8 -6 -4 -2 0])


