% Build Force Curves  (tested for the ascii data from the DI). 
% By Alexei Valiaev  October 2006
% Please send comments or bug reports to Alexei Valiaev at av12 at duke dot edu
% Program requires loading of the matrix with point indices that would define the
% region in which you fit the data. The output of tbe Single Curve Model
% provides you the matrix of point indices

%%%%%%%GENERAL DESCRIPTION %%%%%%%%%%%%
%The program creates the force-separation plot based on ASCII data for Photo detector and Z-piezo position
% Input parameters: ASCII data for Photo detector; ASCII data for Z-piezo position;
% Cantilever spring constant Kc  Output parameters: Separation between surface and substrate in
% nanometers [nm] force in pN (picoNewtons). Program is working in the following order. 
% Load the ASCII files from the dirname 
% User selects the slope in the contact region of the force curve manually or lets the program to select it automatically
% The zero force is defined The Photo detector signal is converted into force based on cantilever stiffness and position of zero force
% The point of zero separation is selected at point the point where the cantilever hits the solid surface.
% The point of zero separation is refined by shifting the constant Compliance region to zero.
% Force and separation are saved in the mat file. The retraction (extensional) data is being saved 
% The new file is loaded and the same steps are applied
% 
% During the file conversion user can either skip the file  and load
% the next, exit the program, proceed further

%Variables:
% dirname - directory with your files 
% convertdir - specify the directory to convert the output data
% Piezo - Piezo detenctor data in nm 
% Photo - photodector data in V
% VphtoVp - the photodetector sensiticvoty in [V/nm]
% Kc - spring constant
% Sep - course value of the separation. Then this value is being refined based on the constant complinace region 
% S - separation
% F- force 


%FUNCTIONS
% function [VphtoVp, xslopeend,xslopestart]= find_slope_retract(Photo,Piezo,Numpoints);
% used to define the value of the slope on retraction

% function[VphtoVp,xslopeend,xslopestart]=find_slope_approach(Photo,Piezo,Numpoints);
% used to define the value of the slope on approach

% [force,zero,x]=find_force(Photo,VphtoVp,numfiles,Kc,Numpoints);
%  finds the value for the zero force by taking the mean value between the approach and retraction directions in the region far from the surface

% function  S=find_sep(force,Sep,numfiles,VphtoVp,Numpoints);   
%is the second  step in finding the separation which shifts the mean value of the constant compliance region to zero separation

% function  pzero = point_zero(Photo,Piezo) - gives us the initial guess for the zero separation

function make_force_curve

clear all
dirname='D:\LAPTOP\DISK_C\DIWORK\2006\28_March\Fridge_temper\JUNK';

convertdir = 'fsep';
mkdir(dirname,convertdir)
curdir = pwd; files = dir(dirname);
numfiles = length(files) - 2; 
j = 0; 
Kc=20;
i=1;
while i<=numfiles
    filename = files(i+2).name;
    eval(['cd ',dirname]);
    Out = load(filename);
    eval(['cd ',curdir]);
    Piezo = Out(:,1);
    Photo = Out(:,2);
    Numpoints=length(Piezo)/2;
    PhotoRetract=Photo(Numpoints+1:2*Numpoints);
    PhotoApproach=Photo(1:Numpoints);
    vtonm = Out(5,3);
    Piezo = (Piezo*1000); % convert to nm scale 
    keyin = input('Use the slope on retraction (rs),Use the slope on approach(as),  ','s');
    if strcmp(keyin,'rs')
        disp('Using the RETRACTION slope')
        [VphtoVp, xslopeend,xslopestart]= find_slope_retract(Photo,Piezo,Numpoints);
    elseif strcmp(keyin,'as')
        disp('Using the APPROACH slope')
        [VphtoVp,xslopeend,xslopestart]= find_slope_approach(Photo,Piezo,Numpoints);
    elseif isempty(keyin)
        disp('Using the RETRACTION slope')
        [VphtoVp, xslopeend,xslopestart]= find_slope_retract(Photo,Piezo) ;
    end
    figure(1)
    [F,zero,x] = find_force(Photo,VphtoVp,numfiles,Kc,Numpoints);
    keyin = input('press a command Skip to next(n) Exit(e),accept this force?(y) ReSelect the force region(r) ','s');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if isempty(keyin)
        keyin = 'anythingwouldwork';
    end
    if strcmp(keyin,'n')
        close(1)
        i=i+1;
        continue;
    elseif strcmp(keyin,'e')
        break;
    elseif strcmp(keyin,'r')
        close(1)
        i=i;
        continue;
    elseif strcmp(keyin,'y')
        pzero=point_zero(Photo,Piezo);
    end
    
    if VphtoVp<0
        Sep = -(Piezo - Piezo(pzero) + (Photo(pzero) - Photo)./-VphtoVp);
        %the first part of the expression is asscociated with the of the piezo,
        %the second part is associated with the cantilever bending.
    else
        Sep = -(Piezo - Piezo(pzero) + (Photo(pzero) - Photo)./VphtoVp);
    end
    [maxpiezo,maxpiezoi] = max(Piezo);
    [maxforce,maxforcei] = max(F);
    figure(2)
    subplot(311);
    hold on
    plot(Piezo(1:Numpoints),F(1:Numpoints),'b') ;
    plot(Piezo(Numpoints+1:2*Numpoints),F(Numpoints+1:2*Numpoints),'k')
    hold off
    hold on
    plot(Piezo(xslopeend),F(xslopeend),'r*');
    plot(Piezo(xslopestart),F(xslopestart),'ro');
    hold off
    title('Force Curve');
    xlabel('(nm)');
    ylabel('Force (pN)');
    str1(1) = {['Zero = ',num2str(zero),' V']};
    str1(2) = {['K_{c} = ',num2str(Kc),' pN/nm']};
    str1(3) = {['Slope = ',num2str((1/VphtoVp(1)/vtonm)),' V_{photo}/V_{piezo}']};
    % not clear enough why 1/VphtoVp(1)/vtonm is valid
    str1(4) = {filename};
    text((maxpiezo*0.1),(maxforce*0.8),str1);
    
    S= find_sep(F,Sep,numfiles,VphtoVp,Numpoints);
    keyin = input('press a command Skip to next(n) Exit(e), ReSelect separation region(rr),accept separation region(yy)  ','s');
    %%%%%%%%%%%%%%%%%%%
    if isempty(keyin)
        keyin = 'anythingwouldwork';
    end
    %%%%%%%%%%%%%%%%%%%%%%
    if strcmp(keyin,'n')
        i=i+1;
        close(2)
        continue;
    elseif strcmp(keyin,'e')
        break;
    end
    while strcmp(keyin,'rr')
        S= find_sep(F,Sep,numfiles,VphtoVp,Numpoints);
        keyin = input('press a command Skip to next(n) Exit(e), ReSelect separation region(rr),accept separation region(yy)  ','s');
    end
    F=F(Numpoints+1:2*Numpoints);  % We will save only the retraction (extensional) part of the curve    
    
    %%save defined force and separation vectors into the Matlab file%%
    if strcmp(keyin,'yy')
        fnl=length(filename)-4;
        Filename=filename(2:fnl);
        savename=[dirname,'\',convertdir,'\','forSep',Filename];    
        Sname=['S',filename(2:fnl)];
        Fname=['F',filename(2:fnl)];
        eval([Sname,' = S;']);
        eval([Fname,' = F;']);
        eval(['save(savename,''',Sname,''',''',Fname,''')'])
        i=i+1;
    end
end

function [F,zero,x] = find_force(Photo,VphtoVp,numfiles,Kc,Numpoints);
i=Numpoints+1:2*Numpoints;
I=-i;
A=-sort(I);
PhotoRetract=Photo(A);
PhotoApproach=Photo(1:Numpoints);
figure(1)
hold on
plot(PhotoApproach,'r');
plot(PhotoRetract,'k');
hold off
hold on
disp('zoom in')
zoom on
waitforbuttonpress;  
% Initially, the list of points is empty.
xy = [];
n = 0;
disp('select the region to define zero force')
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
but = 1;
while but == 1
    [xi,yi,but] = ginput(1);
    xi=round(xi);
    yi=round(yi);
    plot(xi,PhotoApproach(xi),'*',xi,PhotoRetract(xi),'ro')
    n = n+1;
    x(:,n) = [xi];      
end;
zoom out
hold off    

photoApproach= PhotoApproach(x(1):x(2));
photoRetract= PhotoRetract(x(1):x(2));
Mean_photoApproach = mean(photoApproach);
Mean_photoRetract= mean(photoRetract);
zero= (Mean_photoApproach+Mean_photoRetract)/2; 
%%%%%%%%%%%%%%%%%%%%%%
if VphtoVp<0
    F =(1/-VphtoVp(1))*Kc*(Photo-zero);
else
    F =(1/VphtoVp(1))*Kc*(Photo-zero);
end
[maxforce,maxforcei] = max(F);
zerostr = {['Zero = ',num2str(zero),' V']}
%%%%%%%%%%

function [VphtoVp, xslopeend,xslopestart]= find_slope_retract(Photo,Piezo,Numpoints);   
figure(1);
Photo=Photo(Numpoints:2*Numpoints);
plot(Photo,'k') ;
hold on
xy = [];
n = 0;
disp('select points to define retraction slope')
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
but = 1;
while but == 1
    [xi,yi,but] = ginput(1);
    xi=round(xi);
    yi=round(yi);
    plot(xi,Photo(xi),'*')
    n = n+1;
    x(:,n) = [xi];
end
hold off    
VphtoVp = polyfit(Piezo(x(1):x(2)),Photo(x(1):x(2)),1); % fitting 1st order polynomial
VphtoVp = VphtoVp(1);
xslopestart=x(1);
xslopeend=x(2);
close(1);

function [VphtoVp,xslopeend,xslopestart]= find_slope_approach(Photo,Piezo,Numpoints)  
figure(1)
Photo=Photo(1:Numpoints);
plot(Photo,'r') 
hold on
xy = [];
n = 0;
disp('select points to define approach slope ')
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
but = 1;
while but == 1
    [xi,yi,but] = ginput(1);
    xi=round(xi);
    yi=round(yi);
    plot(xi,Photo(xi),'*')
    n = n+1;
    x(:,n) = [xi];
end;
hold off    
VphtoVp = polyfit(Piezo(x(1):x(2)),Photo(x(1):x(2)),1); 
VphtoVp = VphtoVp(1);
xslopestart=x(1);
xslopeend=x(2);
close(1)

function S= find_sep(F,Sep,numfiles,VphtoVp,Numpoints)
SepRetract=Sep(Numpoints:2*Numpoints);
if VphtoVp<0
    SepApproach=SepRetract;
else 
    SepApproach=Sep(1:Numpoints);
end

SepRetract=Sep(Numpoints+1:2*Numpoints);
forceApproach=F(1:Numpoints);
forceRetract=F(Numpoints+1:2*Numpoints);
figure(2)
subplot(312);
plot(SepApproach)
hold on
disp('zoom in')
zoom on
waitforbuttonpress;  
xy = [];
n = 0;
disp('select the region for constant compliance')
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
but = 1;
while but == 1
    [xi,yi,but] = ginput(1);
    xi=round(xi);
    yi=round(yi);
    plot(xi,SepApproach(xi),'ro')
    n = n+1;
    x(:,n) = [xi];
end;
hold off    
% Here we shift the constant compliance region to zero 
SepApproach_cut=SepApproach(x(1):x(2));   
Smean=mean(SepApproach_cut);
Smeanstr = {['Smean = ',num2str(Smean),' nm']}
Afull=Smean*ones(size(SepRetract));
S=-Afull+SepRetract;
figure(2)
subplot(313);
plot(SepRetract,forceRetract,'k',S,forceRetract,'r');


function pzero = point_zero(Photo,Piezo)
figure(1)
plot(Photo)
hold on
disp('zoom in')
zoom on
waitforbuttonpress;  
disp('Left click zero.')
[xg,yg] = ginput(1);
plot(xg,yg,'ro')
xg=round(xg);
yg=round(yg);
pzero=xg;
hold off
close(1)
