clear all;
filename = '10x50ff.txt';

rawdata = load(filename);
N = length(rawdata(:,1));

data =[];
for i = 1:N
    if((rawdata(i,1)>0)&&(rawdata(i,2)>0))
        data = [data; rawdata(i,:)];
    end
end

data = sortrows(data,1);
N = length(data(:,1));
plot(data(:,1),data(:,2),'.');
xlabel('Depth (nm)','fontsize',12);
ylabel('Force (uN)','fontsize',12);
title('Force-displacement curve of 10x50 fixed-fixed bridge','fontsize',12);

p = polyfit(data(:,1),data(:,2),3); 

yfit = polyval(p,data(:,1));
hold on;
plot(data(:,1),yfit,'r','linewidth',2);

function_test = @(x,a) a(1)*x.^3 + a(2)*x + a(3);
a0 = [p(1) p(3) p(4)];
dYda = {@(x,a) 1}; sgn = 0;
[a,aerr,chisq,yfit,corr] = levmar(data(:,1),data(:,2),10*ones(N,1),...
    function_test,a0,dYda,sgn);

yfit2 = function_test(data(:,1),a);
plot(data(:,1),yfit2,'g','linewidth',2);


% analysis
W = 6.48;
L = 49.44;

% From the cubic term
W = W*1e-6;
L = L*1e-6;
t = 1.96e-6;

cc = p(1)*1e21; % cubic coefficient in standard SI (N/m^3)
cc2 = a(1)*1e21;

E = 8*L^3/(pi^4*W*t)*cc;
E2 = 8*L^3/(pi^4*W*t)*cc2;
E
E2