%%
%% 5.5  The quantum mechanical scattering problem (Computational example)
% Here we calculate the reflection and transmission of a plane quantum wave at a 1D Woods-Saxon potential.
% This is done by using an approach that considers the incoming, the reflected and the transmitted wave parts.
% Schrödinger's equation is solved numerically in the potential domain.
% We obtain the coefficients for the reflection and transmission through continuous and differentiable adaptation of all wave components.
%

clear all
syms x
syms r a V0 b m hbar alpha positive
syms k E V A B R T
syms psi(x) chi(x) ua(x) us(x)
syms usl usr vsl vsr ual uar val var real
Par=[r==1 a==0.1 V0==2 b==4 m==1 hbar==1 k==sqrt(2*E)]
X0=2;
X=2*X0*linspace(-1,1);
vars=[psi(x) chi(x)];

%
%% 1 Time independent Schrödinger equation SE and potential (Woods-Saxon)
% E\psi=-{\hbar^2 \over 2m} \ {\partial^2\psi\over\partial x^2}+V\psi with V={V_0\over 1+e^{|x|-R \over a}}

e1=E*psi(x)==V*psi-hbar^2*diff(psi(x),x,2)/2/m

%

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

% Graphic  representation

VP=sube(V,Par)
fplot(VP,X0*[-1 1],'LineWidth',2,'Color','r')

%
%% 2 Approach to the wave function
% Wave function: \psi(x)=\left\{{
% e^{ikx}+Re^{-ikx},\qquad\ {\rm für}\ x<-X_0\atop
% Au_s(x)+Bu_a(x),\  {\rm für}\ |x|\le X_0\atop
% Te^{ikx}, \qquad \qquad \quad   {\rm für}\  x>X_0\
% Derivative: \frac{\partial}{\partial x}\psi(x)=\left\{{
% \frac{\partial}{\partial x}
% (e^{ikx}+Re^{-ikx}),\qquad\ {\rm für}\ x<-X_0\atop
% Av_s(x)+Bv_a(x),\qquad\ \  {\rm für}\ |x|\le X_0\atop
% \frac{\partial}{\partial x}Te^{ikx}, \qquad \qquad \qquad   {\rm für}\  x>X_0
% here means:
% \bullet u_s, u_a, v_s, v_a: symmetric and antisymmetric solutions of the Schrödinger equation
%      in the potential domain and their derivatives
% \bullet B,A,R,T: coefficients to fulfill the connection conditions
% \bullet k=\frac{\sqrt{2mE}}{\hbar}: wave number in the potential-free area

psi3=piecewise(x<=-X0,exp(1i*k*x)+R*exp(-1i*k*x),...
    abs(x)<X0,A*us(x)+B*ua(x),x>=X0,T*exp(1i*k*x))

% Derivative

psi3_x=diff(psi3,x)

%
%% 3 Numerical solution in the scattering range\Bigl(|x|\le X_0\Bigl)
% A symmetric and an antisymmetric solution span the solution space.
% Equations of motion (corresponds to SE): \[{\left\{ \begin{array}{l}
% {{u'}_{\left\{ {s,a} \right\}}} = {v_{\left\{ {s,a} \right\}}}\\
% {{v'}_{\left\{ {s,a} \right\}}} = \frac{{2m}}{{{\hbar ^2}}}\left( {V - E} \right){u_{\left\{ {s,a} \right\}}}
% \end{array} \right.}\] with  \[{\left\{ \begin{array}{l}
% {u_s}\left( 0 \right) = 1,\quad {v_s}\left( 0 \right) = 0,{\rm{ }}\left( {{\rm{symmetric}}} \right)\\
% {u_a}\left( 0 \right) = 0,\quad {v_a}\left( 0 \right) = 1,{\rm{ }}\left( {{\rm{antisymmetric}}} \right)
% \end{array} \right.}\]

isolate(e1,diff(psi(x),x,2));
[diff(psi(x),x)==chi(x) diff(chi(x),x)==rhs(ans)]
e4=sube(subs(ans),Par)

% thereby (and by the function 'NuSol' at the script end)

[t,uv_s]=NumSol(rhs(e4),2.5,linspace(0,X0,20),[1 0],vars);
plot([0 X0],[0 0],'k',t,uv_s)
[t,uv_a]=NumSol(rhs(e4),2.5,linspace(0,X0,20),[0 1],vars);
plot([0 X0],[0 0],'k',t,uv_a)

%
%% 4 Values at the connection points
% \[\left( {x =  \pm {X_0} \to usl,usr,vsl,vsr,ual,uar,val,var } \right)\]
% (u symmetric left, u symmetric right, ...)

uv=Connect(uv_s,uv_a); %function at the script end
vpa(transpose(uv),4)

%
%% 5   B A R T
% Coefficients of the wave function from continuity and differentiability at the connection points
% Continuity at x=\pm X_0: $\left\{ {\begin{array}{*{20}{c}}
%   {\mathop {\lim }\limits_{x \to  - X_0^{\left(  -  \right)}} \psi \left( x \right) = \mathop {\lim }\limits_{x \to  - X_0^{\left(  +  \right)}} \psi \left( x \right)} \\
%   {\mathop {\lim }\limits_{x \to X_0^{\left(  -  \right)}} \psi \left( x \right) = \mathop {\lim }\limits_{x \to X_0^{\left(  +  \right)}} \psi \left( x \right)}
% \end{array}} \right.$
% Differentiability at x=\pm X_0: $\left\{ {\begin{array}{*{20}{c}}
%   {\mathop {\lim }\limits_{x \to  - X_0^{\left(  -  \right)}} \frac{d}{{dx}}\psi \left( x \right) = \mathop {\lim }\limits_{x \to  - X_0^{\left(  +  \right)}} \frac{d}{{dx}}\psi \left( x \right)} \\
%   {\mathop {\lim }\limits_{x \to X_0^{\left(  -  \right)}} \frac{d}{{dx}}\psi \left( x \right) = \mathop {\lim }\limits_{x \to X_0^{\left(  +  \right)}} \frac{d}{{dx}}\psi \left( x \right)}
% \end{array}} \right.$
%

psi6=sube(psi3,Par)
psi6_x=sube(psi3_x,Par)

% Continuity left, right, differentiability right, left

BART=GetBART(psi6,psi6_x,uv,X0,10/9);
transpose(double([BART.B BART.A BART.R BART.T]))

%
%% 6 Reflection and transmission
% \bullet Reflection probability: |R(E)|^2
% \bullet Transmission probability: |T(E)|^2
% \bullet Phase of reflection: \arg(R(E))
% \bullet Phase of transmission: \arg(T(E))

EE=linspace(0,10,20);
EE=linspace(0,10,30);
TT=[];
RR=[];
for E=EE
    [t,psi_s]=NumSol(rhs(e4),E,[0 X0],[1 0],vars);
    [t,psi_a]=NumSol(rhs(e4),E,[0 X0],[0 1],vars);
    [uv]=Connect(psi_s,psi_a);
    BART=GetBART(psi6,psi6_x,uv,X0,E);
    TT=[TT double(BART.T)];
    RR=[RR double(BART.R)];
end
%Normalization
Nor=sqrt(abs(TT).^2+abs(RR).^2);
TT=TT./Nor;
RR=RR./Nor;
figure
subplot(1,2,1)
plot(EE,abs(TT).^2,'g',EE,abs(RR).^2,'r')
axis square
xlabel('Energy E')
ylabel('abs^2')
lh=legend('|T|^2','|R|^2');
subplot(1,2,2)
plot(EE,angle(TT),'g',EE,angle(RR),'R',[0 EE(end)],[0 0],':k')
axis square
xlabel('Energy E')
ylabel('phi')
lh=legend('arg(T)','arg(R)');

%
%% 7 Probability density
% Probability density: \rho(x)=|\psi(x)|^2
% Energies: E\in\left\{ {
% {1.9}\limits_{\rm Tunneling},
% {2.3}\limits_{\rm Partial\ reflektion},
% {3.25}\limits_{\rm Total\ transmission},
% {4.25}\limits_{\rm Partial\ reflection}
% } \right\}

psi8=sube(psi6,[us(x)==sym('uss') ua(x)==sym('uaa')] );
psi8_x=sube(psi6_x,[diff(us(x),x)==sym('ussx') diff(ua(x),x)==sym('uaax')] );

EE=[1.9 2.25 3.25 4.25];
for n=1:size(EE,2) %this takes some time!
    E=EE(n);
    subplot(2,2,n)
    PloRho(rhs(e4),psi8,psi8_x,E,X0,vars,sube(V,Par),sube(V0,Par))
end

%
% ________________________________________________________________________
%% Functions  for
% \cdot\ \Numerical solution in the potential area
% \cdot\ \Determination of connection conditions
% \cdot\ \Determination of the coefficients BART
% \cdot\ \Graphic  representation of the measurement density \rho

function [t,y]=NumSol(e4,E,X,X0,vars)
sube(e4,sym('E')==E);
odeFunction(ans,vars);
[t,y]=ode23s(ans,X,X0,[1 0]);
end

function uv=Connect(psi_s,psi_a)
syms usl usr vsl vsr ual uar val
uv=[...
    usl==psi_s(end,1) usr==psi_s(end,1)...
    vsl==-psi_s(end,2) vsr==psi_s(end,2) ...
    ual==-psi_a(end,1) uar==psi_a(end,1)...
    val==psi_a(end,2) sym('var')==psi_a(end,2) ...
    ];
end

function BART=GetBART(psif,psi_xf,uv,X0,E)
psi=sube(psif,sym('E')==E);
psi_x=sube(psi_xf,sym('E')==E);
[...
    limit(psi,sym('x'),-X0,'left')==sym('A')*sube(sym('usl'),uv)+sym('B')*sube(sym('ual'),uv)...
    limit(psi,sym('x'),X0,'right')==sym('A')*sube(sym('usr'),uv)+sym('B')*sube(sym('uar'),uv)...
    limit(psi_x,sym('x'),-X0,'left')==sym('A')*sube(sym('vsl'),uv)+sym('B')*sube(sym('val'),uv)...
    limit(psi_x,sym('x'),X0,'right')==sym('A')*sube(sym('vsr'),uv)+sym('B')*sube(sym('var'),uv)...
    ];
BART=solve(ans);
end

function PloRho(e4,psi6,psi6_x,E,X0,vars,V,V0)
N=14;
N=25;

[t,psis]=NumSol(e4,E,[0 X0],[1 0],vars);
[t,psia]=NumSol(e4,E,[0 X0],[0 1],vars);
uv=Connect(psis,psia);
BART=GetBART(psi6,psi6_x,uv,X0,E);
[t,psis]=NumSol(e4,E,X0*linspace(-1,1,N),[psis(end,1) -psis(end,2)],vars);
[t,psia]=NumSol(e4,E,X0*linspace(-1,1,N),[-psia(end,1) psia(end,2)],vars);
fplot(V,[-4 4],'r','LineWidth',2)
hold on
abs(BART.A*psis(:,1)+BART.B*psia(:,1)).^2;
plot(t,ans,'g')
sube(psi6,['B'==BART.B,'A'==BART.A,'R'==BART.R,'T'==BART.T]);
sube(ans,[sym('E')==E sym('uss')==0 sym('uaa')==0]);
f=abs(ans)^2;
fplot(f,[-4 -2],'g')
fplot(f,[2 4],'g')
hold off
axis([-4 4 0 4])

xlabel('Position x')
th=title(sprintf('Energy E=%3.2f',E));
lh=legend(sprintf('V_0=%0.0f',V0),'\rho');
end
