function v=voigt(x,pos,Gausswidth,voigtalpha)
% Unit height Voigt profile function. x is the independent variable
% (energy, wavelength, etc), Gausswidth is the Gaussian (Doppler) width
% and voigtalpha is ratio of the Gausswidth to the LorentzWidth (pressure
% width). Version 5, June, 2022
% 
% Example 1
% Calculate and plot a Voigt profile and compare to pure Gaussiaqn and pure Lorentzian
% x=1:1000;
% position=500;
% Gausswidth=50;
% alpha=1;
% V=voigt(x,position,Gausswidth,alpha);
% w=halfwidth(x,V);
% G=gaussian(x,position,w);
% L=lorentzian(x,position,w);
% plot(x,V,x,G,x,L,'linewidth',2) 
% legend('Voigt','Gaussian','Lorentzian')
%
% Example 2
% Plots five Voigt profiles with alphas from 0.001 (mostly Gaussian) to 20
% (mostly Lorentzian) and displays a table of 
% x=2:.01:4;
% figure(1)
% clf
% position=3;
% hold on
% disp('    trial     GaussW     Alpha    LorW      fwhm      area')
% for trial=1:5
%     gw=[.3 .2 .1 .05 .02]; % Gaussian width
%     al=[.001 .03 .1 3 10]; % alpha
%     v=voigt(x,position-min(x),gw(trial),al(trial));
%     plot(x,v);
%     fwhm=halfwidth(x,v);
%     va=trapz(x,v);
%     lw=gw(trial)*al(trial); % Lorentzian width
%     disp([trial gw(trial) al(trial) lw fwhm va])
% end
% hold off
% legend('alpha=.001','alpha=.03','alpha=.1','alpha=3','alpha=10')
%

LorentzWidth=Gausswidth.*voigtalpha;
xx=x-min(x);
if LorentzWidth<0, LorentzWidth=-LorentzWidth;end
if Gausswidth<0, Gausswidth=-Gausswidth;end
dx=xx(2)-xx(1); % x increment
ex=[xx-max(xx)-dx xx xx+max(xx)]; % Extended x
gau=gaussian(ex,0,Gausswidth);
lor=lorentzian(ex,pos,LorentzWidth);
VoigtConv=ifft(fft(gau).*fft(lor))./sum(lor);
g=VoigtConv./max(VoigtConv);
oex=ex-max(x)+min(x);
outrange=val2ind(oex,0):val2ind(oex,max(x))-dx;
v=g(outrange+1);
end

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);
end