% THIS FUNCTION USES THE CORRECT ERROR CALCULATION FOR THE WEIGHTED FIT 
%CORRECTION MADE ONB MAY 02 2006
function built_force_curve
close all
clear all

%dirname='D:\LAPTOP\DISK_C\DIWORK\2005\February_20th_2005\ELP180_pbs\cav_NEW_10July2006\fsep_selected';
%dirname='D:\LAPTOP\DISK_C\DIWORK\2003\av150303\test_083';
%dirname='D:\LAPTOP\DISK_C\DIWORK\2003\av150303_s3\test_009s3';
%dirname='D:\LAPTOP\DISK_C\DIWORK\2003\av150303_s3\test_038s3';
%dirname='D:\LAPTOP\DISK_C\DIWORK\2003\av150303_s3\test_025s3'; % BAD EXAMPLE 
%dirname='D:\LAPTOP\DISK_C\DIWORK\2003\av150303_s3\test_031s3';
dirname='D:\LAPTOP\DIWORK\2003\15March_2003\test_02';

p1_init =[0.25 200]; 
T=22;
point_high_force=300;
point_low_force=60;
num_draw=1;
sep_residual_info_matrix=zeros(1,2);

%%%load the file one by one into the matlab%%%
curdir = pwd;
files = dir(dirname);
numfiles = length(files) - 2;
i=1;
while i<=numfiles
    if i==1;
        p1_matrix_kuhn=zeros(1,1);
    end
    close all
    filename = files(i+2).name;
    eval(['cd ',dirname]); 
    load(filename);
    eval(['cd ',curdir]);
    endmat=length(filename)-4; 
    eval(['SepNew=SepNew',filename(7:endmat),';']);
    eval(['force=force',filename(7:endmat),';']);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%end load files%%%%%%%%%%
    %force=0.85*force;
    figure(1) 
    plot(force);
    grid on
    title('ForceVSpoints');
    xlabel('(points)');
    ylabel('Force (pN)');
    str1(1) = {['File = ',filename]};
    text((max(force)*0.2),(max(force)*0.2),str1);  
    %%%%%%%%%%%%%%%%%%%%%%
    keyin = input('Skip to next(n) Exit(e) Continue(Anything else)  ','s');
    if strcmp(keyin,'n') 
        i=i+1;
        close(1)
        continue;
    elseif strcmp(keyin,'e') %to exit the program before you've gone through all the curves
        break;
    end   
    % Find ADHESION and fitting points %%%%%%%
    xp=find_points(force);
    x_adh=xp(1);
    x_low=xp(2);
    x_high=xp(3);
    %Find contribution from the DRAG FORCE %%%%%%%%%%
    
    %%%[force,force_wdrag,Force_new,dForce,Sep_Fit] = find_drag_contrib(xp,SepNew,force);
    
    [x,high_force,low_force]=find_points_moving_window_02(force,point_high_force,point_low_force,xp(1))
    qq=1;
    Y=zeros(qq,3);
    U=zeros(qq,2);    
    G=zeros(qq,2);
    n=length(x)  
    for k=1:n/2  
        datapoints=length(force);  
        x=x+xp(1)-1; %-1 should be in this relation I checked that
        b=[x(1) x(2)];     
        Force_selected = force(b(1):b(2));
        Sep_selected = SepNew(b(1):b(2));
        Force_new=Force_selected;       
        %  [uncert_kuhn_bootstrap,uncert_contour_bootstrap,sigma_square,bb,p1_kernel]=bootstrap_uncertainity(Sep_selected,Force_new,T,p1_init,num_draw,i)
     %   [uncert_kuhn_bootstrap,uncert_contour_bootstrap,sigma_square,Vd_f,p1_kernel,p1_matrix_kuhn,sep_residual]=bootstrap_uncertainity(Sep_selected,Force_new,T,p1_init,num_draw,i,p1_matrix_kuhn);
     p1_kernel=[0.4495  132.4710];
     Y_kernel(k,:)=p1_kernel;                 
        %%%%%%%%%%%%%%Fit the data%%%%%%%%%%%%%%%%
        ff1_mod=linspace(0,-Force_new(end),length(Sep_selected))';
        ZZ1_kernel_mod = fjc(ff1_mod,p1_kernel,T);
        %%%%%%%%%%%%%%%%Calculate Separation residuals%%%%%%%%%%%%%%%%%%%
        sep_residual=calc_sep_residual_FJC_new_appr_noplot(p1_kernel,Force_new,Sep_selected);
        force_residual=calc_force_residual_no_smooth(Force_new,ZZ1_kernel_mod,ff1_mod,Sep_selected);
        %%%%%%%%%%%%%%%%%%%%%Save FJC fit results and the residuals%%%%%%
        if i==1;
            uncert_kuhn_contour_bootstrap=zeros(1,2);
            uncert_kuhn_contour_bootstrap=[uncert_kuhn_contour_bootstrap;uncert_kuhn_bootstrap uncert_contour_bootstrap];
            ZZ1_kernel_mod_matrix=zeros(length(force),1);  
            ff1_mod_matrix=zeros(length(force),1);  
            ff1_mod_matrix(:,i)=[zeros(b(1)-1,1); ff1_mod; zeros(length(SepNew)-b(2),1)];
            ZZ1_kernel_mod_matrix(:,i)=[zeros(b(1)-1,1); ZZ1_kernel_mod; zeros(length(SepNew)-b(2),1)];
            force_residual_matrix=zeros(length(force),1);    
            force_residual_matrix(:,i)=[zeros(b(1)-1,1); force_residual; zeros(length(SepNew)-b(2),1)];
            sep_residual_matrix=zeros(length(force),1);    
            sep_residual_matrix(:,i)=[zeros(b(1)-1,1); sep_residual; zeros(length(SepNew)-b(2),1)];
            force_residual_matrix(:,i)=[zeros(b(1)-1,1); force_residual; zeros(length(SepNew)-b(2),1)];
        else  
            uncert_kuhn_contour_bootstrap=[uncert_kuhn_contour_bootstrap;uncert_kuhn_bootstrap uncert_contour_bootstrap];
            sep_residual_matrix(:,i)=[zeros(b(1)-1,1); sep_residual; zeros(length(SepNew)-b(2),1)];
            ZZ1_kernel_mod_matrix(:,i)=[zeros(b(1)-1,1); ZZ1_kernel_mod; zeros(length(SepNew)-b(2),1)];
            ff1_mod_matrix(:,i)=[zeros(b(1)-1,1); ff1_mod; zeros(length(SepNew)-b(2),1)];
        end    
        %%%%%%%%%%%%%%%%%%%%%%end save fit results%%%%%%%%%%%%%%%%%%%%%
        figure(2)
        hold on
        %        plot(SepNew,-force,'k', ZZ1_kernel_mod,ff1_mod,'r','LineWidth',1);
        plot(SepNew,-force,'k');
        plot(Sep_selected,-Force_new,'k');
        plot(ZZ1_kernel_mod,ff1_mod,'r','LineWidth',1);
        hold off
        xlabel('Separation [nm]','FontSize',20,'Color',[1 0.1 0.3])
        ylabel('Force [pN]','FontSize',20,'Color',[1 0.1 0.3])    
        set(gca,'FontSize',24)
        str1(4) = {['p1_kernel = ',num2str(p1_kernel)]};
        stand_dev_fit=2*sqrt(sigma_square.*Vd_f);
        ff_vector=[0:700]';
        x_vector=zeros(701,1);
        sep_residual_smooth = smooth(sep_residual,0.4,'lowess');
        
        %Something is wrong in calculation of sep_residual_mod
        figure(3)
        hold on
        plot(-Force_new,sep_residual,'r*');
        plot(-Force_new,sep_residual_smooth,'b*');
        plot(ff_vector,x_vector,'k')
        plot(-Force_new,stand_dev_fit,'k*');
        plot(-Force_new,-stand_dev_fit,'k*');
        %plot(zz1_inverse,ff1_inverse)
        hold off
        title('sep_residual');
        xlabel('Force [pN]','FontSize',20,'Color',[1 0.1 0.3])
        ylabel('Sep residuals [nm]','FontSize',20,'Color',[1 0.1 0.3])    
        set(gca,'FontSize',24)
    end 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    keyin = input('Reselect points(r), Next figure(n),Exit(e)','s');
    if strcmp(keyin,'n') 
        i=i+1;
        close all 
        continue;
    elseif strcmp(keyin,'r')
        %%%%%%%%%%%%
        continue;
    elseif strcmp(keyin,'e') %to exit the program before you've gone through all the curves
        break;
    elseif isempty(keyin)
        keyin = 'anythingwouldwork';
        %%%%%%%%%%Save the fit parameters to the vector%%
        Short_fname_num=find_file_name(k,filename);
        if i==1
            Apar_kernel=zeros(1,3);
            points=zeros(1,2);
            Apar_kernel=[Apar_kernel;Y_kernel Short_fname_num];
            points=[points;x'];
        else 
            Apar_kernel=[Apar_kernel;Y_kernel Short_fname_num];
            points=[points;x'];
        end   
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        disp('closing figure')
        i=i+1;    
    end 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    figure(k+11)
    plot(ZZ1_kernel_mod,ff1_mod,'r',SepNew,-force,'b',Sep_selected,-Force_new,'k','LineWidth',1);
    axis([0 max(SepNew) -50 max(-force)])
end
close all;
%%%%%%%THIS PART BUILT CUSTOM DESIGNED BOX PLOTS%%%
uncert_kuhn_contour_bootstrap(1,:)=[];
Apar_kernel(1,:)=[];
points(1,:)=[];
%plot the boxplots 
[num_boostraps num_curves] = size(p1_matrix_kuhn);
Grid=repmat(1:num_curves,[num_boostraps,1]);
p1_original=Apar_kernel(:,1);
uncert_kuhn_bootstrap=uncert_kuhn_contour_bootstrap(:,1);
%%
figure(2)
hold on
a=1;
for i=1:num_curves
    conf_int_max=(p1_original(a)+2*uncert_kuhn_bootstrap(a)); 
    conf_int_min=(p1_original(a)-2*uncert_kuhn_bootstrap(a));
    % now we will draw the short lines indicating the conf. intrevals 
    line_conf_int=[a-0.2,a+0.2];
    plot(Grid(:,a),p1_matrix_kuhn(:,i),'b*',line_conf_int,conf_int_min*ones(1,2),'r',line_conf_int,conf_int_max*ones(1,2),'r');
    a=a+1;
end
hold off
set(gca,'FontSize',24);
xlabel('Curve','FontSize',20,'Color',[1 0.1 0.3])
ylabel('Kuhn Length [nm]','FontSize',20,'Color',[1 0.1 0.3])   
%%%%%%%%%%%%%%Save the data in mat files to the directory param%%%%%%%%%%%
convertdir = 'param_bootstr_weight'; 
mkdir(dirname,convertdir)
savename=[dirname,'\',convertdir,'\','p_bootstr_weight_Langev_',filename(1:(end-8))]; 
p1_matrix_kuhn_name=['p1_matrix_kuhn',filename(9:(endmat-4))];
Kernelname=['pFJC_original_pars',filename(9:(endmat-4))];
force_residual_name=['force_residual',filename(9:(endmat-4))];   
sep_residual_name=['sep_residual',filename(9:(endmat-4))];   
Pname=['points',filename(9:(endmat-4))];
ZZ1_kernel_mod_matrix_name=['FJC_fit_extension_matrix',filename(9:(endmat-4))];
ff1_mod_matrix_name=['FJC_fit_force_matrix',filename(9:(endmat-4))];
uncert_kuhn_contour_bootstrap_name=['kuhn_contour_bootstrap',filename(9:(endmat-4))];

eval([p1_matrix_kuhn_name,' = p1_matrix_kuhn;']);
eval([Kernelname,' = Apar_kernel;']);
eval([force_residual_name,' = force_residual_matrix;']);
eval([sep_residual_name,' = sep_residual_matrix;']);
eval([Pname,' = points;']);
eval([ZZ1_kernel_mod_matrix_name,' = ZZ1_kernel_mod_matrix;']);
eval([ff1_mod_matrix_name,' = ff1_mod_matrix;']);
eval([uncert_kuhn_contour_bootstrap_name,' = uncert_kuhn_contour_bootstrap;']);
eval(['save(savename,''',p1_matrix_kuhn_name,''',''',Kernelname,''',''',Pname,''',''',force_residual_name,''',''',sep_residual_name,''',''', ZZ1_kernel_mod_matrix_name,''',''',ff1_mod_matrix_name,''',''',uncert_kuhn_contour_bootstrap_name,''')'])

%%%%%%%%%%%%%%%Built the histogramms to plot Kuhn and Contour lengths%%%%%%%%
File_ind = Apar_kernel(:,3);
Apar_kernel(:,3)=[];
[p_sort,index] = sort(Apar_kernel,1);
File_ind_kuhn = File_ind(index(:,1));   %%%%%%%%%%in the followign 3 lines we found the order of the curves corresponsding to particular index.
File_ind_contour = File_ind(index(:,2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1=0.1:0.1:1;
x2=0:50:500;
kuhn_group=histc(p_sort(:,1),x1);
contour_group=histc(p_sort(:,2),x2);
%%%%%%%%%%%%%%%%saving curve_number_bin_kuhn;curve_number_bin_contour;

%Calculates uncertainities based on the bootstrap method
function [uncert_kuhn_bootstrap,uncert_contour_bootstrap,sigma_square,Vd_f,p1_kernel,p1_matrix_kuhn,sep_residual]=bootstrap_uncertainity(Sep_selected,Force_new,T,p1_init,num_draw,i,p1_matrix_kuhn);
[bootstat,bootsam] = bootstrp(num_draw,'mean',Sep_selected); % now bootstrap_test is the fake function because bootstrp rtequires soem function handle  
sep_perms=Sep_selected(bootsam);
force_perms=Force_new(bootsam);
sep_perms_sort=sort(sep_perms);
force_perms_sort=sort(force_perms);
rr=1;
p1_kernel_matrix=zeros(1,2);
p1_matrix=zeros(1,2);
sigma_square_matr=zeros(1,1);
sep_resid_weight_matr=zeros(length(Force_new),1);
%%%%%%This part will weight the original data before the bootstrap will start its job
while rr<=num_draw
    %%%%%%%%%%
    ss=1;
    delta_p1_kernel=1;
    while ss<=50 
        if  abs(delta_p1_kernel)<=0.01
            break;
        else 
            if ss==1;
                [p1_kernel_bootstr,fval,exitflag,output] = fminsearch(@fjcmin_noweight,p1_init,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),sep_perms(:,rr),-force_perms(:,rr),T);
            else
                [p1_kernel_bootstr,fval,exitflag,output] = fminsearch(@fjcmin_resid_weight,p1_init,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),sep_perms(:,rr),-force_perms(:,rr),T,sep_resid_weight,sigma_square);
                delta_p1_kernel=p1_kernel_matrix(ss)-p1_kernel_matrix(ss-1);    
            end
            p1_kernel_matrix=[p1_kernel_matrix; p1_kernel_bootstr];            
            [sep_residual,sep_resid_weight,sigma_square,sigma_square_hat,Vd_f]=calc_sep_residual_weight(p1_kernel_bootstr,force_perms(:,rr),sep_perms(:,rr));
            %  [sep_residual,sep_resid_weight,sigma_square,bb]=calc_sep_residual_weight(p1_kernel_bootstr,force_perms(:,rr),sep_perms(:,rr)); % Up to here it is OK
        end 
        ss=ss+1;
    end
    %%%%%%%%%%%%%p1_matrix will keep the values of the parameters only on the last iteration of weighting
    p1_matrix=[p1_matrix; p1_kernel_bootstr];
    % p1_matrix_kuhn(rr,:,i)=p1_kernel_bootstr(1);
    p1_matrix_kuhn(rr,i)=p1_kernel_bootstr(1);
    rr=rr+1;
end            
%%%%%%%%%%%%%%%%HERE I CLACULATE THE ORIGINAL DATA %%%

ss=1;
delta_p1_kernel=1;
while ss<=50 
    if  abs(delta_p1_kernel)<=0.01
        break;
    else 
        if ss==1;
            [p1_kernel,fval,exitflag,output] = fminsearch(@fjcmin_noweight,p1_init,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),Sep_selected,-Force_new,T);
        else
            [p1_kernel,fval,exitflag,output] = fminsearch(@fjcmin_resid_weight,p1_init,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),Sep_selected,-Force_new,T,sep_resid_weight,sigma_square);
            %             figure(2)
            %             hold on
            %             plot(-Force_new,sep_resid_weight,'*')
            %             hold off
            delta_p1_kernel=p1_kernel_matrix(ss)-p1_kernel_matrix(ss-1);    
        end
        % [sep_residual,sep_resid_weight,sigma_square,bb]=calc_sep_residual_weight(p1_kernel,Force_new,Sep_selected);    
        [sep_residual,sep_resid_weight,sigma_square,sigma_square_hat,Vd_f]=calc_sep_residual_weight(p1_kernel,Force_new,Sep_selected);
        p1_kernel_matrix=[p1_kernel_matrix; p1_kernel];     
        sep_resid_weight_matr=[sep_resid_weight_matr,sep_resid_weight];
        sigma_square_matr=[sigma_square_matr,sigma_square];
    end 
    ss=ss+1;
end

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

p1_kernel_matrix(1,:)=[];
p1_matrix(1,:)=[];
%%%%%%Lets estimate the uncertaintity using bootstrap%%%
kuhn_mean_bootstrap=sum(p1_matrix(:,1))/num_draw;
contour_mean_bootstrap=sum(p1_matrix(:,2))/num_draw;
for pp=1:num_draw
    BB_kuhn(pp)=(p1_matrix(pp,1)-kuhn_mean_bootstrap)^2;
    BB_contour(pp)=(p1_matrix(pp,2)-contour_mean_bootstrap)^2;
end 
uncert_kuhn_bootstrap=((1/(num_draw-1))*sum(BB_kuhn))^0.5;
uncert_contour_bootstrap=((1/(num_draw-1))*sum(BB_contour))^0.5; 


%%%%%%%%%%%%%%%%%%%%%%
function [sep_residual,sep_resid_weight,sigma_square,sigma_square_hat,Vd_f]=calc_sep_residual_weight(p1_kernel,Force_new,Sep_selected);
%function [sep_residual,sep_resid_weight,sigma_square,bb]=calc_sep_residual_weight(p1_kernel,Force_new,Sep_selected);
%Force_new=-Force_new;
T=22;
KT = 4.1*((T+273.15)/298);
num_point=1;
sep_residual=zeros(1,1);
Index_matrix=zeros(1,1);
p1=p1_kernel;
%Problem is in the calculation of sep_residual_element
while num_point < length(Force_new);
    f=-Force_new(num_point);
    z = p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1)));
    sep_residual_element=Sep_selected(num_point)-z; 
    sep_residual=[sep_residual;sep_residual_element]; 
    num_point=num_point+1;
end

b=p1(1)*(-Force_new)./KT;
%  # Langevin(beta)
Lng_b=(coth(b)-1./b);
Vd_f=p1(2)*p1(1)*(1-Lng_b.^2-2*Lng_b./b); % THIS IS CORRECT
sigma_square_hat=(1/(length(Force_new)-2))*sum(sep_residual.^2);
sigma_square=(1/(length(Force_new)-2))*sum(sep_residual.^2./Vd_f); % This is CORRECT 
sep_resid_weight=(1./Vd_f); % redefined weights (NEW corrected on May 2005)


function xp=find_points(force)
hold on
disp('zoom in')
zoom on
waitforbuttonpress;  
% Initially, the list of points is empty.
xy = [];
n = 0;
disp('select the region that contain junk adhesion force')
disp('select the approx region for the fit')
but = 1;
while but == 1
    [xpi,ypi,but] = ginput(1);
    xpi=round(xpi);
    ypi=round(ypi);
    plot(xpi,force(xpi),'ro')
    n= n+1;
    xp(:,n) = [xpi];
    yp(:,n) = [ypi];
end;
hold off    

function  [x,high_force,low_force]=find_points_moving_window_02(force,point_high_force,point_low_force,xp)
x=zeros(2,1);
end_point=length(force)
force=force(xp:end_point);
force=-force;
%Here we will find the high force point 
[qq,ww]=find(force>=point_high_force);
find_force=force(qq);
x(2)=qq(1);
high_force=find_force(1);
%%Here we will find the low force point using the moving window approach
window_size=10;
window_increment=1;
ii=x(2);
try
    while ii>0
        force_mean=mean(force((ii-window_size):ii));
        if force_mean>point_low_force
            ii=ii-window_increment;
            continue;
        else
            %%%%%%%%%%%%%%%%%%%%%Some additional predictor window 
            low_force=force(ii-window_increment);
            x(1)=ii-window_increment;
            break;
        end
    end 
    
catch 
    disp('blin')
end 

%%Fuction introduces the drag force correction
function [force,force_wdrag,Force_new,dForce,Sep_Fit] = find_drag_contrib(xp,SepNew,force);
time=1;
Numpoints=length(SepNew);
time_step=1/Numpoints;
TIME=time_step:time_step:time;
delta=TIME(2)-TIME(1);

b=[xp(2) xp(3)];
T_new=(TIME(b(1):b(2)))';

nn=5;
[pp,s]=polyfit(T_new,SepNew(b(1):b(2)),nn);
Sep_Fit=polyval(pp,T_new);
Sep_difference=Sep_Fit-SepNew(b(1):b(2));
stand_dev_sep=std(Sep_difference);
DSep_Fit=gradient(Sep_Fit,delta);  %%%%%%%%%%%FOR THE VELOcITY OF 1 MICRON/SEC
%DSep_Fit=gradient(Sep_Fit,2*delta); %%%%%%%%%%%FOR THE VELOcITY OF 0.5 MICRON/SEC
%dForce=0.0242*DSep_Fit; %%%%%%for the TM in hexodecane 
dForce=0.0069*DSep_Fit; %%%%%%for the TM in water 
%dForce=zeros(length((b(1):b(2))),1); % NO drag force 
Force_new=force(b(1):b(2))+dForce;

force_wdrag=force;
force=[force(1:b(1)-1);Force_new;force(b(2)+1:end)];
figure(8)
hold on
plot(SepNew(xp(2):xp(3)),Force_new,'r');
plot(SepNew(xp(2):xp(3)),force_wdrag(xp(2):xp(3)),'b');
hold off
%hold on
%plot(force_test,'r')
%plot(force)
%hold off

%fuction weights the data points based on the kernel density estimate
%Error function to the FJC model to the force extension data 

function err = fjcmin_resid_weight(p1,x,f,T,sep_resid_weight,sigma_square)
%weighted_force = f./sep_resid_weight;
KT = 4.1*((T+273.15)/298); % pNnm (4.114E-21 J at 298 K J=Nm=1E21 pNnm 298 - 273.15 = 24.85 C = 76.73 F)
z = zeros(length(f),1);
z = p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1)));
err = sum(sep_resid_weight.*(z-x).^2-sigma_square*log(sep_resid_weight)); % CORRECT
%err = norm(sep_resid_weight.*z-sep_resid_weight.*x); % OLD WEIGHTING 

function err = fjcmin_noweight(p1,x,f,T)
KT = 4.1*((T+273.15)/298); % pNnm (4.114E-21 J at 298 K J=Nm=1E21 pNnm 298 - 273.15 = 24.85 C = 76.73 F)
z = zeros(length(f),1);
z = p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1)));
err = norm(z-x);

%Fit the force extension data based on the pars of the FJC model
function z = fjc(f,p1,T)
KT = 4.1*((T+273.15)/298); % pNnm (4.114E-21 J at 298 K J=Nm=1E21 pNnm 298 - 273.15 = 24.85 C = 76.73 F)
z = zeros(length(f),1);
z =p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1)));

%Calculates the separation residuals 
function sep_residual=calc_sep_residual_FJC_new_appr_noplot(p1_kernel,Force_new,Sep_selected);
T=22;
KT = 4.1*((T+273.15)/298);
num_point=1;
sep_residual=zeros(1,1);
Index_matrix=zeros(1,1);
p1=p1_kernel;
while num_point < length(Force_new);
    f=-Force_new(num_point);
    z = p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1)));
    %z = p1(3)+p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1))); % Sep offset 
    %z=p1(2)*(coth((f+p1(3)).*p1(1)/KT)-KT./((f+p1(3))*p1(1))); %force offset
    sep_residual_element=Sep_selected(num_point)-z; 
    sep_residual=[sep_residual;sep_residual_element]; 
    num_point=num_point+1;
end

%function calculates the force residuals
function force_residual=calc_force_residual_no_smooth(Force_new,ZZ1_kernel_mod,ff1_mod,Sep_selected);
num_point=1;
force_residual=zeros(1,1);
Index_matrix=zeros(1,1);
Sep_selected_matrix=ones(length(ZZ1_kernel_mod),1);
while num_point < length(Force_new);
    [Y,ind] = min(abs((ZZ1_kernel_mod-Sep_selected_matrix*Sep_selected(num_point)))); 
    force_residual_element=ff1_mod(ind)+Force_new(num_point); 
    force_residual=[force_residual;force_residual_element]; 
    Index_matrix=[Index_matrix;ind]; 
    num_point=num_point+1;
end

function Short_fname_num=find_file_name(k,filename);
if k==1;
    Short_fname=filename(end-6:end-4)
    Short_fname_num = str2num(Short_fname);
elseif k==2
    Short_fname=filename(end-6:end-4);
    Short_fname_num = str2num(Short_fname);
    Short_fname_num_k2=str2num(strcat(Short_fname,'.2'));
    Short_fname_num = [Short_fname_num; Short_fname_num_k2];
elseif k==3
    Short_fname=filename(end-6:end-4);
    Short_fname_num = str2num(Short_fname);
    Short_fname_num_k3=str2num(strcat(Short_fname,'.3'));
    Short_fname_num = [Short_fname_num; Short_fname_num_k3];
end 

