function [C11, C12, C44] = fit_C11C12C44( C11C12file, C44file )

data1 = load( C11C12file );

eps_11 = data1(:,1);
sig_11 = data1(:,2); % stress component sigma_{11} in MPa
sig_22 = data1(:,3); % stress component sigma_{22} in MPa
Wn     = data1(:,4); % strain energy density in eV/A^3

data2 = load( C44file );

eps_12 = data2(:,1);
sig_12 = data2(:,2);
Ws     = data2(:,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot the lattice energy vs lattice constant
% and find the equilibrium lattice constant, a_0
[P, S] = polyfit(eps_11,sig_11,1); C11 = P(1);
[P, S] = polyfit(eps_11,sig_22,1); C12 = P(1);
[P, S] = polyfit(eps_11,Wn,2);     C11_e = P(1)*160.2e3*2;   

[P, S] = polyfit(eps_12,sig_12,1); C44 = P(1);
[P, S] = polyfit(eps_12,Ws,2);     C44_e = P(1)*160.2e3*2; 

disp(sprintf('-----------------------------------------------------'))
disp(sprintf('From Virial stress calculation,'))
disp(sprintf(['C11 = %10.8e, C12 = %10.8e, C44 = %10.8e in MPa.'],C11,C12,C44))
disp(sprintf('bulk modulus B = (C11 + 2*C12)/3 = %10.8e in MPa.',(C11+2*C12)/3))         
disp(sprintf('-----------------------------------------------------'))
disp(sprintf('From energy calculation,'))
disp(sprintf('the elastic constants C11 = %10.8e, C44 = %10.8ein MPa.',...
             C11_e, C44_e))

figure(1)
plot(eps_11, sig_11, 'b*-', eps_11, sig_22,'r*-');

figure(2)
plot(eps_12, sig_12, 'g*-')

figure(3)
plot(eps_11, Wn, 'b*-', eps_12, Ws, 'r*-')