function [EMeasArea1,EMeasArea2,valleymindex,valleymx]=EqualPerpDrop(x,y,m1,m2,plotting)
% function [EMeasArea1,EMeasArea2,valleymindex,valleymx]=EqualPerpDrop(x,y,m1,m2,plotting)
% Function for performing the "equalization" method for locating the
% perpendicular drop separation point for measuring the areas of two
% overlapping peaks. 
%
% Inputs: x and y are the independent and dependent variables of a segment
% of signal containing two overlapping peaks. m1 and m2 are the estimated
% peak positions based on visual estimation of the signal. plotting=1 to
% plot the signal and geometrical construction; plotting=0 for no plot.
% Outputs: EMeasArea1 and EMeasArea2 are the measures areas of the two peaks.
% valleymindex and valleymx are the vector index and x value of the
% separation point.
%
% Example: see EqualPerpDropTest.m

% Calculate three straight-line segments that touche the estimated maxima
% of the two peaks. The quotient of the original signal, divided by this
% line, results in a temporarily normalized signal that has more nearly
% equal peak heights. The effect of this treatment is to deepen the valley
% between the peaks, so that it remains distinct for lower values of the
% second peak height. This is used only for the purpose of determining its
% minimum and then is discarded. Using that minimum as the separation point
% between the peaks, the perpendicular drop areas are then calculated on
% the original observed signal.
x1=m1;
x2=m2;
index1=val2ind(x,m1);
index2=val2ind(x,m2);
y1=y(index1);
y2=y(index2);
cline=(x-x1).*(y2-y1)./(x2-x1)+y1; % Two-point straight line between the two peaks.
for z=1:length(x)
    if cline(z)>y1,cline(z)=y1;end % Flat line at the height of peak 1
    if cline(z)<y2,cline(z)=y2;end % Flat line at the height of peak 2
end

my=y./cline;
index1=val2ind(x,m1);
index2=val2ind(x,m2);
valleymy=min(my(index1:index2));
valleymindex=val2ind(my,valleymy);
valleymx=x(valleymindex);

if plotting
    plot(x,y,x,cline,':',x,y./cline,[valleymx(1) valleymx(1)],[0 valleymy],'k--','LineWidth',2)
    title('Overlapping peak area measurement by geometrical height equalization')
    legend('Observed signal','cline','Equalization','valley point')
    grid
end

EMeasArea1=trapz(x(1:valleymindex),my(1:valleymindex).*cline(1:valleymindex));
EMeasArea2=trapz(x(valleymindex:index2),my(valleymindex:index2).*cline(valleymindex:index2))+trapz(x(index2:length(x)),my(index2:length(x)).*cline(index2:length(x)));

function [index,closestval]=val2ind(x,val)
% Returns the index and the value of the element of vector x that is closest to val
% If more than one element is equally close, returns vectors of indicies and values
% Tom O'Haver (toh@umd.edu) October 2006
% Examples: If x=[1 2 4 3 5 9 6 4 5 3 1], then val2ind(x,6)=7 and val2ind(x,5.1)=[5 9]
% [indices values]=val2ind(x,3.3) returns indices = [4 10] and values = [3 3]
dif=abs(x-val);
index=find((dif-min(dif))==0);
closestval=x(index);

