%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Author: Rachel Cannara, Engineering Physics Department, University of
%Wisconsin-Madison, 1500 Engineering Dr., Madison, WI 53706.
%This function/script is authorized for use in government and academic
%research laboratories and non-profit institutions 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 that may arise 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.
%Original Date posted: 7/8/2004
%Latest update: 10/7/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function F_v_L = friction_v_load(file_name)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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.
%Make sure that the argument has the DI format, e.g., 'RCB2804B.480'. The output text file 
%contains a (2 x image_size) 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.
%
%Before running, move/copy the DI file(s) and following .m files to the same folder:
%friction_v_load.m, get_image_data.m, extract_num.m, di_Header_find.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Read in the DI file and separate it into three matrices (images):
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');
%if (offset_input == [])
    plot(norm)
    okay = 1;
    col1 = 'Deflection';
    offset = 0;
    load_in_volts = norm;
    while (okay == 1)
        offset = offset + input('Input offset value to be subtracted from norm (none = 0): ');
        if (offset ~= 0)
            load_in_volts = norm - offset;
            plot(load_in_volts)
            col1 = 'Load';
            okay = input('Do you want to refine the offset? (yes = 1) ');
        else
            okay = 0;
        end
    end
    %else
    %load_in_volts = norm - offset_input;
    %plot(load_in_volts)
    %col1 = 'Load';
    %end

%Multiply the load data by the normal force calibration, if it is known (otherwise
%the load units remain in volts):
norm_cal = [];
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): '); %Entering 1 or nothing maintains load units in volts.
end
if (norm_cal ~= 1)
    load = norm_cal*load_in_volts;
    units1 = '[nN]';
    pull_off_force_nN = norm_cal*min(load_in_volts)
else
    load = load_in_volts;
    units1 = '[V]';
    pull_off_force_V = min(load_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 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) ':'])
        disp('See plot for (lateral offset) image (rotate the image for a top view.)')
        okay = 1;
        while (okay == 1)
            mesh(.5*(trace + retrace))
            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)
                            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
    disp('Calculated lateral offset and halfwidth from each of the two sloped features.')
    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);
else
    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...')
        image_size
        scan_size_nm
        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)
            mesh(.5*(trace-retrace))
            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 = .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*(trace-retrace);
        col2 = ' Energy';
    else 
        diff = .5*(trace-retrace);
        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 or return): '); %1 or nothing maintains friction force in volts
    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 = [load;friction];
    plot(load,friction,'ok')
    label = [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, according to user calibration inputs:
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\r',F_v_L);
fclose(fid);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                  THE END                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%