% Developed by Matthew Brukman and Matt Marcus
% Departments of Engineering Physics and Physics
% Univerisity of Wisconsin -- Madison
% June, 2004


function Output = getfd_new_oldNanoscope(file_name)

%These commands find the line numbers in the header file that contain
%the important information
pos_spl = di_header_find(file_name,'\Samps/line');
pos_data = di_header_find(file_name,'\Data offset');
scal_data = di_header_find(file_name,'\@4:Z scale: V [Sens.');
pos_senszscan =  di_header_find(file_name,'\@Sens. Zscan');
pos_ramp = di_header_find(file_name,'Ramp size Zsweep');
pos_sensdef = di_header_find(file_name,'Sens. Deflection');

%Open the DI file, move to the various line numbers, and read the numbers
%therein to extact the values mentioned above.

fid = fopen(file_name,'r');
fseek( fid , pos_spl(2), -1 );
line = fgets(fid);
spl = extract_num(line);
line = fgets(fid);
linno = extract_num(line);
fseek(fid, pos_senszscan, -1);
line = fgets(fid);
senszscan = extract_num( line );
fseek(fid, pos_ramp, -1);
line = fgets(fid);
ramp = extract_num( line );
fseek(fid,pos_data(1),-1);
line = fgets(fid);
imag_pos = extract_num(line);   
fseek(fid,scal_data(1),-1);
line = fgets(fid);
scaling = extract_num(line);   

hscale = senszscan*ramp*2^16/spl;

% Go to 1st pixel data and start reading
fseek(fid,imag_pos,-1);
% Convert to PSD (normal force) volts
A = scaling*fread(fid,[1 2*spl] ,'int16');

%Split extend/retract data and convert to nm  
va = A(1:spl);
vr = A( (spl+1) : (2*spl) );
B = hscale*[1:spl]; 

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

% Detemine whether or not a portion of the in-contact curve is flat (ie:
% did the tip deflect so much that the laser moved beyond the detector?)
% and find the index of the end of the flat part (end_plateau)

% First smooth the data to clean up some noise (NOTE: this is only to
% determine if there's a plateau. The smoothened data is not analyzed in
% any other way).
Smooth_va = va;
smooth_range = 5;
for ii = 2:spl-1
    if ii - smooth_range < 1
        imin = 1;
        imax = ii + smooth_range;
    elseif ii + smooth_range > spl
        imin = ii - smooth_range;
        imax = spl;
    else
        imin = ii - smooth_range;
        imax = ii + smooth_range;
    end
    
    p = PolyFitPlus([B(imin:ii-1)'; B(ii+1:imax)'], [va(imin:ii-1)'; va(ii+1:imax)'], 2);

    Smooth_va(ii) = p(1)*B(ii)^2+p(2)*B(ii)+p(3);
end

va_diff = diff(Smooth_va*100)+ 100;

end_plateau = find(va_diff < 99.6,1); % end_plateau is the va index of the end of the plateau

if end_plateau > 1
    end_plateau = end_plateau + 3;
end

clear Smooth_va

plot(B,va,'b')
%%% approach curve
hold on
plot(B,vr,'r')
%%% retract curve
plot(B(end_plateau),va(end_plateau),'go')
%%% a green circle is placed at the end of the plateau

%Define pulloff voltage as the difference between the minimum PSD voltage
%and the PSD signal at zero deflection.
[minv, i] = min(vr);

zero_polyfit = PolyFitPlus(B(int32((spl+i)/2):spl), vr(int32((spl+i)/2):spl),1);                        
zero_defl_val = zero_polyfit(1)*B(i) + zero_polyfit(2);
Vpo = zero_defl_val-minv;

plot(B(i),vr(i),'ro')
%%% a red circle is placed at the pull-off pt.
hold on
plot(B(i),zero_defl_val,'ko','MarkerSize', 3)
hold on
%%% a black circle is placed at a pt. that represents out-of-contact (zero deflection)

%Fit line to the contact region of approach trace
slope_fit_end_idx = find(va(end_plateau:spl)<= (va(end_plateau)+zero_defl_val)/2,1) + end_plateau - 1;
p= PolyFitPlus( B(end_plateau:slope_fit_end_idx), va(end_plateau:slope_fit_end_idx), 1); % Fit a polynomial to the appropriate data range

plot(B(slope_fit_end_idx),va(slope_fit_end_idx),'bo') 
%%% a blue circle is placed at the snap-in pt.
hold off
pause(0.1)

st=fclose(fid);

if end_plateau > 1
    % If a plateau occurred in the data, extrapolate to find the deflection at max. extension
    Output = [-p(1) Vpo -Vpo/p(1) -(p(1)*B(1)+p(2)-zero_defl_val)/p(1)];
else
    Output = [-p(1) Vpo -Vpo/p(1) -(vr(1)-zero_defl_val)/p(1)];
end
% Thus, the output contains:
% [slope of the in-contact line; pull-off voltage; pull-off deflection (nm); max. in-contact deflection (nm)]

end

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

function [p,S,mu] = PolyFitPlus(x,y,n)
%POLYFIT Fit polynomial to data.
%This function is the same built-in Matlab polyfit function, but the errors
%have been suppressed. 

if ~isequal(size(x),size(y))
    error('MATLAB:polyfit:XYSizeMismatch',...
          'X and Y vectors must be the same size.')
end

x = x(:);
y = y(:);

if nargout > 2
   mu = [mean(x); std(x)];
   x = (x - mu(1))/mu(2);
end

% Construct Vandermonde matrix.
V(:,n+1) = ones(length(x),1,class(x));
for j = n:-1:1
   V(:,j) = x.*V(:,j+1);
end

% Solve least squares problem.
[Q,R] = qr(V,0);

p = R\(Q'*y);    % Same as p = V\y;

r = y - V*p;
p = p.';          % Polynomial coefficients are row vectors by convention.

% S is a structure containing three elements: the triangular factor from a
% QR decomposition of the Vandermonde matrix, the degrees of freedom and
% the norm of the residuals.
S.R = R;
S.df = max(0,length(y) - (n+1));
S.normr = norm(r);
end


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);
end


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

end