function [C11,C12] = fit_C11EB(datafile,a0,B)
data = load(datafile);
    a = data(:,1);
    E = data(:,3);
    V = data(:,1)*a0^3;
    
    [P1, S1] = polyfit((a-1), E./a0^3, 2);
    C11 = 2*P1(1) * 160.2;
    C12 = 1/2*(3*B - C11);
    disp(sprintf('C11 = %.6f (GPa)',C11));
    disp(sprintf('C12 = %.6f (GPa)',C12));
    
    aaxis = [0:0.01:1]*(max(a)-min(a))+min(a)-1;
    Efit = polyval(P1,aaxis);
    
    figure(1)
    plot(a-1,E./a0^3,'o',aaxis,Efit,'-');