function [C44] = fit_C44EB(datafile,a0,C11n,C12n)
data = load(datafile);

    a = data(:,1);
    E = data(:,3);
    
    [P1, S1] = polyfit((a-1), E./a0^3/2, 2);
    C1111prime = 2*P1(1) * 160.2;
    
    C1111prime
    
    aaxis = [0:0.01:1]*(max(a)-min(a))+min(a)-1;
    Efit = polyval(P1,aaxis);
    
    figure(1)
    plot(a-1,E./a0^3/2,'o',aaxis,Efit,'-');
    
    syms C11 C12 C44
    
xyz0 = [ 1 0 0; % original coordinate
         0 1 0;
         0 0 1 ];
xyz1 = [ 1 1 0; % new coordinate
        -1 1 0;
         0 0 1 ];
T = zeros(3,3);

for i=1:3
    for j=1:3
        T(i,j) = dot(xyz0(i,:),xyz1(j,:)/norm(xyz1(j,:)));
    end
end

p=1; q=1; r=1; s=1;
A = T(1,p)*T(1,q)*T(1,r)*T(1,s) + T(2,p)*T(2,q)*T(2,r)*T(2,s) + ...
                                  T(3,p)*T(3,q)*T(3,r)*T(3,s);
B = T(1,p)*T(1,q)*T(2,r)*T(2,s) + T(1,p)*T(1,q)*T(3,r)*T(3,s) + ...
                                  T(2,p)*T(2,q)*T(3,r)*T(3,s);
C = T(2,p)*T(3,q)*T(2,r)*T(3,s) + T(1,p)*T(3,q)*T(1,r)*T(3,s) + ...
                                  T(1,p)*T(2,q)*T(1,r)*T(2,s);
CC = A*C11 + 2*B*C12 + 4*C*C44;
cc = subs(CC, C11, C11n); cc = subs(cc, C12, C12n);

A, B, C

disp(sprintf('C44 = %.6f (GPa)',double(solve(cc-C1111prime,C44))));