close all; clear; clc;

%%
% Atomic units
m = 1; % Electron mass
omega = 1; % Hartree per hbar (per radian)
L = 15; % Domain size
xx = linspace(-L,L,1e3); dx = xx(2)-xx(1);
V = @(x) m/2*omega.^2*x.^2;
Ns = 50; % Number of states to calculate
load('U_and_E_harmonic_N50.mat');

Psi2 = interp1(xx+6,Psi,xx,'linear',0);
F = zeros(1,size(Psi,2));
for dGi = 1:size(Psi,2)
    F(dGi) = abs(sum(Psi(:,1).*Psi2(:,dGi)).*dx);
end
F2 = sum(bsxfun(@times,Psi(:,1),Psi2))*dx;

%% Plot
cm = lines;
figure(1); clf;
a = subplot(1,2,1); hold on;
pad = 50; tau = 2;
for j = 1:30
    region = (find(V(xx-6)-E(j)-E(1)<0,1,'first')-pad):min((find(V(xx-6)-E(j)<0,1,'last')+pad),numel(xx));
    color = cm(2,:).*exp(-(j-1)/tau) + (1-exp(-(j-1)/tau)).*[1 1 1];
    patch(xx(region),0.5*Psi2(region,j)+E(j)-E(1)-5,color,'edgecolor','none','facealpha',0.2);
    plot(xx(region),0.5*Psi2(region,j)+E(j)-E(1)-5,'color',color);
    region = (find(V(xx)-E(j)-E(1)<0,1,'first')-pad):min((find(V(xx)-E(j)<0,1,'last')+pad),numel(xx));
    color = cm(1,:).*exp(-(j-1)/tau) + (1-exp(-(j-1)/tau)).*[1 1 1];
    patch(xx(region),0.5*Psi(region,j)+E(j)-E(1),color,'edgecolor','none','facealpha',0.2);
    plot(xx(region),0.5*Psi(region,j)+E(j)-E(1),'color',color);
end
plot(xx,V(xx)-E(1),'k-');
plot(xx,V(xx-6)-E(1)-5,'k-');

xlim([-4 12]); ylim([-6 6]);
xlabel('GMC'); ylabel('Energy / E_h');
grid on; box on;
axis square;

a = subplot(1,2,2); hold on;
for e = unique(sort([0.1 1:1:2]))
    F = zeros(Ns,Ns);
    Psi2 = interp1(xx+6,Psi,xx,'linear',0);
    for j = 1:Ns
        for i = 1:(Ns+1-j)
            F(i,j) = abs(sum(Psi(:,j).*Psi2(:,i+j-1)).*dx);
        end
        [~,I] = min(abs(V(xx)-(V(xx-6)+1-j)));
        L(j) = V(xx(I));
    end
    E1 = exp(-(0:Ns-1)/e); E1 = E1'./sum(E1);
    plot(1:Ns,sgolayfilt(F*E1,1,3),'.-');
end
l = plot(1:Ns,max(F(:,1))*exp(-L),'k.-');
uistack(l,'bottom');
%plot(1:Ns,FC,'.-'); hold on;
plot([1 1]*10,[eps 1],'k--');
plot([1 1]*18,[eps 1],'k--');
plot([1 1]*30,[eps 1],'k--');
xlabel('\DeltaG / E_h');
ylabel('Franck-Condon overlap');
grid on; box on;
axis square;
a.YScale = 'log';
xlim([2 40]); ylim([1e-3 1]);

%% Plot zoom of the overlaps
Psi2 = interp1(xx+6,Psi,xx,'linear',0);
F = zeros(Ns,Ns);
for j = 1:Ns
    for i = 1:(Ns+1-j)
        F(i,j) = abs(sum(Psi(:,j).*Psi2(:,i+j-1)).*dx);
    end
end
E1 = exp(-(0:Ns-1)/e); E1 = E1'./sum(E1);

cm1 = plasma(40);
purple = [0.759048870650256   0.234615862262369   0.506923752247699];
orange = [0.979440072684787   0.608377519323362   0.237494528945057];
cm2 = lines;
pad = 50; tau = 1;
figure(1); clf;
a = subplot(1,3,1); hold on;
plot(xx,V(xx)-E(1),'--','color',cm2(1,:),'LineWidth',1);
j = 1;
region = (find(V(xx)-E(1)-E(1)<0,1,'first')-pad):min((find(V(xx)-E(1)-E(1)<0,1,'last')+pad),numel(xx));
plot(xx(region),Psi(region,1)+E(1)-E(1),'-','color',cm2(1,:))
patch(xx(region),Psi(region,1)+E(1)-E(1),cm2(1,:),'edgecolor','none','facealpha',0.1);
plot(xx,V(xx-6)-8-E(1),'--','color',cm2(2,:),'LineWidth',1);
gray = [1 1 1];
for j = 1:40
    region = (find(V(xx-6)-E(j)-E(1)<0,1,'first')-pad):min((find(V(xx-6)-E(j)<0,1,'last')+pad),numel(xx));
    color = exp(-abs(9-j)/tau)*cm2(2,:)+(1-exp(-abs(9-j)/tau))*gray;
    patch(xx(region),Psi2(region,j)+E(j)-E(1)-8,color,'edgecolor','none','facealpha',0.2);
    plot(xx(region),Psi2(region,j)+E(j)-E(1)-8,'color',color);
end
grid on; box on;
xlim([-2.5 3.5]); ylim([-2.2 2.2]);
axis square;
ylabel('Energy / hbar\omega');

a = subplot(1,3,2); hold on;
plot(xx,V(xx)-E(1),'--','color',cm2(1,:),'LineWidth',1);
j = 1;
region = (find(V(xx)-E(1)-E(1)<0,1,'first')-pad):min((find(V(xx)-E(1)-E(1)<0,1,'last')+pad),numel(xx));
plot(xx(region),Psi(region,1)+E(1)-E(1),'-','color',cm2(1,:))
patch(xx(region),Psi(region,1)+E(1)-E(1),cm2(1,:),'edgecolor','none','facealpha',0.1);
plot(xx,V(xx-6)-18-E(1),'--','color',purple,'LineWidth',1);
gray = [1 1 1];
for j = 1:40
    region = (find(V(xx-6)-E(j)-E(1)<0,1,'first')-pad):min((find(V(xx-6)-E(j)<0,1,'last')+pad),numel(xx));
    color = exp(-abs(19-j)/tau)*purple+(1-exp(-abs(19-j)/tau))*gray;
    patch(xx(region),Psi2(region,j)+E(j)-E(1)-18,color,'edgecolor','none','facealpha',0.2);
    plot(xx(region),Psi2(region,j)+E(j)-E(1)-18,'color',color);
end
grid on; box on;
xlim([-2.5 3.5]); ylim([-2.2 2.2]);
axis square;
a.YTickLabel = '';
xlabel('GMC');

a = subplot(1,3,3); hold on;
plot(xx,V(xx)-E(1),'--','color',cm2(1,:),'LineWidth',1);
j = 1;
region = (find(V(xx)-E(1)-E(1)<0,1,'first')-pad):min((find(V(xx)-E(1)-E(1)<0,1,'last')+pad),numel(xx));
plot(xx(region),Psi(region,1)+E(1)-E(1),'-','color',cm2(1,:))
patch(xx(region),Psi(region,1)+E(1)-E(1),cm2(1,:),'edgecolor','none','facealpha',0.1);
plot(xx,V(xx-6)-30-E(1),'--','color',orange,'LineWidth',1);
gray = [1 1 1];
for j = 1:40
    region = (find(V(xx-6)-E(j)-E(1)<0,1,'first')-pad):min((find(V(xx-6)-E(j)<0,1,'last')+pad),numel(xx));
    color = exp(-abs(31-j)/tau)*orange+(1-exp(-abs(31-j)/tau))*gray;
    patch(xx(region),Psi2(region,j)+E(j)-E(1)-30,color,'edgecolor','none','facealpha',0.2);
    plot(xx(region),Psi2(region,j)+E(j)-E(1)-30,'color',color);
end
grid on; box on;
xlim([-2.5 3.5]); ylim([-2.2 2.2]);
axis square;
a.YTickLabel = '';