function []=fcf_2mode()
%
%
%     SYNTAX
%        fcf_2mode
%
%        E0 is the energy of the electronic origin in cm-1
%        s1 is the Huang-Rhys factor for mode 1, must be greater than 0
%        hw1 is the vibration frrequency (in cm-1) for mode 1, must be greater than 0
%        s2 is the Huang-Rhys factor for mode 2, must be greater than 0
%        hw2 is the vibration frrequency (in cm-1) for mode 2, must be greater than 0
%
%
choice = menu('Select type of Spectrum','Absorption','Luminescence');
if choice <= 1
    fact=1;
else
    fact=-1;
end
%
modes= menu('How many distoring modes?','1-Mode','2-Modes');
if modes <= 1
    modes=1;
else
    modes=2;
end
%
flip=1;
dirplot={'normal','reverse'};
%
options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';
%
defans={'0','0','0','0','0'};
prompt={'Electronic Origin (in cm^{-1})', ...
        'S_{HR} for mode 1', ...
        'Frequency for mode 1 (in cm^{-1})', ...
        'S_{HR} for mode 2', ...
        'Frequency for mode 2 (in cm^{-1})'};
name='Calculation Parameters';
numlines=1;
%
%
keepgoing = -1;
while keepgoing < 0;
    
answer=inputdlg(prompt(1:2.*modes+1),name,numlines,defans(1:2.*modes+1),options);
if ~isempty(answer)
    defans(1:2.*modes+1)=answer;
end
%
E0=str2num(char(answer(1)));
s1=str2num(char(answer(2)));
hw1=str2num(char(answer(3)));
s2=0;
hw2=0;
if modes == 2
  s2=str2num(char(answer(4)));
  hw2=str2num(char(answer(5)));
end
%
I1=exp(-s1);
I2=exp(-s2);
%
test=-1;
n=0;
while test <0
    n=n+1;
    I1=[I1; I1(end).*s1./n];
    comp=I1(end)./max(I1);
    if comp < 1e-3
        test=1;
    end
end
%
test=-1;
m=0;
while (test < 0) & (modes == 2)
    m=m+1;
    I2=[I2, I2(end).*s2./m];
    comp=I2(end)./max(I2);
    if comp < 1e-3
        test=1;
    end
end
%
I12=I1*I2;
%
len1=length(I1);
len2=length(I2);
%
E1=0:1:len1-1;
E1=E1.*hw1.*fact;
E2=0:1:len2-1;
E2=E2.*hw2.*fact;
%
E12=(diag(E1)*ones(len1,len2))+(diag(E2)*ones(len2,len1))';
E12=E0+E12;
%
I12v=reshape(I12,len1.*len2,1);
E12v=reshape(E12,len1.*len2,1);
%
[E12v, Isort]=sort(E12v);
I12v=I12v(Isort);
%
[E12vu, Iu]=unique(E12v);
if length(E12vu)<length(E12v)
    fprintf(' Overlapping bars! \n')
    E12v=E12vu;
    len12=length(Iu);
    I12vu(1)=sum(I12v(1:Iu(1)));
    for j=2:len12
        I12vu(j)=sum(I12v(Iu(j-1)+1:Iu(j)));
    end
    I12v=I12vu;
end     
%
%
h1=bar(E12v./1e3,I12v,0.5);
h2=gca;
set(h2,'Xdir',char(dirplot(flip)));
xlabel('E (10^3 cm^{-1})');
ylabel('Intensity (a.u.)');
%
%
repmenu=-1;
while repmenu < 0
choice = menu('Options','Flip plot left/right','Calculate another','Exit');
if (choice <= 0) | (choice == 3)
    keepgoing=1;
    repmenu=1;
elseif choice == 1
    keepgoing=-1;
    repmenu=-1;
    flip=3-flip;
    set(h2,'Xdir',char(dirplot(flip)));
else
    keepgoing = -1;
    repmenu=1;
end
%
end
%
end
%
end



