function friction_v_load(file_name)
clc
close all
% clear all
fclose('all')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Use this function to calibrate friction or to calculate friction versus load from a DI file
%containing three images: (1)deflection, (2)lateral trace, and (3)lateral retrace.
%The output text file contains a (image_size x 3) matrix and header indicating the proper units. This
%function also calculates and displays the pull-off force (in volts or nN), which
%can be recorded manually. If the normal and/or lateral force calibrations are
%known, input them when prompted.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Read in the DI file and separate it into three matrices (images):
[filename, pathname] = uigetfile('*.*', 'Pick all AFM files in order','*.*','MultiSelect','off');

file_name = [pathname,filename];
A = get_image_data(file_name);
image_size = length(A);
deflection = A(:,:,1);
trace = A(:,:,2);
retrace = A(:,:,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ANALYZING LOAD DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Calculate, plot, and offset the average line-by-line deflection signal to get
%load (in volts):
norm = mean(deflection');

plot(norm)
okay = 1;
col1 = 'Deflection';
offset = 0;
defl_in_volts = norm;
while (okay == 1)
    offset = offset + input('Input offset value to be subtracted from norm (none = 0): ');
    if (offset ~= 0)
        defl_in_volts = norm - offset;
        plot(defl_in_volts)
        okay = input('Do you want to refine the offset? (yes = 1) ');
    else
        okay = 0;
    end
end


%Multiply the load data by the normal force calibration, if it is known (otherwise
%the load units remain in volts):
norm_cal = 1;
cal = input('Perform lateral calibration using this image? (yes = 1) ');
if (cal == 1)
    disp('Normal and lateral forces will remain in volts for lateral calibration routine.')
    units2 = '[V]';
else
    norm_cal = input('Enter the normal force calibration in nN/V (none = 1 or return): ');
end
if (norm_cal ~= 1)
    load = norm_cal*defl_in_volts;
    units1 = '[nN]';
    pull_off_force_nN = norm_cal*min(defl_in_volts)
else
    load = defl_in_volts;
    units1 = '[V]';
    pull_off_force_V = min(defl_in_volts)
end
plot(load)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ANALYZING FRICTION DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Find the size of the image in nm:
sizepos = di_header_find(file_name,'\Scan size');
fid = fopen(file_name,'r');
fseek(fid,sizepos(1),-1);
line = fgets(fid);
scan_size_nm = extract_num(line);

%Select a region(s) over which to calculate calibration or to average line-by-line friction:
region_nm = scan_size_nm;
if (cal == 1)
    disp('Your image should contain 2 sloped features/facets parallel to the edge of the image.')
    image_size
    scan_size_nm
    edge1 = 1;
    edge2 = image_size;
    for (n = 1:2)
        disp(['Specify facet ' num2str(n) ':'])
        mesh(.5*(trace + retrace))
        view(0,90)
        disp('See plot for (lateral offset) image (rotate the image for a top view.)')
        okay = 1;
        while (okay == 1)
            disp('Input the range of pixels to include (smallest value first.)')
            query1 = round(input('Select the first edge (integer between 1 and image pixel size): '));
            if (query1 <= image_size) & (query1 >= 1)
                edge1(n) = query1;
                stillokay = 1;
                while (stillokay == 1)
                    query2 = round(input('Select the second edge (integer between 1 and image pixel size): '));
                    if (query2 <= image_size) & (query2 >= 1)
                        edge2(n) = query2;
                        region = abs(edge2(n) - edge1(n) + 1);
                        if (region > 1)
                            region_nm = region*scan_size_nm/image_size
                            resize = .5*(trace(:,edge1(n):edge2(n)) + retrace(:,edge1(n):edge2(n)));
                            mesh(resize)
                            view(0,90)
                            okay = input('Do you want to revise this region? (yes = 1) ');
                            stillokay = 0;
                        else
                            disp('Invalid region size.')
                            stillokay = 0;
                        end
                    else
                        disp('Invalid entry. Try again.')
                    end
                end
            else
                disp('Invalid entry. Try again.')
            end
        end
        half_width(n,:) = mean((.5*(trace(:,edge1(n):edge2(n)) - retrace(:,edge1(n):edge2(n))))');
        lat_offset(n,:) = mean((.5*(trace(:,edge1(n):edge2(n)) + retrace(:,edge1(n):edge2(n))))');
    end

    close all

    col2 = ' Halfwidth(1)';
    col3 = ' LatOffset(1)';
    col4 = ' Halfwidth(2)';
    col5 = ' LatOffset(2)';
    F_v_L = [load;half_width(1,:);lat_offset(1,:);half_width(2,:);lat_offset(2,:)];
    label = [col1 units1 col2 units2 col3 units2 col4 units2 col5 units2];
    k = findstr('.',file_name);
    str = file_name(k+1:k+3);
    [pathstr,name,ext,versn] = fileparts(file_name);
    new_file = fullfile(pathstr,[name '_' str '.txt']);
    fid = fopen(new_file,'w');
    fprintf(fid,'%s\r',label);
    fprintf(fid,'%5.8f %5.8f %5.8f %5.8f %5.8f\r',F_v_L);
    fclose(fid);

    quest_input = input('Do you want to calculate S (alpha/beta) now? (yes = 1) ');

    if quest_input == 1

        theta1 = 55*pi/180; % The quoted value for the left slope (slope1)
        theta2 = -55*pi/180; % The quoted value for the right slope (slope2)

        disp(sprintf('\n------------------------------------------------------------\n'));

        theta1_meas = input('Enter the angle (in degrees) of the first slope (a positive value, ideally 55 degrees): ');
        theta2_meas = input('Enter the angle (in degrees) of the second slope (a negative value, ideally -55 degrees): ');

        theta1_meas = theta1_meas*pi/180; % An uncalibrated measure of the angle of slope1
        theta2_meas = theta2_meas*pi/180; % An uncalibrated measure of the angle of slope2

        [minval1, minidx1] = min(load(1:int32(length(load)/2)));
        [minval2, minidx2] = min(load(length(load):-1:int32(length(load)/2)));
        minidx1 = minidx1 + 1;
        minidx2 = length(load)-minidx2;

        load = load(minidx1:minidx2);
        W1 = half_width(1,minidx1:minidx2);
        O1 = lat_offset(1,minidx1:minidx2);
        W2 = half_width(2,minidx1:minidx2);
        O2 = lat_offset(2,minidx1:minidx2);

        figure;
        plot(load, W1,'r.')
        hold on
        plot(load, W2,'b.')
        plot(load, O1,'g.')
        plot(load, O2,'k.')
        xlabel('Load (V or nN)')
        legend('Width of friction loop on positive slope','Width of friction loop on negative slope',...
            'Offset of friction loop on positive slope','Offset of friction loop on negative slope')
        

        p = polyfit(load,W1,1);
        dW1_fit = p(1); % the slope of W1 determined by curve fitting
        disp(['Half-width, slope 1 = ', num2str(dW1_fit)]);

        p = polyfit(load,W2,1);
        dW2_fit = p(1); % the slope of W1 determined by curve fitting
        disp(['Half-width, slope 2 = ', num2str(dW2_fit)]);

        p = polyfit(load,O1,1);
        dO1_fit = p(1); % the slope of W1 determined by curve fitting
        disp(['Offset, slope 1 = ', num2str(dO1_fit)]);

        p = polyfit(load,O2,1);
        dO2_fit = p(1); % the slope of W1 determined by curve fitting
        disp(['Offset, slope 2 = ', num2str(dO2_fit)]);

        s1_meas = tan(theta1_meas); % the uncalibrated slope1
        s2_meas = tan(theta2_meas); % the uncalibrated slope2

        gamma = ((s1_meas -s2_meas) + sqrt((s1_meas-s2_meas)^2 - 4*s1_meas*s2_meas*tan(theta1-theta2)^2))/(2*s1_meas*s2_meas*tan(theta1-theta2));

        disp(['The z-piezo correction factor gamma = ',num2str(gamma)]);

        tilt_angle = atan(gamma*s1_meas)-theta1;
        % This is the tilt angle of the surface normal with respect to the Z axis.
        % If the sample was perfectly flat, the tilt angle = 0.
        disp(['The sample tilt angle = ',num2str(tilt_angle*180/pi)]);

        theta1_true = atan(tan(theta1_meas)*gamma);
        theta2_true = atan(tan(theta2_meas)*gamma);
        % These are the true angles of the slopes (corrected)
        disp(['The measured theta1 = ',num2str(theta1_meas*180/pi)]);
        disp(['The true theta1 = ',num2str(theta1_true*180/pi)]);
        disp(['The measured theta2 = ',num2str(theta2_meas*180/pi)]);
        disp(['The true theta2 = ',num2str(theta2_true*180/pi)]);

        p = dW2_fit/dW1_fit;
        q = (dO1_fit-dO2_fit)/dW2_fit;
        
        for mu1 = 0.05:0.001:0.6
                     
            k = p*mu1/(cos(theta1_true)^2 - (mu1^2)*sin(theta1_true)^2);
            % k is the kappa from the 1996 paper
            % To solve for mu2 in terms of mu1, use the quadratic formula.
            a2 = k*sin(theta2_true)^2;
            b2 = 1;
            c2 = (-k*sin(2*theta2_true)^2)/(4*sin(theta2_true)^2);
            
            mu2 = (-b2 + sqrt(b2^2 - 4*a2*c2))/(2*a2);
            % This equation is in terms of mu1 (symbolic variable)                    
            
            gg = 0.5*((1/mu1 + mu1)*sin(2*theta1_true)*1/p - (1/mu2 + mu2)*sin(2*theta2_true)) - q;
            % This equation has mu1 and mu2 in it, so mu2 needs to be substituted with
            % mu2(mu1) from above. Also, all the other parameters (p, q,
            % theta1_true,...) have to be substituted as well.
            
            if gg < 0
                break
            end
        end
              
        mu1_found = mu1;
        disp(['The friction coefficient on slope 1 = ',num2str(mu1_found)]);
        
        mu2_found = round(mu2*1000)/1000;
        disp(['The friction coefficient on slope 2 = ',num2str(mu2_found)]);
        
        S = mu1_found/(dW1_fit*(cos(theta1_true)^2-(mu1_found^2)*sin(theta1_true)^2));
        disp(['The S (alpha/beta) conversion factor = ',num2str(S)]);       

        disp(sprintf(['\nTo convert your lateral signal T[V] to nN, here''s the equation:\n T[V]*S*k_normal[nN/nm]*defl_sens[nm/V] = T[nN]']));
    end

else
    diff = 0.5*(trace-retrace);
    query1 = input('Eliminate edges from friction trace and retrace? (yes = 1) ');
    if (query1 == 1)
        disp('Calculating friction force (or energy dissipation) from region to be specified...')
        mesh(diff)
        view(0,90)
        disp('See plot for trace-retrace image. Note: rotate the image for a top view.')
        edge1 = 1;
        edge2 = image_size;
        okay = 1;
        while (okay == 1)
            disp('Input the range of pixels to include (smallest value first.)')
            edge1 = input('Select the first edge (>=1): ');
            edge2 = input('Select the second edge (>=1): ');
            region = abs(edge2 - edge1 + 1);
            region_nm = region*scan_size_nm/image_size
            diff = 0.5*(trace(:,edge1:edge2)-retrace(:,edge1:edge2));
            col2 = ' Friction';
            mesh(diff)
            okay = input('Do you want to change the region? (yes = 1) ');
        end
    else
        disp('Calculating friction (or energy dissipation) from entire image...')
    end

    %Choose whether to calculate friction force or energy dissipation:
    query2 = input('Choose friction or energy dissipation (energy = 1): ');
    if (query2 == 1)
        diff = region_nm*diff;
        col2 = ' Energy';
    else
        col2 = ' Friction';
        query2 = 0;
    end

    %Multiply the friction data by the lateral force calibration, if it is known
    %(otherwise the friction units remain in volts or volts*nm):
    lat_cal = input('Enter the lateral calibration factor in nN/V (none = 1): ');
    if (lat_cal ~= 1) & (query2 == 1)
        friction = lat_cal*mean(diff');
        units2 = '[aJ]';
    elseif (lat_cal > 1) & (query2 == 0)
        friction = lat_cal*mean(diff');
        units2 = '[nN]';
    elseif (query2 == 1)
        friction = mean(diff');
        units2 = '[V*nm]';
    else
        friction = mean(diff');
        units2 = '[V]';
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %Combine friction and load into one (2 x image_size) matrix:
    F_v_L = [(1:length(load));load;friction];
    plot(load,friction,'ok')
    label = ['Line ' col1 units1 col2 units2];

    %Save the data into a .txt file named according to the input argument
    %(e.g., calculations for RCB2804B.480 will be saved under RCB2804B_480.txt).
    %The data will be space-delimited with a header containing the proper
    %labels and units:
    k = findstr('.',file_name);
    str = file_name(k+1:k+3);
    [pathstr,name,ext,versn] = fileparts(file_name);
    new_file = fullfile(pathstr,[name '_' str '_friction.txt']);
    fid = fopen(new_file,'w');
    fprintf(fid,'%s\r',label);
    fprintf(fid,'%5.8f %5.8f %5.8f\r',F_v_L);
    fclose(fid);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                  THE END                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 FUNCTIONS                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Author: Matthew M. Marcus, Physics Department, University of
%Wisconsin-Madison, 1150 Engineering Dr., Madison, WI 53706.
%This function/script is authorized for use in government and academic
%research laboratories and non-profit insititutions only. Though this
%function has been tested prior to its posting, it may contain mistakes or
%require improvements. In exchange for use of this free product, we 
%request that its use and any issues relating to it be reported to us. 
%Comments and suggestions are therefore welcome and should be sent to 
%Prof. Robert Carpick <carpick@engr.wisc.edu>, Engineering Physics 
%Department, UW-Madison.
%Date posted: 7/8/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Data = get_image_data(file_name)


scal_data = di_header_find(file_name,'\@2:Z scale: V [Sens.');
pos_spl = di_header_find(file_name,'\Samps');
pos_data = di_header_find(file_name,'\Data offset');
file_type_data = di_header_find(file_name,'Image Data:');


L = length(pos_data);

%image_position = ones(1,L);

fid = fopen(file_name,'r');
fseek(fid, pos_spl(1),-1);
line = fgets(fid);
spl = extract_num(line);
line = fgets(fid);
linno = extract_num(line);

for i = 1:L 
   % figure
   fseek(fid,pos_data(i),-1);
   line = fgets(fid);
   imag_pos = extract_num(line);   
   
   fseek(fid,scal_data(i),-1);
   line = fgets(fid);
   scaling = extract_num(line);   
        
   fseek(fid,imag_pos,-1);
   A = fread(fid, [spl linno],'int16');
   Data(:,:,i) = rot90(scaling*A);
   %image(Data(:,:,i));
   % fseek(fid,file_type_data(i),-1);
   % tl = fgetl(fid);
   % title(tl);
   % axis image;
end

fclose('all');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Author: Matthew M. Marcus, Physics Department, University of
%Wisconsin-Madison, 1150 Engineering Dr., Madison, WI 53706.
%This function/script is authorized for use in government and academic
%research laboratories and non-profit insititutions only. Though this
%function has been tested prior to its posting, it may contain mistakes or
%require improvements. In exchange for use of this free product, we 
%request that its use and any issues relating to it be reported to us. 
%Comments and suggestions are therefore welcome and should be sent to 
%Prof. Robert Carpick <carpick@engr.wisc.edu>, Engineering Physics 
%Department, UW-Madison.
%Date posted: 7/8/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [ position ] = DI_header_find( file_name ,  find_string);
%
% Function DI_header_find finds all occurences of "find_string" in a Digital 
% Instruments file header from "file_name"
% Example 
% [location] = DI_header_find( '033115500.001' , 'Offline Planefit')



function [position ] = di_header_find( file_name ,  find_string);

% Define End of file identifier



% Opend the file given in argument and reference as
% fid.  Also if there was an error output error
% number and error message to screen

fid = fopen(file_name,'r');
[message,errnum] = ferror(fid);
if(errnum)
   fprintf(1,'I/O Error %d \t %s',[errnum,message]);
%   break
end


header_end=0;
eof = 0;
counter = 1;


byte_location = 0;


while( and( ~eof, ~header_end ) )
   
   byte_location = ftell(fid);
   
   
   line = fgets(fid);
   
   if( (-1)==line )
      eof  = 1;
      break
   end
     
   
   if( length( findstr(line,find_string) ) )
      position(counter) = byte_location;
      counter = counter + 1;
   end
   
      
   if length( findstr( line, '\*File list end' ) )
      header_end = 1;
   end
   
      
end

fclose(fid);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Author: Matthew M. Marcus, Physics Department, University of
%Wisconsin-Madison, 1150 Engineering Dr., Madison, WI 53706.
%This function/script is authorized for use in government and academic
%research laboratories and non-profit insititutions only. Though this
%function has been tested prior to its posting, it may contain mistakes or
%require improvements. In exchange for use of this free product, we 
%request that its use and any issues relating to it be reported to us. 
%Comments and suggestions are therefore welcome and should be sent to 
%Prof. Robert Carpick <carpick@engr.wisc.edu>, Engineering Physics 
%Department, UW-Madison.
%Date posted: 7/8/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function val = extract_num(str)

%Ascii table of relevant numbers
%character    ascii code
%  e          101
%  E          69
%  0          48
%  1          49
%  2          50
%  3          51 
%  4          52
%  5          53
%  6          54 
%  7          55 
%  8          56
%  9          57

eos = 0;
R = str;

while(~eos)
   
   [T,R] = strtok(str);
   if( length(R) == 0) 
       eos = 1; 
   end
   I = find( (T>=48) & (T<=57) | 101==T | 69==T | T==173 | T== 45 | T==46 | T==40);
  
   LT = length(T);
   LI = length(I);
   
   if( LI == LT )
      J = find(T~='(');
      val = str2num(T(J));
      break
   end
      str =R;
end