Fit a0B.m

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(2) plot(a,p,'o',aaxis,pfit,'-');