Fit a0B.m
Jump to navigation
Jump to search
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,'-');