function [a0, Ecoh, B] = fit_a0EB( datafile, n, N)

data = load( datafile );

a = data(:,1);
E = data(:,3);
V = a.^3;

[P1, S1] = polyfit(a, E, 2);
[P2, S2] = polyfit(V, E/n*N, 2);

a0 = -P1(2)/(2*P1(1));
V0 = a0^3;
Ecoh = polyval(P1,a0)/n;
B = 2*V0*P2(1) * 160.2; % conversion from eV/A^3 to GPa

disp(sprintf('a0   = %.6f (A)',a0));
disp(sprintf('Ecoh = %.6f (eV)',Ecoh));
disp(sprintf('B    = %.6f (GPa)',B));

aaxis = [0:0.01:1]*(max(a)-min(a))+min(a);
Efit = polyval(P1,aaxis);
vaxis = aaxis.^3;

figure(1)
plot(a,E,'o',aaxis,Efit,'-');