function [vi, thresh, vrAmpl, vnDur] = find_esac(vtEod, thresh)
% James Jun 2016 Nov 6
% returns the index of the E-scan onset
% [input]
% vtEod: Time of EOD
% threshold: E-scan threshold (optional)
% [output]
% vi: index of esac

if nargin==0
    fprintf('Demo mode\n');
    S=load('20160825_t3_Ts.mat');
    find_esac(S.V20160825_t3_Ch1.times);
    return;
end

vrDI = differentiate3(differentiate3(vtEod));
vrDI = [vrDI(1); vrDI(:); vrDI(end)];
if nargin < 2
    thresh = [];
end
minDur = 2; %minimum duration of esac

if isempty(thresh)
%     thresh = -.04e-3;
    thresh = findEscanThresh(vrDI);
    fprintf('auto-thresh: %f\n', thresh);
end

vl = vrDI < thresh;
via = find(diff(vl) > 0) + 1;
vib = find(diff(vl) < 0) + 1;
vib(vib < via(1)) = []; % remove prev transitions
n = min(numel(via), numel(vib));
via = via(1:n);
vib = vib(1:n);

vnLen = vib - via;
vi = via(vnLen >= minDur);

if nargout >= 3
    viAmpl = zeros(size(vi));
    for i=1:numel(vi);
        imin = vi(i);
        while(vrDI(imin) > vrDI(imin+1))
            imin = imin + 1;
        end        
        viAmpl(i) = imin;     
    end
    vrAmpl = vrDI(viAmpl);
end

if nargout >= 4
    vnDur = zeros(size(vi));
    for i=1:numel(vi);
        imin = vi(i);
        while(vrDI(imin) < thresh)
            imin = imin + 1;
        end        
        vnDur(i) = imin - vi(i);     
    end
end

if nargout == 0
    vrX = 1:numel(vrDI);
    figure; hold on;
    plot(vrX, vrDI*1000, 'k.-');
    plot(vrX([1,end]), thresh*[1 1]*1000, 'r-', 'LineWidth', 2);
    plot(vi, vrDI(vi)*1000, 'ro');
    grid on;
    xlabel('Sample #');
    ylabel('\Delta IPI (ms)');
%     axis([0 2000 -5e-1 2e-1]);
end

end %func


function thresh = findEscanThresh(vrDI)

vrX = (0:.000005:.0005);
vrDIp = vrDI(vrDI>=0);
vrDIn = -vrDI(vrDI<=0);

vpDIp = ksdensity(vrDIp, vrX, 'function', 'survivor');
vpDIn = ksdensity(vrDIn, vrX, 'function', 'survivor');
vpDIp(1)=1;
vpDIn(1)=1;

figure; 
% suptitle(csTitle{iPhase});
AX = [];

AX(1) = subplot(311); hold on; 
plot(vrX*1000, vpDIp, 'r', vrX*1000, vpDIn, 'b');
legend('pos', 'neg');
ylabel('Survivor');
set(gca, 'YLim', [1e-4 1])
set(gca, 'YScale', 'log');
AX(end+1) = gca;
grid on; title('Finding ESAC threshold');

AX(2) = subplot(312); hold on; 
plot(vrX*1000, vpDIp./vpDIn, 'k');
ylabel('neg/pos');
set(gca, 'YLim', [0 1.2])
set(gca, 'YTick', 0:.4:1.2);
grid on;

AX(3) = subplot(313); hold on; 
vrSlope = differentiate3(vpDIp./vpDIn);
ylim = [min(vrSlope), max(vrSlope)];
[~,imin] = min(vrSlope);
thresh = vrX(imin); %in msec
plot(vrX*1000, vrSlope, 'r');
set(gca, 'YLim', ylim);
plot(vrX(imin)*[1 1]*1000, ylim, 'r-', 'LineWidth', 2);
ylabel('Slope');
grid on;
title(sprintf('threshold: %f', thresh));

linkaxes(AX, 'x');
set(AX, 'XLim', [0 .2]);
set(AX, 'XTick', 0:.05:.2);
xlabel('|\DeltaIPI| (ms)'); 

thresh = thresh * -1 ; %negative sign

end %func


function vr = differentiate3(vr)
vr = vr(3:end)-vr(1:end-2);
vr = [vr(1); vr(:); vr(end)]/2;
end %func