%%
ref = load('qRef_13_21_32.mat');
filename1 = 'ref';
files = dir('*.mat');
N = numel(files)

[tAxis, ref_s, ref_f, ref_X] = sortData(ref);
ref_s = ref_s - trimmean(ref_s, 0.1);
[w,ws_ref_s] = my_window_function(tAxis'*1e-12,ref_s);
[p,X_ref_s] = my_fft(tAxis,ref_s);
[~,X_w_ref_s] = my_fft(tAxis,ws_ref_s');
figure();
plot(tAxis, ref_s, tAxis, ws_ref_s, tAxis, w.*max(ref_s))

%%
% test = (0:0.5:11);
% knifeposition = flip(test);
%Pmax = 1;
%P0 = -206.26;
%x0 = 8.2;
resultTable = table('Size', [0 6], 'VariableTypes', {'string', 'double', 'double', 'double','double','double'}, 'VariableNames', {'Filename', 'YValue', 'MinValueInRange', 'XIndexMinValue','MaxValueInRange', 'XIndexMaxValue'});
% Get list of data files in directory
files = dir('*.mat');
colors = lines(numel(files)); % create a matrix of colors
T = table();
figure;

for i = 1:numel(files)
    data = load(files(i).name);
    
    timeAxis = data.tAxis;
    Power = data.sAxis;
    Power(isnan(Power))= 0;
    Power = Power - trimmean(Power, 0.1);
    [w,ws_zs] = my_window_function(timeAxis'*1e-12,Power);% 5.3022e-11); %5.1021e-11
    [p,X_zs] = my_fft(timeAxis,Power);
    [p,X_w_zs] = my_fft(timeAxis,ws_zs');
    
    X_w_zs = abs(X_w_zs);
    
    %Referenced windowed intensity
    plot(p,(abs(X_w_zs.^2)./abs(X_w_ref_s.^2)))
    ylim([0,1])
    xlim([0.4,1.4])

    % Interpolate to find the y-axis value at x = 0.5
    yvalue = interp1(p, (abs(X_w_zs.^2)./abs(X_w_ref_s.^2)), 0.5)

    ydata = (abs(X_w_zs.^2)./abs(X_w_ref_s.^2));

    MinRange = (p >= 0.9 & p <= 1);
    MaxRange = (p >= 1 & p <= 1.1);

    [minValueInRange, minIdx] = min(ydata(MinRange));
    [maxValueInRange, maxIdx] = max(ydata(MaxRange));

    xIndexMinValue = p(MinRange);
    xIndexMinValue = xIndexMinValue(minIdx);

    xIndexMaxValue = p(MaxRange);
    xIndexMaxValue = xIndexMaxValue(maxIdx);
    
    newRow = {files(i).name, yvalue, minValueInRange, xIndexMinValue,maxValueInRange,xIndexMaxValue};
    resultTable = [resultTable; newRow];

    % Add legend
    if i == 1
        legend(files(i).name)
        hold on
    else
        legend(files(i).name)
    end
end

% Display the values

% Add legend outside the loop
legend(files.name, 'Location', 'best')
%% Functions

function [timeAxis,sBinned,frequency,FFT] = sortData(dat)
    sBinned = dat.sAxis;
    sBinned(isnan(sBinned))= 0;
    timeAxis = dat.tAxis;
    frequency = dat.f;
    FFT = dat.X;
end

function [f,X] = my_fft(t,s)
    
    nfft = pow2(nextpow2(length(t))+6); %number of points in f axis / zero padding
    T = t(2)-t(1); 
    fs = 1/(T*nfft); 
    X = fft(s',nfft); 
    X = X(1:nfft/2);
    f =(0:nfft/2-1)*fs;           % Frequency
    

end
function [w,ws] = my_window_function(t,s,peak,A,B,offset)
% Inputs:
%   t - time axis (s)
%   s - signal
%   peak (optional) - location of the peak of the THz pulse
%   A (optional) - scaling parameter for changing the steepness
%   B (optional) - scaling the parater for changing the separation
%   offset (optional) - offets the window to alter the peak

if ~exist('peak','var')
     % third parameter does not exist, so default it to something
    [m,ind]=max(s);
    for i=1:length(ind)
        
        peak(i)=t(ind(i),i);
        peak
    end
end
if ~exist('B','var')
     % third parameter does not exist, so default it to something
    B=8e-12;
end
if ~exist('offset','var')
     % third parameter does not exist, so default it to something
    offset=0;%1e-12
end
if ~exist('A','var')
     % third parameter does not exist, so default it to something
    A=0.2e-12;
end
%% Windowing function
w = 1./((1+exp((-t+(peak-B+offset))./A)).*(1+exp((t-(peak+B+offset))./A))); % double fermi-dirac
ws = w.*s;
end
