
function masses_2003(tbl,alpha,kl,lambda)

% 
%  Determine LDM parameters from N and Z >= 8 masses
% 
%  This is perhaps the most versatile program, which
% allows for various terms to be included or excluded 
% by changing the corresponding elements in the array Ind
% 
%  NB Coulomb and symmetry energy can have different forms 
%
%  tbl = 1,2,3 various choices of data sets 
%
%  alpha is the exponent of A for the pairing energy
% 
%  if kl = 0  then the pairing energy is 0 for even-even nuclei
%  if kl = 1  then the pairing energy is 0 for even-odd and odd-even nuclei
%  if kl = 2  then the pairing energy is 0 for odd-odd nuclei
%
%  lambda exponent in Wigner energy V_W exp( - lambda |N-Z|/A)
% 

if tbl == 1

%%    
%%  This set contains 1888 nuclear binding energies as provided by
%%  Pearson (via Bertsch) and used in their hfbcs-1 calculations
%%  Tondeur et al, Phys. Rev. C62, 024308 (2000) and 
%%  Samyn et al,   Nucl. Phys. A700, 1422 (2002).
%%  http://www-astro.ulb.ac.be/
%%

load mass.dat;

Z = mass(:,1);
A = mass(:,2);
N = A - Z;
I = (kl - mod(Z,2) - mod(N,2));
np   = size(N,1);
Eexp = mass(:,4);

elseif     tbl == 2 

%%   
%%  extracted binding energies from Bruxxels group 
%%  http://www-astro.ulb.ac.be/Nucdata/Masses/hfb2-plain
%%  it includes only updated (unpublished) Audi & Wapstra masses
%%  for 2135 nuclei. The calculated binding energies are hfb2 results
%%  of Goriely et al. Phys. Rev. C 66, 024326 (2002).
%%  Errors are differences between calculated and measured binding energies.
%%  In addition this file contains deformations and charge radii.
%% 

load eb_hfb2.dat;
   
N = eb_hfb2(:,1);
Z = eb_hfb2(:,2);
A = eb_hfb2(:,3);
I = (kl - mod(Z,2) - mod(N,2));
Eexp = - eb_hfb2(:,4);
Ecal = - eb_hfb2(:,5);
Err  = Eexp-Ecal;
Rch  = eb_hfb2(:,6);
b2   = eb_hfb2(:,7);
b4   = eb_hfb2(:,8);
n    = size(N);
np   = n(1);

elseif    tbl == 3
    
%%    
%%   extracted Audi & Wapstra 1995 2931 recommended masses and errors 
%%   http://csnwww.in2p3.fr/AMDC/masstables/Ame1995/mass_rmd.mas95
%%   

load mass1995.dat;

Z0   = mass1995(:,2);
A0   = mass1995(:,3);
N0   = mass1995(:,1);
Et0  = mass1995(:,4)/1000;
Eerr = mass1995(:,5)/1000;

zz = find(Z0 >= 8 & N0 >= 8);

Z = Z0(zz);
A = A0(zz);
Eexp = Et0(zz);

N = A - Z;
I = (kl - mod(Z,2) - mod(N,2));
n = size(N);
np = n(1);

end

EE = find(mod(Z,2)==0 & mod(N,2)==0);
OE = find(mod(Z,2)==1 & mod(N,2)==0);
EO = find(mod(Z,2)==0 & mod(N,2)==1);
OO = find(mod(Z,2)==1 & mod(N,2)==1);

Curv = A.^(1/3);                        
Cons = ones(size(A));                    
Cm   = 1./A.^(1/3);                      
CouL = Z.*(Z-1)./A.^(1/3);
T    = abs(N-Z)/2;                  % Isospin
SyM  = 4*T.*(T+1)./A;

Vol  = A;                               %  1

Surf = A.^(2/3);                        %  2

Coul = Z.^2./A.^(1/3);                  %  3
Coul = Z.*(Z-1)./A.^(1/3);

Couls = Coul./Surf;                     %  4

Exch = Z;                               %  5 
Exch = Z.^(4/3)./A.^(1/3);

Ccor = Z.^(2/3);                        %  6


Sym  = 4*T.^2./A;                       %  7
Sym  = 4*T.*(T+1)./A;

Syms = 4*T.^2./A.^(4/3);                %  8 
Syms = 4*T.*(T+1)./A.^(4/3);

Pair = I./A.^(alpha);                   %  9 

Wig  = exp(-lambda*abs(N-Z)./A);        %  10



Be = Eexp;

mm   = [ Vol Surf Coul Couls Exch Ccor Sym Syms Pair Wig ];
Ind =  [ 1   1    1    1     0    0    1   1    1    1   ];

Cff = find(Ind);

 M  =  eye(size(Ind,2)); 
 Ne =  zeros(size(Ind'));
for i1 = Cff
for i2 = Cff
   M(i1,i2) = mm(:,i1)'*mm(:,i2);
end
   Ne(i1, 1)= mm(:,i1)'*Be;
end
C = inv(M);
LDM = C*Ne;

 Ldm = -LDM(Cff)'

Cd = diag(C); 
zz = find(Cd==1);
Cd(zz)=0;
SqCd =  sqrt(diag(diag(Cd)));

 err =  SqCd(Cff)';
 rel =  SqCd(Cff)'./abs(LDM(Cff))';

 SCd = sqrt(diag(C));
 Ab=((1./SCd)*(1./SCd)').*C;
 Ab(Cff,Cff);

Eldm  = mm*LDM; 
Ediff = Eldm - Be;   
rms = sqrt( sum( (Ediff).^2 ) / size(Ediff,1) ) 

subplot(3,1,1)
plot(Z(EE),Ediff(EE),'ro', Z(OE),Ediff(OE),'b.',...
     Z(EO),Ediff(EO),'g.', Z(OO),Ediff(OO),'ms')
grid on
xlabel('Z')
ylabel('E(N,Z)-E_{LDM}(Z,N)')
title('Shell correction energies')
axis([0 120 -15 10])
legend('ee','oe','eo','oo',0)

subplot(3,1,2)
plot(N(EE),Ediff(EE),'ro', N(OE),Ediff(OE),'b.',...
     N(EO),Ediff(EO),'g.', N(OO),Ediff(OO),'ms')
grid on
xlabel('N')
ylabel('E(N,Z)-E_{LDM}(Z,N)')
axis([0 160 -15 10])
legend('ee','oe','eo','oo',0)

subplot(3,1,3)

hist(Ediff,50)
xlabel('E_{LDM} - E_{th} [MeV]')


