function [a0, B] = fit_a0B( datafile )

%data = load('platt.B.dat');
data = load( datafile );

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

[P1, S1] = polyfit(a, p, 1);
[P2, S2] = polyfit(V, p, 1);

a0 = -P1(2)/P1(1);
V0 = a0^3;
B = -V0*P2(1) * 0.1; % conversion from kiloBar to GPa

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

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

figure(1)
plot(a,p,'o',aaxis,pfit,'-');