%%
%% 5.4 The quantum mechanical binding problem (Computational example)
% The Heisenbergian matrix equation results from the time-independent Schrödinger equation by representing the wave function in a base.
% Orthonormal basis (OB): \[{u_k}\left( x \right) = \sqrt {\frac{2}{b}} \sin \left( {k\pi \left( {\frac{x}{b} - \frac{1}{2}} \right)} \right),\quad \left| x \right| < \frac{b}{2}\]
% Wave function: \psi(x)=\sum c_ku_k(x)
% Time-independent Schrödinger equation:\left({\hat T +\hat V}\right)\psi=E\psi
% thereby:  \int dx u_l
% \left({\hat T+\hat V}\right)\psi dx=
% \int dx u_l E\psi dx
% \ \to \
% \sum\limits_l \left({T_{kl}+V_{kl}}\right) c_l=
% Ec_k  with T_{kl}={-\hbar^2\over 2m}\int u_k u_l\prime\prime dx and V_{kl}=\int u_k \hat Vu_l\, dx
% Heisenberg's matrix equation: \left({
% (T_{kl})+(V_{kl})
% }\right)\vec c=
% E\vec c
% Woods-Saxon potential: V^{WS}(x)=-{V_0\over 1+e^{|x|-R \over a}
% Parameter: m=8,\ \hbar=1,\ b=4,\ V_0=2,\ R=1,\ a=0.1,\ \alpha=0.5,\ n\in\{5,10\}
%

clear all
syms x t real
syms R a V0 b m hbar alpha positive
syms k l integer
Par=[R==1 a==0.1 V0==2 b==4 m==8 hbar==1 alpha==0.5]
lim=double(sube(b/2,Par))
X=double(lim*linspace(-1,1));
N=5;

%
%% 1 Wood-Saxon potential
% V(x)={V_0\over e^{-{R-|x| \over a}}+1}

V=-V0/(exp(-(R-abs(x))/a)+1)

% Graphic representation

V_par=sube(V,Par)
fplot(V_par,2*[-1 1],'LineWidth',4,'Color','r')
xlabel('Position x')
ylabel('Potential V')

%
%% 2 Base
% OB: \[{u_k}\left( x \right) = \sqrt {\frac{2}{b}} \sin \left( {k\pi \left( {\frac{x}{b} - \frac{1}{2}} \right)} \right),\quad \left| x \right| < \frac{b}{2}\]
% Base as column vectors: \vec u=\left({
% u_k(x)
% }\right),\ k=1..N
% Norm matrix: N^{Norm}=\left({
% \int\limits_{-b\over 2}^{b\over 2}u_k(x)u_l(x)dx
% }\right),\ k,l=1..N

u=sqrt(2)*sin(pi*k*(x/b-1/2))/sqrt(b)

%

u_vec=subs(u,k,(1:N)')
u_vecp=sube(u_vec,Par);

% Graphic representation

fplot(u_vecp,lim*[-1,1])

% Test norm

assumeAlso(k>0 & l>0)
Norm=int(u*sube(u,k==l),x,-b/2,b/2)

%
%% 3 Kinetic energy
% T_{mat}=\left({
% T_{kl}}\right) =
% \left({ - \frac{\hbar ^2}{2m}
% \int\limits_{-\infty}^{\infty} u_ku_l^{\prime\prime} dx
% }\right)

-hbar^2/2/m*u*diff(sube(u,k==l),x,2);
int(ans,x,-b/2,b/2)
subs(ans,k,1:N);
T_mat=subs(ans,l,(1:N)')
T_mat=double(sube(T_mat,Par))

%
%% 4 Potential energy
% V_{mat}= \left( {{V_{kl}}} \right) = \left( {\int\limits_{-\infty}^{\infty} {{u_k}V{u_l}dx} } \right)


V_mat=vpaintegral(u_vecp*V_par*u_vecp',-lim,lim);
V_mat=double(V_mat)

%
%% 5 Eigenvalue problem

[EV,EW]=eigs(T_mat+V_mat)

% Get sorted eigenvectors

[EW,ind]=sort(diag(EW));
EV=EV(:,ind);

%
%% 6 Graphic representation of the wave function

figure
col=[0 0 1;0 1 0;.7 .7 0;1 0 0;0 1 1];
plot(lim*[-1 1],[0 0],':k',[0 0],[-2 1],':k')
hold on
fplot(V_par,lim*[-1 1],'LineWidth',4,'Color','r')
hnd=fplot(u_vecp'*EV(:,1:3),lim*[-1 1]);
for n=1:3, hnd(n).Color=col(n,:);end
xlabel('Position x')
title('V(x), u_k(x)')
hold off

%
%% 7 Energy spectrum

eps=0.4;
plot([0 N+1],[0 0],':k')
hold on
for n=1:N
    plot([n-eps n+eps],EW(n)*[1 1],'LineWidth',2,'Color',col(n,:))
end
xlabel('State n')
title('Energy spektrum')
hold off

%
%% 8 Animation of the probability density \to'BoundState.gif'
% Wave function: \psi(x,t)=
% \cos(\alpha) \psi^{(1)}e^{-iE_1t/\hbar}+
% \sin(\alpha) \psi^{(2)}e^{-iE_2t/\hbar}
% Density: \rho(x,t)=|\psi(x,t)|^2
% Period: T=\frac{2\pi\hbar }{E_2-E_1}

psi1=u_vecp'*EV(:,1);
psi2=u_vecp'*EV(:,2);
alph=sube(alpha,Par);
ex1=cos(alph)*exp(-1i*EW(1)*t);
ex2=sin(alph)*exp(-1i*EW(2)*t);
rho=matlabFunction(abs(ex1'*psi1+ex2'*psi2).^2);

%

T=linspace(0,2*pi/(EW(2)-EW(1)),100);
T_max=2*pi/(EW(2)-EW(1))
fanimator(@(Ta)plot(X,rho(Ta,X)) ,'AnimationRange',[0 T_max],'FrameRate',5);
axis([-lim lim 0 1.2])
xlabel('Position x')
title('Density \rho')

%

writeAnimation('BoundState.gif','Loopcount',inf)