
clear all

L = 80;
lattsp = 0.005;
mass = 0.5*938.92*lattsp;
ell = 2;
ellz = 0;

Gausswidth = 0.02;
Gaussdepth = 2;

Rwall = 0.18;
nwall = floor(Rwall/lattsp);
if (nwall >= L/2)
    "L is too small for Rwall"
    stop
end
numeigs = nwall/2;

r = [0:L^3-1];
nx = mod(r,L);
ny = mod((r-nx)/L,L);
nz = mod((r-ny*L-nx)/L^2,L);
x = (nx < L/2).*nx + (nx > L/2).*(nx-L);
y = (ny < L/2).*ny + (ny > L/2).*(ny-L);
z = (nz < L/2).*nz + (nz > L/2).*(nz-L);

r2 = min(nx.^2,(L-nx).^2) + min(ny.^2,(L-ny).^2) + min(nz.^2,(L-nz).^2);
rabs = sqrt(r2);
rabs(1) = 1.E-10;

nstep = 0;
dstep = 0.1;
for step = 0:dstep:nwall
    points = find(rabs >= step - 0.5*dstep & rabs < step + 0.5*dstep);
    numpoints = size(points,2);
    if (size(points,2) > 0)
        vpoints = zeros(L^3,1);
        % normalization will be fixed later
        if (ell == 0)
            vpoints(points) = 1;
        elseif (ell == 1 & ellz == 1)
            vpoints(points) = -(x(points)+i*y(points))./rabs(points);
        elseif (ell == 1 & ellz == 0)
            vpoints(points) = z(points)./rabs(points);
        elseif (ell == 1 & ellz == -1)
            vpoints(points) = (x(points)-i*y(points))./rabs(points);
        elseif (ell == 2 & ellz == 2)
            vpoints(points) = (x(points)+i*y(points)).^2./rabs(points).^2;
        elseif (ell == 2 & ellz == 1)
            vpoints(points) = ...
                -(x(points)+i*y(points)).*z(points)./rabs(points).^2;
        elseif (ell == 2 & ellz == 0)
            vpoints(points) = ...
                (2*z(points).^2-x(points).^2-y(points).^2)./rabs(points).^2;
        elseif (ell == 2 & ellz == -1)
            vpoints(points) = ...
                (x(points)-i*y(points)).*z(points)./rabs(points).^2;
        elseif (ell == 2 & ellz == -2)
            vpoints(points) = (x(points)-i*y(points)).^2./rabs(points).^2;
        else
            "ell, ellz not available"
            stop
        end
        norm = sqrt(vpoints'*vpoints);
        if (norm > 0)
            nstep = nstep + 1;  
            projectors(:,nstep) = vpoints/norm;
            dist(nstep) = step;
        end      
    end
end

V = -Gaussdepth*lattsp*exp(-0.5*r2/(Gausswidth/lattsp)^2);
Vwall = 10^6*(rabs > nwall)';

%%%%%%%%%%%%%%%

w0 = 49.D0/36.D0;
w1 = 3.D0/2.D0;
w2 = 3.D0/20.D0;
w3 = 1.D0/90.D0;

% getting the coordinates for nearest neighbors in the x,y,z 
% directions (xp is +1 in the x-direction, xm is -1 in the 
% x-direction, and so on)

r_xp = nz*L^2          + ny*L          + mod(nx+1,L);
r_xpp = nz*L^2          + ny*L          + mod(nx+2,L);
r_xppp = nz*L^2          + ny*L          + mod(nx+3,L);
r_xm = nz*L^2          + ny*L          + mod(nx-1,L);
r_xmm = nz*L^2          + ny*L          + mod(nx-2,L);
r_xmmm = nz*L^2          + ny*L          + mod(nx-3,L);

r_yp = nz*L^2          + mod(ny+1,L)*L + nx;
r_ypp = nz*L^2          + mod(ny+2,L)*L + nx;
r_yppp = nz*L^2          + mod(ny+3,L)*L + nx;
r_ym = nz*L^2          + mod(ny-1,L)*L + nx;
r_ymm = nz*L^2          + mod(ny-2,L)*L + nx;
r_ymmm = nz*L^2          + mod(ny-3,L)*L + nx;

r_zp = mod(nz+1,L)*L^2 + ny*L          + nx;
r_zpp = mod(nz+2,L)*L^2 + ny*L          + nx;
r_zppp = mod(nz+3,L)*L^2 + ny*L          + nx;
r_zm = mod(nz-1,L)*L^2 + ny*L          + nx;
r_zmm = mod(nz-2,L)*L^2 + ny*L          + nx;
r_zmmm = mod(nz-3,L)*L^2 + ny*L          + nx;

Hfree = sparse([1:L^3],[1:L^3],3.0/mass*w0);  

Hfree = Hfree + sparse([1:L^3],r_xp+1,-0.5/mass*w1); 
Hfree = Hfree - sparse([1:L^3],r_xpp+1,-0.5/mass*w2); 
Hfree = Hfree + sparse([1:L^3],r_xppp+1,-0.5/mass*w3); 
Hfree = Hfree + sparse([1:L^3],r_xm+1,-0.5/mass*w1); 
Hfree = Hfree - sparse([1:L^3],r_xmm+1,-0.5/mass*w2); 
Hfree = Hfree + sparse([1:L^3],r_xmmm+1,-0.5/mass*w3); 

Hfree = Hfree + sparse([1:L^3],r_yp+1,-0.5/mass*w1); 
Hfree = Hfree - sparse([1:L^3],r_ypp+1,-0.5/mass*w2); 
Hfree = Hfree + sparse([1:L^3],r_yppp+1,-0.5/mass*w3); 
Hfree = Hfree + sparse([1:L^3],r_ym+1,-0.5/mass*w1); 
Hfree = Hfree - sparse([1:L^3],r_ymm+1,-0.5/mass*w2); 
Hfree = Hfree + sparse([1:L^3],r_ymmm+1,-0.5/mass*w3); 

Hfree = Hfree + sparse([1:L^3],r_zp+1,-0.5/mass*w1); 
Hfree = Hfree - sparse([1:L^3],r_zpp+1,-0.5/mass*w2); 
Hfree = Hfree + sparse([1:L^3],r_zppp+1,-0.5/mass*w3); 
Hfree = Hfree + sparse([1:L^3],r_zm+1,-0.5/mass*w1); 
Hfree = Hfree - sparse([1:L^3],r_zmm+1,-0.5/mass*w2); 
Hfree = Hfree + sparse([1:L^3],r_zmmm+1,-0.5/mass*w3);

Hfree = Hfree + sparse([1:L^3],[1:L^3],Vwall); 
H = Hfree + sparse([1:L^3],[1:L^3],V);  

Hfreerad = projectors'*(Hfree*projectors);
Hrad = projectors'*(H*projectors);
Hfreerad = (Hfreerad + Hfreerad')/2;
Hrad = (Hrad + Hrad')/2;

[vfreerad,dfreerad] = eigs(Hfreerad,numeigs,'sr');
Efreerad = diag(dfreerad);
Efreerad = sort(Efreerad);

[vrad,drad] = eigs(Hrad,numeigs,'sr');
Erad = diag(drad);
Erad = sort(Erad);

format long
[Erad Efreerad]

bound = 0;

for nn = 1:numeigs
    if (Erad(nn) > 0) 
        pfreerad(nn) = sqrt(2*mass*Efreerad(nn));
        prad(nn) = sqrt(2*mass*Erad(nn)); 
        phase(nn) = ...
            mod(atan(besselj(ell+0.5,prad(nn)*nwall) ...
            /bessely(ell+0.5,prad(nn)*nwall)) ...
            - atan(besselj(ell+0.5,pfreerad(nn)*nwall) ...
            /bessely(ell+0.5,pfreerad(nn)*nwall)),pi);
    else
        bound = bound + 1;
    end
end
    
scatt_states = [bound+1:numeigs];
scatter(prad(scatt_states)'/lattsp,phase(scatt_states)*180/pi')
xlabel('momentum (MeV)')
ylabel('phase shift (deg)')    

% To check with full Hamiltonian results:

% numeigsfull = 60;
% energiesfree = eigs(Hfree,numeigsfull,'sa');
% energies = eigs(H,numeigsfull,'sa');




