%% Convolution
% % Convolve interpolated spectrum with a gaussian
% params = load(gauss_params); % Load gaussian window
% wl_step = mean(diff(wl_crop)); % wavelength step interval
% std_dev = params.k / wl_step;
% wsize = round(3 * std_dev);
% win = -wsize:wsize;
% gauss_win = params.amp * exp(-(win).^2./(2 * std_dev.^2));
% spec_crop_convolve = conv(spec_crop,gauss_win,'same');
%% Interpolation
% Interpolate spectrum to 100x sampling rate
wl_interp100x = linspace(min(wl_crop),max(wl_crop),numel(wl_crop) * 100);
spec_crop_interp100x = interp1(wl_crop,spec_crop,wl_interp100x);
%% Integration
% Integrate interpolated spectrum over cropped wavelengths
spec_crop_sum = cumsum(spec_crop_interp100x);
%% Find threshold wavelengths
% Find max integrated power for normalization scaling
spec_max = max(spec_crop_sum);
% Normalize spectrum sum to max integrated power
spec_crop_sum_norm = spec_crop_sum / spec_max;
% Initialize wl_thresh
wl_thresh = zeros(size(thresh));
% Iterate through all power thresholds
for ii = 1:length(thresh)
    thresh_ii = thresh(ii); % extract ii-th power threshold
    % The next line extracts the index of the power value closest to the
    % power threshold
    i_p = find(spec_crop_sum_norm > thresh_ii, 1, 'first');
    % This index is used to determine the next wl_thresh
    wl_thresh(ii) = wl_interp100x(i_p);
end
%% Plot results
figure(11)
plot(wl_interp100x,spec_crop_sum_norm)
title('Auto-cropped and noise corrected cumulative sum')
hold on
xlim([min(wl_interp100x),max(wl_interp100x)])
ylim([-max(spec_crop_sum_norm)*0.05,max(spec_crop_sum_norm)*1.05])
str = {};
for jj = 1:length(thresh)
    thresh_jj = thresh(jj);
    wl_thresh_jj = wl_thresh(jj);
    plot([min(wl_interp100x),wl_thresh_jj],[thresh_jj, thresh_jj],'--g')
    plot([wl_thresh_jj,wl_thresh_jj],[-1000, thresh_jj],'--r')
    str{end + 1} = sprintf('%0.0f%% power = %0.2fnm',thresh(jj)*100,wl_thresh_jj);
end
hold off
dim = [0.7 0.5 0.3 0.3];
annotation('textbox',dim,'String',str,'FitBoxToText','on','BackgroundColor','white');