% E_calcAFM
%
% Authors: Mike Rosenbluth, Ailey Crow, Win Pin Ng 
% (Fletcher Lab, UC Berkeley)
% ---------------------------------------------------------------
%
% Updated February 18, 2013. Now able to process multiple curves. 
%
% Updated July, 2009. This curve fitting algorithm is now
% compatible with data taken on the Veeco Catalyst AFM. Enjoy.
%
% Updated August 14, 2004.  The curve fitting algorithm was
% changed so that z0 and d0 (contact point) are fit rather than
% selected by the user.  The user selects the contact point guess,
% the upper limit of the data to include, and the lower limit on
% the approach.
%
% Updated August 6, 2004.  This program finds the young's modulus
% for a given curve when given voltage vs. distance data, the
% spring constant of the cantilever (k), and the optical lever
% sensitivity (ols). 
%
%-------------------------------------------_--------------------
%
% E_calcAFM Computes elasticity from AFM extension-deflection data
%
% [stiffness] = E_calcAFM(k, Def_sens, Eguess)
% where k = cantilever stiffness in N/m, Def_sens = ols in nm/V,
% Eguess in kPa
%
% Example:
% E = E_calc(.01, 40, 0.5)

function stiffness = E_calcAFM(k, Def_sens, Eguess)
% k in N/m, Def_sens in nm/V, Eguess in kPa

% setup directories
curdir=cd;
matlabdir = uigetdir(curdir,'Select Matlab files directory');
datadir = uigetdir(curdir,'Select AFM data directory');
if (exist('datadir') == 0)
    datadir = uigetdir('C:\', 'Select Data Dir');
end

cd(datadir)

[AFMfilename, AFMpathname] = uigetfile( ...
    {'*.txt','txt files (*.txt)'}, ...
    'Pick AFM data files to analyze', ...
    'MultiSelect', 'on');

% constants
alpha = 20*(pi/180);
nu = .5;

C = 3*tan(alpha)/(4*k*1000*(1-nu.^2));

% load voltage vs. distance data
% [fname,pname] = uigetfile('*.TXT','Select AFM Data File (please)');
if iscell(AFMfilename)
    totalAFMfiles = length(AFMfilename);
else
    totalAFMfiles = 1;
end

for j = 1:totalAFMfiles
    if iscell(AFMfilename)
        fname = AFMfilename{j};
    else
        fname = AFMfilename;
    end
    pname = AFMpathname;
    
    name = strcat(pname, fname);
    datafile = importdata(name, '\t');
    cd(matlabdir)
    
    % User selects columns containing data
    distcol = 1;
    z = -datafile.data(:,distcol);
    defcol = 2;
    v = datafile.data(:,defcol);
    unitcheck = 2;
    
    if(unitcheck==1)
        d=(v.*Def_sens).*10^-3;
    elseif(unitcheck==2)
        d=v.*10^-3;
    elseif(unitcheck==3)
        d=v./k.*10^-6;
    else
        disp('User needs to learn how to count to 3.')
    end
    
    repeat = 1;
    
    while (repeat)
        clf
        figure(1)    
        plot (z, d, 'ro');
        f1 = figure(1);
        xlabel('piezo height (um)');
        ylabel('cantilever deflection (um)');
        
        % pick contact point
        title('Guess the contact point.')
        [z0, d0] = ginput(1);
        z0counterindex = find(z<z0);
        z0counter = z0counterindex(1);
        bkgndpts = 1:(z0counter-1);
        
        % do background correction m = (y1-y0)/(x1-x0)
        P = polyfit(z(bkgndpts),d(bkgndpts),1);
        slope = P(1);
        cintersect = P(2);
        
        dbkgnd = slope.*z + cintersect - d(1);
        dnew = d - dbkgnd;
        
        clf
        plot (z, dnew, 'ro');

        hold on
        plot(z,d,'k')
        xlabel('piezo height (um)');
        ylabel('cantilever deflection (um)');
        
        d = dnew;
        totalextension = z(end);
        
        title('Click on the lower limit (end of flat line) and then the upper limit (top of defl) ')
        [zrange,drange] = ginput(2);
        z1counterindex = find(z<zrange(1));
        z1counter = z1counterindex(1);
        z2counterindex = find(z<zrange(2));
        z2counter = z2counterindex(1);
        zspan = z(z0counter) - z(z2counter);
        dspan = d(z2counter) - d(z0counter);
        
        starting(1) = z0;                    % guess for z0
        starting(2) = d0;                    % guess for d0
        starting(3) = Eguess;                % guess for Young's modulus in kPa, Mike had .5
        
        zsmall = z(z1counter:z2counter);
        dsmall = d(z1counter:z2counter);
        
        
        % FIT!
        
        options=optimset('MaxIter', 1000,  'MaxFunEvals', 100000);
        
        fit_params=lsqnonlin(@hertz1, starting, [], [], options, dsmall, zsmall, C);
        
        starting(1) = fit_params(1);
        starting(2) = fit_params(2);
        starting(3) = fit_params(3);
        
        % set new range for repeated fit
        z2counterindex = find(d>(starting(2) + dspan));
        z2counter = z2counterindex(1);
        
        % set data range to curve fit
        zsmall = z(z1counter:z2counter);
        dsmall = d(z1counter:z2counter);
        
        % curve fit to find Young's Modulus
        
        options=optimset('MaxIter', 1000,  'MaxFunEvals', 100000);
        
        fit_params=lsqnonlin(@hertz1, starting, [], [], options, dsmall, zsmall, C);
        
        E = fit_params(3);
        z0fit = fit_params(1);
        d0fit = fit_params(2);
        E_actual = E;
        
        zcounternew = find(z<z0fit);
        dfit = d0fit + z0fit - z + 1./(2.*C.*E) - ((z0fit-z)./(C.*E)+ 1./(4.*C.^2.*E.^2)).^.5;
        dfit(1:zcounternew(1)) = d0fit;
        
        sse = sum((d(z1counter:z2counter) - dfit(z1counter:z2counter)).^2);
        
        clf
        plot (z, d, 'ro', z(z1counter:z2counter), dfit(z1counter:z2counter), 'b-', z0fit, d0fit, 'go');
        xlabel('piezo height (um)');
        ylabel('cantilever deflection (um)');
        title(fname);
        text(z0fit+0.2*abs(z0fit),d0fit+0.2*abs(d0fit),['Eguess = ' num2str(E) ' kPa']);
        
        button = questdlg('Good fit?','Check for good fit','Yes');
        switch button
            case 'No'
                repeat = 1;
            otherwise
                repeat = 0;
        end
        
    end  

% This ends the while loop so that the correct Efit will be recorded in the text file
    
    datasummary = [E_actual sse totalextension];

    % Write data to a textfile in the following format:
    % filename Efit(kPa) totalextension zrange(um) drange sse
    cd(curdir)
    fid = fopen('datafile.txt','a+');
    fprintf(fid,'\r\n');
    fprintf(fid,'%25s %9.3f %4.0f %14.3f %14.3f %15.3f %15.3f',...
        fname(1:(end-4)),datasummary, zspan, dspan,(zspan-dspan));
    fclose(fid);
    
    stiffness = datasummary;
    
end

fid = fopen('datafile.txt','a+');
fprintf(fid,'\r\n');
fprintf(fid,' Filename                E(kPa)   sse totalext(um)  zpiezoext(um) deflection(um) indentation(um)');
fclose(fid);


% hertz1.m
% ---------------------------------------------------------------
% This function performs a curve fit on deflection vs. piezo
% position.  It returns

function diff=hertz1(params, d_actual, z_actual, C)

Q = params(3);
z0 = params(1);
d0 = params(2);


dfit = zeros(length(z_actual),1);
for i = 1:length(z_actual)
    if z_actual(i) >= z0
        dfit(i) = d0;
    else
        dfit(i) = d0 + z0 - z_actual(i) + 1./(2.*C.*Q) - ((z0-z_actual(i))./(C.*Q)+ 1./(4.*C.^2.*Q.^2)).^.5;
    end
end

diff = (dfit - d_actual);