% Demonstration of the "equalization" method for locating the perpendicular
% drop separation point for measuring the areas of two overlappting peaks.
% Requires halfwidth.m, fastsmooth.m, ExpBroaden.m, and EqualPerpDrop.m
% functions
%
format compact
clear

% Adjustable parameters
increment=.01; % Difference between adjacent x-values
x=1:increment:100; % Independent variable (time, etc)
TimeConstant=2; % Exponential broadening time constant, x units. Try .2 - 3
noise=0.0; % Try values from 0.1 or 0.001 (10% to 0.1%)
Peak1Variation=0; % fractional random variation of peak 1 amplitude
a1=2; % Pre-broadening height of peak 1 
Peak2Height=1; % Pre-broadening height of peak 2 
m1=30; % Pre-broadening position of peak 1 (Make this larger to decrease resolution)
m2=45; % Pre-broadening position of peak 2 (Make this smaller to decrease resolution)
s1=5; % Pre-broadening width (sigma) of peak 1
s2=5; % Pre-broadening width (sigma) of peak 2
SmoothRatio=.2; % Smooth width relative to peak width; increase if very noisy.
plotting=1; % 1 to plot; 0 for no plotting
% User inputs required by the EqualPerpDrop function
EstimatedPeakPositon1=32; % Estimated position of peak 1 in observed signal
EstimatedPeakPositon2=45.8; % Estimated position of peak 2 in observed signal

% Calculate simulated signal
e=2.718;
a1=a1+Peak1Variation.*randn();
a2=Peak2Height;
g1=a1.*e.^(-(x - m1).^2/(2.*s1.^2))/(sqrt(2.*pi).*s1);
g1=ExpBroaden(g1',-TimeConstant./increment)'; % Add exponential broadening
g1=g1./max(g1);
w1=halfwidth(x,g1); % Post-broadening FWHM of peak 1 (used to compute resolution)
g2=a2.*e.^(-(x - m2).^2/(2.*s2.^2))/(sqrt(2.*pi).*s2);
g2=ExpBroaden(g2',-TimeConstant./increment)'; % Add exponential broadening
g2=g2./max(g2);
w2=halfwidth(x,g2); % Post-broadening FWHM of peak 2 (used to compute resolution)
g1=g1*a1;
g2=a2.*g2;
Resolution=1.18*((m2-m1)./(w1+w2)); % Resolution definned by https://www.shimadzu.com/an/hplc/support/lib/lctalk/resol-1.html
Area1=trapz(x,g1); % True area of peak 1
Area2=trapz(x,g2);% True area of peak 2
gs=g1+g2; % Add the two peaks together
SmoothWidth=SmoothRatio*s1./increment; % Set the appropriate smooth width to minimze distortion
gs=gs+noise.*randn(size(x)); % Add white noise
gs=fastsmooth(gs,SmoothWidth,3); % Slightly smooth to get more reliable max and min points

% Call the EqualPerpDrop function
figure(1)
[EMeasArea1,EMeasArea2,valleymindex,valleymx]=EqualPerpDrop(x,gs,EstimatedPeakPositon1,EstimatedPeakPositon2,plotting);
% disp([100.*(EMeasArea1-Area1)/Area1 100.*(EMeasArea2-Area2)/Area2])
if plotting
xlabel( [' Resolution = ' num2str(Resolution) '   TimeConstant = ' num2str(TimeConstant) '     Noise = ' num2str(100.*noise/a1) '%    AreaErrors = [' num2str([100.*(EMeasArea1-Area1)/Area1 100.*(EMeasArea2-Area2)/Area2]) ']%' ])
end