function [bandwidths]=getSTCSpread(IC,d_prime_mat,slope,dprime_values)
%% This is a function modified from LED_ChRmine_analysis:part03_STCspreadMJ_AD_FE.m
% This function calculates a Spread of excitation metric based on the
% d_prime_mat that is interpolated sing the MATLAB contourlines function.
%Input:
%dPrime-Mat: 32xS double matrix for the 32 electrodes used for recording
%and S comparisons used in the the dPrime analysis. This can be the
%"dPrime_cumsum" matrix generated by the calculateDPrimeMSV function. 
% slope: tonotopic slope for the analyzed animal, this is the correction
% factor with which all of the results will be multiplied. 
% dprime_values: dPrime values for which the SoE should be calculated.
% Cave: this can not be "1" for this analysis method, as through the
% interpolation the Spread would by definition be zero at that point. 
% Output: bandwidths: Matrix containg the SOE values at the specified
% dprime_values
slope_ID=slope;
bandwidths=NaN(1,size(dprime_values,2));
   
        for mm = 1:length(dprime_values)

            dprime_high  = dprime_values(mm);


            cur_data = d_prime_mat;


            f0 = figure;
            clear C1 h1 C2 h2 s1 s2
            minval_dprime1 = NaN;
            minval_dprime2 = NaN;

            try [C1, h1] = contourf(cur_data,[1 1]); s1 = getcontourlines(C1); end
            try [C2, h2] = contourf(cur_data,[dprime_high dprime_high]); s2 = getcontourlines(C2); end
            close(f0)

            f1 = figure('color', [1 1 1], 'units','normalized','position',[.1 .1 0.35 .7]); hold on;
            title(['dprime = ', num2str(dprime_high)])

            subplot(1, 2, 1); hold on;
            contourf(cur_data, [ 1:0.5:ceil(max(max(cur_data)))]);
            colormap jet;

            %                caxis(gca,[-1 max(max(max(cur_data)))])
            set(gca,'YTick', 1:1:32)
            %         set(gca,'XScale', 'log', 'XTickLabel', [])
            %         set(gca,'YTick',  data.analysis_part_02_metadata.all_frequencies([1 9 end-8 end]), 'YTickLabel',   data.analysis_part_02_metadata.all_frequencies([1 9 end-8 end])./1000)

            %     colorbar
            %     set(gca,'XTick', 1:size(data.analysis_part_02_metadata.all_intensities,1),'XTickLabel',  data.analysis_part_02_metadata.all_intensities);
            %

            %             title(['CH ',num2str(ii)])
            xlabel('Laser Intensity [%]')
            ylabel('Channel [N°]')

            if exist('s2')

                i_dprime1  = 0;
                for zz = 1:size(s1, 2);
                    if s1(zz).x(1) == s1(zz).x(end) && s1(zz).y(1) == s1(zz).y(end) % contour is closed
                    else i_dprime1 = i_dprime1 + 1;
                        dprime1(i_dprime1).x = s1(zz).x; dprime1(i_dprime1).y = s1(zz).y;
                    end
                end

                channels = [1:32]'; dprime1_by_channel = zeros(1, 32)';
                try
                    for yy = 1:size(dprime1, 2);
                        cur_thresholds = dprime1(yy);
                        for i_chan = 1:size(channels)
                            cur_channel = channels(i_chan);
                            if nnz(ismember(cur_thresholds.y, cur_channel)) > 0
                                idx_threshold = find(ismember(cur_thresholds.y, cur_channel));
                                cur_threshold = cur_thresholds.x(idx_threshold);

                                if length(cur_threshold) > 1
                                    [minval, minloc] = min(cur_threshold);
                                    cur_threshold = cur_threshold(minloc);
                                    dprime1_by_channel(cur_channel) = cur_threshold;
                                else
                                    dprime1_by_channel(cur_channel) = cur_threshold;
                                end
                            end
                        end
                    end
                catch
                end
                dprime1_by_channel(dprime1_by_channel == 0) = NaN;



                i_dprime2  = 0;
                for zz = 1:size(s2, 2);
                    if s2(zz).x(1) == s2(zz).x(end) && s2(zz).y(1) == s2(zz).y(end) % contour is closed
                    else  i_dprime2 = i_dprime2 + 1;
                        dprime2(i_dprime2).x = s2(zz).x; dprime2(i_dprime2).y = s2(zz).y;
                    end
                end

                channels = [1:32]'; dprime2_by_channel = zeros(1, 32)';


                if exist('dprime2','var')


                    for yy = 1:size(dprime2, 2);
                        cur_thresholds = dprime2(yy);
                        for i_chan = 1:size(channels)
                            cur_channel = channels(i_chan);
                            if nnz(ismember(cur_thresholds.y, cur_channel)) > 0
                                idx_threshold = find(ismember(cur_thresholds.y, cur_channel));
                                cur_threshold = cur_thresholds.x(idx_threshold);

                                if length(cur_threshold) > 1
                                    [minval, minloc] = min(cur_threshold);
                                    cur_threshold = cur_threshold(minloc);
                                    dprime2_by_channel(cur_channel) = cur_threshold;
                                else
                                    dprime2_by_channel(cur_channel) = cur_threshold;
                                end
                            end
                        end
                    end

                    dprime2_by_channel(dprime2_by_channel == 0) = NaN;

                    %                 xlim([1 32])
                    subplot(1, 2, 2); hold on;
                    plot(dprime1_by_channel, channels, 'g*--');
                    plot(dprime2_by_channel, channels, 'r*--');
                    legend('dprime = 1',['dprime = ',num2str(dprime_high)],'AutoUpdate','off','Location','SouthWest')
                    ylim([1 32]); %xlim([0 end])
                    set(gca,'YTick',[1:2:32],'YTickLabel',[1:2:32])
                    ylim([1 32]); %xlim([0 100])
                    ylabel('Electrode n°')

                    [minval_dprime1, minloc_dprime1] = min(dprime1_by_channel);
                    plot(minval_dprime1, minloc_dprime1, 'go');
                    minval_dprime2 = dprime2_by_channel(minloc_dprime1);
                else
                    text(15, 31, 'no d-prime > 2 for best electrode');
                    clear  dprime1 dprime2
                end

                if isnan(minval_dprime2)
                    text(15, 31, 'no d-prime > 2 for best electrode');
                    clear  dprime1 dprime2
                else
                    plot(minval_dprime2, minloc_dprime1, 'ro');
                    plot([minval_dprime2 minval_dprime2], [1 32], 'k--')

                    dprime1_below_dprime2 = find(dprime1_by_channel < minval_dprime2);
                    lower_electrode_in_band = min(dprime1_below_dprime2);
                    lower_electrode_off_band = lower_electrode_in_band-1;
                    upper_electrode_in_band = max(dprime1_below_dprime2);
                    upper_electrode_off_band = upper_electrode_in_band+1;

                    if  upper_electrode_off_band <= 32 && lower_electrode_off_band >= 1
                        try
                            lower_electrode_interp = interp1([dprime1_by_channel(lower_electrode_off_band) dprime1_by_channel(lower_electrode_in_band)],[lower_electrode_off_band lower_electrode_in_band],minval_dprime2);
                            upper_electrode_interp = interp1([dprime1_by_channel(upper_electrode_off_band) dprime1_by_channel(upper_electrode_in_band)],[upper_electrode_off_band upper_electrode_in_band],minval_dprime2);

                            plot([minval_dprime2 dprime1_by_channel(lower_electrode_in_band:upper_electrode_in_band)' minval_dprime2], [lower_electrode_interp lower_electrode_in_band:upper_electrode_in_band upper_electrode_interp], 'b-', 'LineWidth', 2);
                            plot(minval_dprime2,lower_electrode_interp , 'b*'); plot(minval_dprime2,upper_electrode_interp , 'b*');

                            if (upper_electrode_interp-lower_electrode_interp) > 0
                                bandwidth = (upper_electrode_interp-lower_electrode_interp).*50;% here -1
                            else
                                bandwidth = upper_electrode_interp-lower_electrode_interp.*50;
                            end

                            text(3, 31, ['Bandwidth = ', num2str(round(bandwidth)) ' \mum']);
                            bandwidths (mm) = bandwidth*slope_ID/1000;
                            clear dprime1 dprime2
                        end
                    elseif  upper_electrode_off_band <= 32 && lower_electrode_off_band < 1
                        try
                            lower_electrode_interp = interp1([dprime1_by_channel(1) dprime1_by_channel(2)],[1 2],minval_dprime2);
                            upper_electrode_interp = interp1([dprime1_by_channel(upper_electrode_off_band) dprime1_by_channel(upper_electrode_in_band)],[upper_electrode_off_band upper_electrode_in_band],minval_dprime2);

                            plot([minval_dprime2 dprime1_by_channel(lower_electrode_in_band:upper_electrode_in_band)' minval_dprime2], [lower_electrode_interp lower_electrode_in_band:upper_electrode_in_band upper_electrode_interp], 'b-', 'LineWidth', 2);
                            plot(minval_dprime2,lower_electrode_interp , 'b*'); plot(minval_dprime2,upper_electrode_interp , 'b*');

                            if (upper_electrode_interp-lower_electrode_interp) > 0
                                bandwidth = (upper_electrode_interp-1).*50;% here -1
                            else
                                bandwidth = (upper_electrode_interp-1).*50;
                            end

                            text(15, 31, ['bandwidth = ', num2str(round(bandwidth)) ' \mum']);
                            bandwidths (mm) = bandwidth*slope_ID/1000;
                            clear dprime1 dprime2
                        end
                        %                     text(3, 31, 'no lower limit detected');
                        %                      clear  dprime1 dprime2
                    elseif upper_electrode_off_band > 32 && lower_electrode_off_band >= 1
                        try
                            lower_electrode_interp = interp1([dprime1_by_channel(lower_electrode_off_band) dprime1_by_channel(lower_electrode_in_band)],[lower_electrode_off_band lower_electrode_in_band],minval_dprime2);
                            upper_electrode_interp = interp1([dprime1_by_channel(end-1) dprime1_by_channel(end)],[size(dprime1_by_channel,1)-1 size(dprime1_by_channel,1)],minval_dprime2);

                            plot([minval_dprime2 dprime1_by_channel(lower_electrode_in_band:upper_electrode_in_band)' minval_dprime2], [lower_electrode_interp lower_electrode_in_band:upper_electrode_in_band upper_electrode_interp], 'b-', 'LineWidth', 2);
                            plot(minval_dprime2,lower_electrode_interp , 'b*'); plot(minval_dprime2,upper_electrode_interp , 'b*');

                            if (upper_electrode_interp-lower_electrode_interp) > 0
                                bandwidth = (32-lower_electrode_interp).*50;%here-1
                            else
                                bandwidth = (32-lower_electrode_interp).*50;
                            end

                            text(15, 31, ['bandwidth = ', num2str(round(bandwidth)) ' \mum']);
                            bandwidths (mm) = bandwidth*slope_ID/1000;
                            clear dprime1 dprime2
                        
                        end
                        %                     text(15, 31, 'no upper limit detected');
                        %                      clear  dprime1 dprime2
                    elseif upper_electrode_off_band > 32 && lower_electrode_off_band < 1
                        text(15, 31, 'no limits detected');
                        clear  dprime1 dprime2
                    else
                        bandwidths(mm)=NaN;
                    end



                end
            else

                subplot(1, 2, 2); hold on;
                text(15, 16, ['no d-prime = ' num2str(dprime_high) ' in this file'])
                set(gca,'YTick',[1:2:32],'YTickLabel',[1:2:32])
                ylim([1 32]); %xlim([0 100])
                ylabel('Electrode n°')
                bandwidths(mm)=NaN;


            end
            close(f1)
            
        end
end
