
clear all

% The physical parameters in units of MeV raised to the appropriate power

lattsp = 0.001;
Rwall = 1.000;
ell = 2;

Gausswidth = 0.02;
Gaussdepth = 2;

mass = 0.5*938.92*lattsp;
nwall = floor(Rwall/lattsp);
numeigs = floor(nwall/4);

r = [0:nwall];
centrifugal = (ell*(ell+1))./(2*mass*max(r.^2,0.5));
V = -Gaussdepth*lattsp*exp(-0.5*r.^2/(Gausswidth/lattsp)^2);

Hfree = sparse(r+1,r+1,1.0/mass);
Hfree = Hfree + sparse(mod(r+1,nwall+1)+1,r+1,-0.5/mass);
Hfree = Hfree + sparse(mod(r-1,nwall+1)+1,r+1,-0.5/mass);
Hfree = Hfree + sparse(r+1,r+1,centrifugal);

% Must force psi(0) to be zero
Hfree(1,1) = Hfree(1,1) + 10^9;
Hfree(nwall+1,nwall+1) = Hfree(nwall+1,nwall+1) + 10^9;

H = Hfree + sparse(r+1,r+1,V);

Efree = eigs(Hfree,numeigs,'SA');
E = eigs(H,numeigs,'SA');

bound = 0;

for nn = 1:numeigs
    if (E(nn) > 0) 
        pfree(nn) = sqrt(2*mass*Efree(nn));
        p(nn) = sqrt(2*mass*E(nn)); 
        phase(nn) = ...
            mod(atan(besselj(ell+0.5,p(nn)*nwall) ...
            /bessely(ell+0.5,p(nn)*nwall)) ...
            - atan(besselj(ell+0.5,pfree(nn)*nwall) ...
            /bessely(ell+0.5,pfree(nn)*nwall)),pi);
    else
        bound = bound + 1;
    end
end

scatt_states = [bound+1:numeigs];
scatter(p(scatt_states)'/lattsp,phase(scatt_states)*180/pi')
xlabel('momentum (MeV)')
ylabel('phase shift (deg)')

