% MultipleBroadenedPeak
% Matlab/Octave function comparing area of a Gaussian peak to the area of
% the same peak after multipme steps of exponential broadening by
% convolution.
format compact
format short g
clear
tau=1; % time constant of exponential broadening 
x=0:.05:100;
xo=6; % peak position
y=exp(-(x-xo).^2); % Gaussian of unit height, peak at x=0; FWHM = 1.6651
% y=exp(-(x-xo)./tau);
b=exp(-(x-xo)./tau); % Exponential broadening function
% Optional convolution by a zero-centered Gaussian
% lx=length(x);
% hlx=round(lx/2);
% cwidth=3.5;  % Halfwidth of Gaussian broadening function
% b=[gaussian(x(1:hlx),0,cwidth) gaussian((x(hlx+1:lx)),max(x),cwidth)];   % Gaussian convolution function, c

NumCycles=50;
yc=ifft(fft(y).*fft(b)./sum(b)); % Apply broadening 1 to y
for cycle=2:NumCycles
  yc=ifft(fft(yc).*fft(b)./sum(b)); % Apply broadening 1 to yc
end

OriginalArea=trapz(x,y)  % Trapezoidal numerical integration.
BroadenedArea1=trapz(x,yc)  % Trapezoidal numerical integration.

figure(1)
clf
plot(x,y,x,yc,'linewidth',2)
ylabel('Detector signal')
xlabel(['tau = ' num2str(tau) '     Number of cycles = ' num2str(NumCycles) ])
title('Multiple successive exponential convolutions')
legend('Original peak','Broadened by fft(b)')
text(0,0.98.*max(y),['                          <=== Area = ' num2str(OriginalArea,7)])
text(x(val2ind(yc,max(yc))),max(yc),['     <== Area = ' num2str(BroadenedArea1,7)])
grid

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
%
% Example 1: 
% 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]
%
% Example 2:
% [indices values]=val2ind(x,3.3) returns indices = [4 10] and values = [3 3]
%
% Example 3:
% Find the value of x between 0 and 6 for which sin(x) is maximum.
% x=0:.001:6;
% y=sin(x);
% x0=x(val2ind(y,max(y)))
% (returns 1.5708, which as expected is pi/2).
%
% Copyright (c) 2012, Thomas C. O'Haver
% 
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to permit
% persons to whom the Software is furnished to do so, subject to the
% following conditions:
% 
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
% NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
% OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
% USE OR OTHER DEALINGS IN THE SOFTWARE.

dif=abs(x-val);
index=find((dif-min(dif))==0);
closestval=x(index);
end

