
%
% construct single-particle spectrum in a cubic box and
% determine magic numbers
% 

clear

Emax = 500;                   % this is chosen arbitrarily large enough

n=(1:30);                     % quantum numbers in one cartesian direction
                              % maximum value chosen arbitrarily
                              
[x,y,z]=meshgrid(n,n,n);      % makes triplets of quantum numbers

en=(x.^2+y.^2+z.^2)*pi^2;            % single-particle energies up to a factor

nd = find(en<Emax);           % select all energies up to a Emax

energy=sort(en(nd));
nn = size(energy);
Nn = (1:1:nn(1));             % 2*Nn number of particles up to this level 
                              % including spin degenaracy

                              
Nw  = energy.^1.5/(6*pi^2) - 3*energy/(8*pi) + 3*sqrt(energy)/(4*pi) - 1/8;   
% this formula gives the total number of particles (no spin degeneracy
% here) with energies smaller than ``energy''

ENw = energy.*Nw - ...
       (energy.^2.5/(15*pi^2) - 3*energy.^2/(16*pi) +  ...
        3*energy.^1.5/(6*pi) - energy/8);          

% ENw would be a pretty good approximation to the total energy of the particles    
% evaluated as sum of single particle energies (remember, we use the independent 
% particle model) if there would be no quantum effects (no energy
% quantization)

ENn = cumsum(energy); % This is simply the total energy of the particles
                      % evaluated simply as a sum of single-particle energies 

% Nw is the so called Weyl approximation to the total number of particles
% in a cubic cavity up to a given single particle energy in the absence 
% quantization


subplot(2,1,1)

plot(2*Nn,energy,'b.', 2*Nw,energy,'r')

% notice I have included the spin degeneracy as a factor of 2

grid on
ylabel('sp energy')
xlabel('Total number of particles')
legend('Quantum/Exact','Weyl/Semiclassics',0)

subplot(2,1,2)
title('Shell correction energy')

% Shell correction energy isd the difference between the actual energy (sum
% of actual single particle energies) and what one would obtain in the 
% absence of quantization. If this energy difference is negative
% quantum/shell effects stabilaze the system. Thus this energy is the most
% negative for magic numbers


n13 = (2*Nn).^(1/3);
N13 = find(n13>3);
plot(n13,2*(ENn-ENw),'b.')

% The factor 2 is for spin degeneracy
% I plot as a function of N^{1/3} simply because for 3D systems this is
% appropriate, remember the radius of the system scales with N this way


xlabel('N_{tot}^{1/3}')
ylabel('E_{shell correction}')
grid on
