%%
%% 5.3 Gaussian wave packet as solution of the free Schrödinger equation (Computational example)
% A Gaussian wave packet is formed by the superposition of plane waves g_p with a Gaussian momentum distribution f(p) (see below).
% Free Schrödinger equation (FSE): -{\hbar\over i}\,{\partial\over \partial t}\psi(x,t)=
% -{\hbar^2\over 2m}\, {\partial^2\over \partial x^2}\psi(x,t)
% Wave function of the Gaussian wave packet: \psi(x,t)=\int\limits_{-\infty}^\infty f(p)g_p(x,t)dp
%

clear all
syms hbar m b x x0 X p p0 t positive
Par=[m==1 hbar==1 b==1 x0==-5 p0==5/2]

%
%% 1 Plane waves
% Plane wave as a solution of the FSE: g_p=\sqrt {1\over 2\pi\hbar }e^{iS/\hbar}, \ {\rm where} \ S=p(x-x_0)-{p^2 \over2m}t

syms g_p f S
g_p=sqrt(1/(2*pi*hbar))*exp(i*S/hbar)
S=p*(x-x0)-p^2*t/(2*m)
g_p=sube(g_p,'S'==S)

%
%% 2 Momentum distribution of the Gaussian wave packet
% Gaussian momentum distribution:f=\sqrt{b\over\hbar\sqrt\pi}e^{-b^2{(p-p_0)^2\over 2\hbar^2}

f=sqrt(b)/sqrt(hbar*sqrt(sym(pi)))*exp(-b^2*(p-p0)^2/(2*hbar^2))

%
%% 3 Integration (Fourier transform)
% \[{\psi \left( {x,t} \right) = \int\limits_{ - \infty }^\infty  {f\left( p \right){g_p}\left( {x,t} \right)dp} }\]\[ = \int\limits_{ - \infty }^\infty  {\left( {f\left( p \right){g_p}\left( {x,t} \right){e^{ - i\frac{{\left( {x - {x_0}} \right)p}}{\hbar }}}} \right){e^{iXp}}dp,\quad {\text{where }}X = \frac{{x - {x_0}}}{\hbar }} \]
% Argument

f*g_p*exp(-i*p*(x-x0)/hbar)
expand(ans)
sube(ans,[m==1]);
simplify(ans);
2*pi*ifourier(ans,p,X);
psi=sube(ans,X==(x-x0)/hbar)
psi=simplify(psi)

%
%% 4 Density
% Density: \rho=\psi^*\psi=|\psi|^2

rho=abs(psi)^2;
rho=simplify(rho,'Steps',80)

%
%% 5 Plots

t=0:2:20;
x=linspace(-10,40,200);
sube(rho,Par)
RHO=matlabFunction(ans);
plot([0 0],[0 0.6],'k')
hold on
for n=1:size(t,2)
    plot(x,RHO(t(n),x))
    txt{n}=sprintf('t=%2.f',t(n));
end
axis([-10 40 0 0.6])
xlabel('Position x')
ylabel('\rho(x,t)')
legend(txt)

%
%% 6 Heisenberg's uncertainty relation \Delta x\Delta p \ge{\hbar \over2}
% It states that the position and momentum of an action quantum cannot be measured arbitrarily exactly. This is examined here for the Gaussian wave packet.
% Mean value of position: \[{\left\langle x \right\rangle  = \int\limits_{ - \infty }^\infty  {x\rho dx \equiv \bar x} }\]
% Variance of position: \[{{{\left( {\Delta x} \right)}^2} = \left\langle {{{\left( {x - \bar x} \right)}^2}} \right\rangle  = \left\langle {{x^2}} \right\rangle  - 2\bar x\left\langle x \right\rangle  + {{\bar x}^2} = \left\langle {{x^2}} \right\rangle  - {{\bar x}^2}}\]
% _____________________________________________________
% Eigenvalue of the momentum: \[{{p_{EW}} = \left( {\hat p\psi } \right)\frac{1}{\psi } = \left( {\frac{\hbar }{i}\frac{{\partial \psi }}{{\partial x}}} \right)\frac{1}{\psi } = {p_{c0}} + {p_{c1}}x}\]
% Mean value of the momentum: \[{\bar p = \left\langle {\hat p} \right\rangle  = \int\limits_{ - \infty }^\infty  {{p_{EW}}\rho dx}  = {p_{c0}} + {p_{c1}}\bar x}\]
% Variance of the momentum: \[{\left( {\Delta p} \right)^2} = \left\langle {{{\hat p}^2}} \right\rangle  - {{\bar p}^2} = \left\langle {{{\left| {{p_{c0}} + {p_{c1}}x} \right|}^2}} \right\rangle  - {\left| {{p_{c0}} + {p_{c1}}\bar x} \right|^2}\]
%                                        \[ = \left( {{{\left| {{p_{c0}}} \right|}^2} + 2{\mathop{\rm Re}\nolimits} \left( {{p_{c0}}{p_{c1}}} \right)\bar x + {{\left| {{p_{c1}}} \right|}^2}\left\langle {{x^2}} \right\rangle } \right) - \left( {{{\left| {{p_{c0}}} \right|}^2} + 2{\mathop{\rm Re}\nolimits} \left( {{p_{c0}}{p_{c1}}} \right)\bar x + {{\left| {{p_{c1}}} \right|}^2}{{\bar x}^2}} \right)\]
%                                        \[ = {\left| {{p_{c1}}} \right|^2}{\left( {\Delta x} \right)^2}\]

syms x t positive
int(x*rho,x,-inf,inf);
x_m=simplify(ans)
int(x^2*rho,x,-inf,inf);
dx2=simplify(simplify(ans)-x_m^2)

%

p_EW=hbar/i*diff(psi,x)/psi;
p_EW=expand(simplify(p_EW))
[s,q]=coeffs(p_EW,x)
p_c0=s(2)
p_c1=s(1)
p_mean=simplify(p_c0+p_c1*x_m)
dp2=simplify(abs(p_c1)^2,'Steps',20)*dx2

%

figure
dxdp=simplify(sqrt(dx2*dp2))
sube(dxdp,Par)
fplot(ans,[0,4])
set(gcf,'DefaultTextInterpreter','Latex')
axis([0 4 0 2]);xlabel('${\rm Time}\ t$');ylabel('$\Delta x\, \Delta p\ [\hbar]$')

%
%% 7 Animation of the interference \to 'Gauss_Interference.gif'
% A beam with action quanta is split and superimposed by deflections again.
% For simplification, the measuring probability is calculated only after the deflections.
% Three cases are distinguished in the wave function \psi_{ab} (see below):
% \bullet No which-way detector (WWD): a=b=\sqrt{1/2}
% \bullet WWD does not respond: a=1,\ b=0
% \bullet WWD responds: a=0,\ b=1
%
% Wave function: \[{{\psi _{ab}}\left( {x,t} \right) = a \cdot {{\left. {\psi \left( {x,t} \right)} \right|}_{{x_0} =  - 5,{p_0} = \frac{5}{2}}} + b \cdot {{\left. {\psi \left( {x,t} \right)} \right|}_{{x_0} = 5,{p_0} =  - \frac{5}{2}}}}\]
% and \[{\left( {a,b} \right) \in \left\{ {
% \left( {\sqrt {1/2} ,\sqrt {1/2} } \right)\limits_\text{no WWD},
% \left( {1,0} \right)\limits_\text{does not respond},
% \left( {0,1} \right)\limits_\text{responds}}
%  \right\}}\]
%

psi7=sube(psi,[hbar==1,b==1]) 
psia=sube(psi7,[x0==-5 p0==5/2 ]);
psib=sube(psi7,[x0==5 p0==-5/2]);
rhoa=matlabFunction(abs(psia)^2);
rhob=matlabFunction(abs(psib)^2);
rhoab=matlabFunction(abs(psia+psib)^2/2);
XMa=20;
YMa=0.8;
x=linspace(-XMa,XMa,200);

%

figure
subplot(3,1,1)
hold on
fanimator(@(Ta)plot(x(:),rhoa(Ta,x(:))) ,'AnimationRange',[0 10],'FrameRate',25);
fanimator(@(Ta)text(0,.45,sprintf('Time t=%2.1f',Ta),'HorizontalAlignment','center'),'AnimationRange',[0 10],'FrameRate',25);
ylabel('Density \rho_a(x)')
axis([-XMa XMa 0 YMa/2])
subplot(3,1,2)
fanimator(@(Ta)plot(x(:),rhoab(Ta,x(:))) ,'AnimationRange',[0 10],'FrameRate',25);
ylabel('Density \rho_a_b(x)')
axis([-XMa XMa 0 YMa/2])
subplot(3,1,3)
fanimator(@(Ta)plot(x(:),rhob(Ta,x(:))) ,'AnimationRange',[0 10],'FrameRate',25);
xlabel('Position x')
ylabel('Density \rho_b(x)')
axis([-XMa XMa 0 YMa/2])

%

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