% Multiple Curve Model  Fitting 
% 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

% Variables 
% dirname - the name of the directory where you store the force curves 
% kuhn_i - the initiatial values for the Kuhn length parameter of the FJC
% model (Value of kuhn_i does not change 
% kuhn - the Kuhn value that we increment during our stagewise minimization.
% dk - is the Kuhn length increment that we use in the minimization. The
% incremenet can be large or smaller than given depending on accuracy you
% want and computer power. Smaller increment will give better accuracy
% but takes it takes longer to simulate 
% T - is the value of the room temperature
% NP - number of points in your force curve 
% N_iter - number of iterations used in the stagewise minimization 
% iterN - is the current iteration number 
%  ii -  index of the forces in the SF_M matrix
% F - is the force vector from your force curve
% S - is the distance vector from your force curve
% FS- force in the selected region of the force curve i.e., the segment between the minimum and maximum force
% SS - separation in the selected region of the force curve

% FUNCTIONS
% find_drag - the program subtracts the drag force contribution

function Fit_FJC_fix_kuhn
clear all
close all
dirname='D:\LAPTOP\DISK_C\DIWORK\2003\15March_2003\test';

curdir = pwd;
files = dir(dirname);
numfiles = length(files) - 2;

load points150303_test_program
points=points150303;

T=25;
indL=7;
kuhn_i=0.35;
kuhn=kuhn_i;
dk=0.002; 
NP=4096;

ii=1;
nd=5;
N_iter=2;
iterN=1;

kuhn_bM=zeros(1,1);
SF_M=zeros(NP,1);

while ii<=numfiles
    filename = files(ii+2).name;
    eval(['cd ',dirname]); 
    load(filename);
    eval(['cd ',curdir]);
    endmat=length(filename)-4; 
    eval(['S=SepNew',filename(indL:endmat),';']);
    eval(['F=force',filename(indL:endmat),';']);
    
    % Here we put the force and separation from the your force curves into one big matrix
    S=sort(S);
    SF=[S,F];
    SF_M=[SF_M, SF];
    ii=ii+1;
end
SF_M(:,1)=[];

while  iterN <= N_iter; 
    ii=1;
    rr=1;
    %  ii is the index of the force curve
    while ii < size(SF_M,2)
        close all
        S=SF_M(:,ii);
        F=SF_M(:,ii+1);
        b=[points(rr,1) points(rr,2)]; 
        
        [F,fwdrag,FS,dF,S_Fit] = find_drag(b,S,F,NP);
        FS = F(b(1):b(2));
        SS = S(b(1):b(2));
        SS= sort(SS);
        L_i =200; 
        
        if ii==1;
            sigma_square_matr=zeros(1,1);
            err_M=zeros(1,1);
        end
        [L,kuhn_bM,sigma_square_matr,err_M]=fit_data_FJC_weight_Langevin(SS,FS,ii,kuhn,kuhn_bM,sigma_square_matr,err_M);
        
        L_kuhn_matr=[kuhn,L];      
        ff1_mod=linspace(0,-FS(end),length(SS))';
        ZZ1_mod = fjc_fix_kuhn(ff1_mod,L,kuhn);
        warning off MATLAB:divideByZero
        %%%%%%%%%%%%%%%%Calculate Separation residuals%%%%%%%%%%%%%%%%%%%
        sep_residual=calc_sep_residual_FJC_new_appr_noplot(p1_bootstr_matr,FS,SS);        
        if rr==1 & iterN==1
            Normalized_chi2_sep=zeros(1,1); 
            Normalized_chi2_sep_matrix=zeros(1,1);
            sum_sep_residual_matrix=zeros(1,1);
            Normalized_chi2_sep= calculate_chi2(sep_residual,L);
            %Normalized_chi2_sep_matrix(iterN,rr)=Normalized_chi2_sep;
            Normalized_chi2_sep_matrix(rr)=Normalized_chi2_sep;
            Apar_kernel=zeros(1,2,1);  
            Apar_kernel(rr,:,iterN)=p1_bootstr_matr;%[Apar_kernel; p1_bootstr_matr];
        else   
            [Normalized_chi2_sep,sum_sep_residual]= calculate_chi2(sep_residual,L);
            Normalized_chi2_sep_matrix(rr)=Normalized_chi2_sep;
            sum_sep_residual_matrix(rr)=sum_sep_residual;
            Apar_kernel(rr,:,iterN)=p1_bootstr_matr;
        end 
        
        ii=ii+2;    
        rr=rr+1;
    end  
    if iterN==1;
        err_MA=zeros(1,1);
    end
    %Here is the error function I am MINIMIZING 
    err_MA=[err_MA;sum(sigma_square_matr)];
    
    %Sum of squared residuals are saved at each iteration (coolumn) and at
    %each boostrap (row)
    
    iterN=iterN+1;
    kuhn=kuhn+dk;
    kuhn_matrix(:,iterN)= kuhn;
    disp('OK')
end  
err_MA(1,:)=[];
%[min_chi2_sep,ind_min_chi2_sep]=min(sigma_square_matr_all_curves);%min(sep_chi2_matrix);
[min_err_val,ind_err_val]=min(err_MA);

%BELOW SEVERAL LINES CACULATE THE KUHN LENGTH THAT CORRSSPOINDS TO THE
%MINIMUM SUM OF SQUARED RESIDUALS
[min_err_val,ind_err_val]=min(err_MA);
kuhn_err_value=kuhn_i+dk*ind_err_val;
%Finished finding the orinial data parameters%%%%%%%%%

%BELOW ARE THE BOOTSTRAPS CALCULATIONS%%%
%%BOOTSTRPING AND ESTIMATING CONFIDENCE INTERVALS%%%

perm_vec=[1:2:6]';

%Changing nd will affect the number of iterations and influence the
%execution time 
%%Bootsam containes the indexes randomly drawn nd times from the perm_vec
[bootstat,bootsam] = bootstrp(nd,'mean',perm_vec); % now bootstrap_test is the fake function because bootstrp rtequires soem function handle  
rand_perms=perm_vec(bootsam); %%%the nd columns of the random draws obtained from perm_vec draws
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Changing N_iter will affect the number of iterations and influence the execution time 
%%%Will determine the range over which the kuhn length 
sep_chi2_matrix=zeros(N_iter,nd);
p1_kernel_matrix=zeros(N_iter,3);
sep_chi2_matrix_bootstrap=zeros(N_iter,N_iter,nd);

kuhn_bM=zeros(1,1);
L_kuhn=zeros(numfiles,N_iter,nd);
kuhn_matrix=zeros(1,1,1);
hh=1;
while hh<=nd
    %%%while ss will generate the permuated matrix of the the separations and forces
    ss=1;
    sep_force_matrix_perm=zeros(length(F),1);
    points_perm=zeros(1,2);
    while ss<=numfiles %nd
        sep_force_matrix_perm=[sep_force_matrix_perm, SF_M(:,rand_perms(ss,hh)),SF_M(:,rand_perms(ss,hh)+1)];
        points_perm=[points_perm; points(0.5*(rand_perms(ss,hh)+1),1) points(0.5*(rand_perms(ss,hh)+1),2)];
        ss=ss+1;
    end 
    sep_force_matrix_perm(:,1)=[];
    points_perm(1,:)=[];
    iterN=1;
    kuhn=kuhn_i;
    while  iterN <= N_iter; 
        ii=1;
        rr=1;
        %  ii is the index of the force curve
        while ii <= size(sep_force_matrix_perm,2)
            close all
            S=sep_force_matrix_perm(:,ii);
            F=sep_force_matrix_perm(:,ii+1);
            b=[points_perm(rr,1) points_perm(rr,2)]; 
            Force_selected = F(b(1):b(2));
            SS = S(b(1):b(2));
            FS=Force_selected;
            L_i =200; 
            T=22;        
            if ii==1;
                sigma_square_matr=zeros(1,1);
            end
            [L,kuhn_bM,sigma_square_matr,err_M]=fit_data_FJC_weight_Langevin(SS,FS,ii,kuhn,kuhn_bM,sigma_square_matr,err_M);
            
            p1_bootstr_matr=[kuhn,L];      
            ff1_mod=linspace(0,-FS(end),length(SS))';
            ZZ1_mod = fjc_fix_kuhn(ff1_mod,L,kuhn);
            warning off MATLAB:divideByZero
            %%%%%%%%%%%%%%%%Calculate Separation residuals%%%%%%%%%%%%%%%%%%%
            sep_residual=calc_sep_residual_FJC_new_appr_noplot(p1_bootstr_matr,FS,SS);        
            if rr==1 
                Normalized_chi2_sep=zeros(1,1); 
                Normalized_chi2_sep_matrix=zeros(1,1);
                [Normalized_chi2_sep,sum_sep_residual]= calculate_chi2(sep_residual,L);
                Normalized_chi2_sep_matrix(rr)=Normalized_chi2_sep;
                sum_sep_residual_matrix(rr)=sum_sep_residual;
                Apar_kernel=zeros(1,2);  
                Apar_kernel=[Apar_kernel; p1_bootstr_matr];
            else   
                [Normalized_chi2_sep,sum_sep_residual]= calculate_chi2(sep_residual,L);
                
                Normalized_chi2_sep_matrix(rr)=Normalized_chi2_sep;
                sum_sep_residual_matrix(rr)=sum_sep_residual;
                Apar_kernel=[Apar_kernel; p1_bootstr_matr];  
            end 
            L_kuhn(rr,iterN,hh)= Apar_kernel(rr+1,2);  
            ii=ii+2;    
            rr=rr+1;
        end  
        
        err_value_matr_all_curves_boot(iterN,hh)=sum(sigma_square_matr);%[err_value_matr_all_curves; sum(sigma_square_matr)];
        
        iterN=iterN+1;
        kuhn=kuhn+dk;
        kuhn_matrix(:,iterN,hh)= kuhn;
        disp('OK')
    end  
    hh=hh+1; %Go to the next bootstrap 
    iter=num2str(hh); disp(iter)
end 
CL=clock; digits(4); Hour=vpa(CL(4)); Minute=vpa(CL(5));
%%%%%%%%%%%%%%%%%%old line %%%%%%%%%%%%%%5
[min_chi2_sep,ind_min_chi2_sep]=min(err_value_matr_all_curves_boot);
%We find the value of the kuhn lenght that minimize the squared residuals at each iteration%%%
for ff=1:length(ind_min_chi2_sep)
    kuhn_min(ff)=kuhn_i+dk*(ind_min_chi2_sep(ff)-1);
end
%%Find the mean value for the kuhn length. This is a number 
mean_kuhn=sum(kuhn_min)/length(kuhn_min);
%Define the term needed to stimate the the bootstrap uncertinity
for pp=1:nd
    BB_kuhn(pp)=(kuhn_min(1,pp)-mean_kuhn)^2;
end 
%points=points';
uncert_kuhn_bootstrap=((1/(nd-1))*sum(BB_kuhn))^0.5;

%%Calculating uncertainity in the contour length%%
%%Determine the minimum over the interations 
%The euncertinity will be calculated for each individual curve. I.e the
%final result is a vector in length equal to number of curves. L_kuhn is
%the 3 D matrix. We sum over iterations i.e over the second dimension to
%get the 2D matrix. row is the curve number, column is the bootstrap
for tt=1:size(L_kuhn,3)
    contour_min(:,tt)=L_kuhn(:,ind_min_chi2_sep(tt),tt);
end 
%%%Sum over the curve number%%%
for gg=1:size(L_kuhn,1)
    sum_contour_min(gg)=sum(contour_min(gg,:));
end
%mean contour length for each individual curve%%
mean_contour=sum_contour_min/size(L_kuhn,2);
dd=1;
%calculate the term nessesary to estimate the uncertinity matrix
%%HERE I HAVE SOME PROBLEM....

for dd=1:size(L_kuhn,1)
    BB_contour(:,dd)=(contour_min(:,dd)-mean_contour').^2;
end 
dd=1;    
% uncert_contour_bootstrap is a vector in length equal to number of curves 
for dd=1:size(L_kuhn,1)
    uncert_contour_bootstrap(dd)=((1/(nd-1))*sum(BB_contour(dd,:)))^0.5;
end

%%%%%%%%%%%%Calculate the total number of points on all curves%%%%%%%%%%%%%%%%%%%%%
Np=zeros(length(points),1);
for i=1:length(points)
    Np(i)=points(i,2)-points(i,1);
end 
sum_Num_points=sum(Np);



convertdir = 'par_fix_Kuhn_weight_Langevin'; 
mkdir(dirname,convertdir)
savename=[dirname,'\',convertdir,'\','pFJC_fix_Kuhn_Langevin',filename(7:(end-8))]; 
Kuhn_matrix=['L_kuhn',filename(9:(endmat-4))];
sep_chi2=['sep_chi2_matrix',filename(9:(endmat-4))];
Kuhn_Min=['kuhn_min',filename(9:(endmat-4))];
Kuhn_Init=['kuhn_i',filename(9:(endmat-4))];   
Uncert_Kuhn=['uncert_kuhn',filename(9:(endmat-4))];
Mean_Kuhn=['mean_kuhn',filename(9:(endmat-4))];
Uncert_Contour=['uncert_contour',filename(9:(endmat-4))];  
Points=['points',filename(9:(endmat-4))];
SumNumPoints=['SumNumPoints',filename(9:(endmat-4))];
Delta_Kuhn=['dk',filename(9:(endmat-4))];
% HERE THE FORCE AND SEPARATION AT RAPTURE POINT IS RECORDED AT THE MAT FILE

eval([sep_chi2,' = sep_chi2_matrix;']);
eval([Kuhn_matrix,' = L_kuhn;']);
eval([Kuhn_Min,' = kuhn_min;']);
eval([Kuhn_Init,' = kuhn_i;']);
eval([Mean_Kuhn,' = mean_kuhn;']);
eval([Uncert_Kuhn,' = uncert_kuhn_bootstrap;']);
eval([Uncert_Contour,' =  uncert_contour_bootstrap;']);
eval([Points,' = points;']);
eval([SumNumPoints,' = sum_Num_points;']);
eval([Delta_Kuhn,' = dk;']);
eval(['save(savename,''',sep_chi2,''',''',Kuhn_Min,''',''',Kuhn_Init,''',''',Uncert_Kuhn,''',''',Mean_Kuhn,''',''',Uncert_Contour,''',''',Kuhn_matrix,''',''',Points,''',''',Delta_Kuhn,''')']);
%%%%%%%%%%%%%%%
%%%%%%%%%%Plotiing the distribution of Kuhn segment lengths %%%%%%%%%%%
kuhn_sort=sort(kuhn_min);
interval_start=min(kuhn_min)-0.1;
interval_end=max(kuhn_min)+0.1;

interval=interval_start:dk:interval_end;
kuhn_group=histc(kuhn_sort,interval);
figure(1)
bar(interval,kuhn_group,'histc'),
title('Kuhn Length','FontSize',12,'Color',[0.5 0.1 0.1])
xlabel('Kuhn Length [nm]','FontSize',20,'Color',[1 0.1 0.3])
ylabel('Ratio of occurrences','FontSize',20,'Color',[1 0.1 0.3])
set(gca,'FontSize',24)


%%%Lets built the gaussian with the corresponding mean and the confidence interval%%
sigma=uncert_kuhn_bootstrap;
mu=mean_kuhn;
W=(interval-mu).^2./(2.*sigma^2);
Gaus=1/(sigma*sqrt(2*pi)).*exp(-W);
conf_interval=[mu-sigma*1.96;mu+sigma*1.96]
bound_low=ones(1,10)*conf_interval(1);
bound_high=ones(1,10)*conf_interval(2);
y_low=0:2:18;
y_high=0:2:18;
figure(2);
hold on
bar(interval,kuhn_group,'histc','b')
plot(interval,Gaus,'r')
plot(bound_low,y_low,'g','LineWidth',2)
plot(bound_high,y_high,'m','LineWidth',2)
hold off
grid on
title('Kuhn Length','FontSize',12,'Color',[0.5 0.1 0.1])
xlabel('Kuhn Length [nm]','FontSize',20,'Color',[1 0.1 0.3])
ylabel('Ratio of occurrences','FontSize',20,'Color',[1 0.1 0.3])
set(gca,'FontSize',24)
axis([0.3 0.7 0 100])


function [L,kuhn_bM,sigma_square_matr,err_M]=fit_data_FJC_weight_Langevin(SS,FS,ii,kuhn,kuhn_bM,sigma_square_matr,err_M);
L_i=200;
ss=1;
delta_kuhn=1;
T=22; 
while ss<=50 
    if  abs(delta_kuhn)<=0.01
        break;
    else 
        if ss==1;
            [L,fval,exitflag,output] = fminsearch(@fjcmin_fix_kuhn_no_weight,L_i,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),SS,-FS,kuhn);
        else
            [L,fval,exitflag,output] = fminsearch(@fjcmin_fix_kuhn_weight,L_i,optimset('Tolx',1e-10,'MaxFunEvals',3000,'MaxIter',7000),SS,-FS,sep_resid_weight,kuhn,sep_resid_weight,sigma_square);
            delta_kuhn=kuhn_bM(ss)-kuhn_bM(ss-1);    
        end
        kuhn_bM=[kuhn_bM;L];
        ff1_mod=linspace(0,-FS(end),length(SS))';
        ZZ1_mod = fjc_fix_kuhn(ff1_mod,L,kuhn);
        warning off MATLAB:divideByZero         
        %%%%%%%%%%%%%%%%Calculate Separation residuals%%%%%%%%%%%%%%%%%%%                
        [sep_residual,sep_resid_weight,sigma_square,sigma_square_hat,Vd_f,err_value]=calc_sep_residual_weight(L,FS,SS,kuhn);
        %Lets add in some control parameters 
    end 
    ss=ss+1;
end
err_M=[err_M;err_value];
sigma_square_matr=[sigma_square_matr;sigma_square];


function [sep_residual,sep_resid_weight,sigma_square_mod,sigma_square_hat,Vd_f,err_value]=calc_sep_residual_weight(L,FS,SS,kuhn);
%function [sep_residual,sep_resid_weight,sigma_square,sigma_square_hat,Vd_f,Rss]=calc_sep_residual_weight(L,FS,SS,kuhn);
T=22;
KT = 4.1*((T+273.15)/298);
num_point=1;
sep_residual=zeros(1,1);
Index_matrix=zeros(1,1);
while num_point < length(FS);
    f=-FS(num_point);
    z =L*(coth(f.*kuhn/KT)-KT./(f*kuhn));
    sep_residual_element=SS(num_point)-z; 
    sep_residual=[sep_residual;sep_residual_element]; 
    num_point=num_point+1;
end
%%Weight sep residuals
%%%%%%%%%%%%%%%%%%%%OLD WEIGHTING %%%%%%%%%%%%%%%%%%%%%%%%%%%
%sigma_square=sum((1/(length(FS)-2))*sep_residual.^2);
%L_B=(coth(-FS.*kuhn./KT)-KT./(-FS.*kuhn));
%aa=KT./(-FS.*kuhn);
%coef=2.*aa.*L_B;
%bb=(L.*kuhn.*(1-L_B.^2-coef));
%sep_resid_weight=(1./(sigma_square*bb));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

b=kuhn*(-FS)./KT;
%  # Langevin(beta)
Lng_b=(coth(b)-1./b);
Vd_f=L*kuhn*(1-Lng_b.^2-2*Lng_b./b);
sigma_square_hat=(1/(length(FS)-2))*sum(sep_residual.^2);
sigma_square=(1/(length(FS)-2))*sum(sep_residual.^2./Vd_f);
sep_resid_weight=(1./Vd_f);
sigma_square_mod=sum((sep_residual.^2.*sep_resid_weight)-sigma_square_hat*log(sep_resid_weight));
%%%%%%%%%%%%Rss is the same as SSR -Sum of Swuared weighted residuals%%%%%%%%%%%%%%%%
%HERE THE CORRECT TVERSION OF THE ERROR FUNCTION IS USED. THIS CORRECTION
%IS INTRODUCED IN MAY 2005
err_value=(1/sigma_square)*sum(sep_resid_weight.*sep_residual.^2-sigma_square*(log(sep_resid_weight)-log(sigma_square)));



function err = fjcmin_fix_kuhn_weight(p1,x,f,sep_resid_weight,kuhn,sep_residual,sigma_square)
%err = fjcmin_fix_kuhn_weight(p1,x,f,sep_resid_weight,kuhn,sigma_square);
T=22;
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(1)*(coth(f.*kuhn/KT)-KT./(f*kuhn));

%sigma_square=(1/(length(f)-2))*sum(sep_residual.^2.*sep_resid_weight);
%%%OLD ERROR FUNCTION%%%
%err = norm(sep_resid_weight.*z-sep_resid_weight.*x);
%err = sum(sep_resid_weight.*(z-x).^2-sigma_square*log(sep_resid_weight));
err = (1/sigma_square)*sum(sep_resid_weight.*(z-x).^2-sigma_square*(log(sep_resid_weight)-log(sigma_square)));

function err = fjcmin_fix_kuhn_no_weight(p1,x,f,kuhn)
T=22;
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(1)*(coth(f.*kuhn/KT)-KT./(f*kuhn));
err = norm(z-x);

function z = fjc_fix_kuhn(f,p1,kuhn);
T=22;
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(1)*(coth(f.*kuhn/KT)-KT./(f*kuhn));

function sep_residual=calc_sep_residual_FJC_new_appr_noplot(p1_bootstr_matr,FS,SS);
T=22;
KT = 4.1*((T+273.15)/298);
num_point=1;
sep_residual=zeros(1,1);
Index_matrix=zeros(1,1);
p1=p1_bootstr_matr;
while num_point < length(FS);
    f=-FS(num_point);
    z = p1(2)*(coth(f.*p1(1)/KT)-KT./(f*p1(1)));
    sep_residual_element=SS(num_point)-z; 
    sep_residual=[sep_residual;sep_residual_element]; 
    num_point=num_point+1;
end 


function [Normalized_chi2_sep,sum_sep_residual]= calculate_chi2(sep_residual,L);
T=22;
KT = 4.1*((T+273.15)/298);
%p1=L;
sep_residual_square=sep_residual'*sep_residual;
sum_sep_residual=sum(abs(sep_residual));
Normalized_chi2_sep=sep_residual_square/length(sep_residual);

function   [F,fwdrag,FS,dF,S_Fit] = find_drag(x,S,F,NP);
% find_drag - the program subtracts the drag force contribution
% time - the number of seconds it takes to go 1 micron distance. This value will depend on your pulling speed 
% fm- the mean value of force in the window
% x - is the vector that contains the index of the points corresponding to minimum and maximum force
% dF - is the drag force contribution. Note that multiplier 0.0069 is for a
% specific contilever. You will need to find this coefficient for the
% cantilever you will use.

% The cantilever retraction velocity (speed) relative to the surface will
% change depending on force exerted on the ecantilever. We thus defined a velocity at each point of
%the force separation data by differentiating separation with respect
%to time. One can estimate this derivative based on the actual
%separation data. However, sensor noise influences the estimate from the actual data.
%Instead we fit the data selected 5{th} order polynomial. We found that 
%this polynomial described the data reasonably well. Then we differentiated
%that polynomial to find relative velosity.
time=1;
time_step=1/NP;
t_vec=time_step:time_step:time;
delta=t_vec(2)-t_vec(1);
Tn=(t_vec(x(1):x(2)))';

nn=5;
warning off MATLAB:polyfit:RepeatedPointsOrRescale
[pp,s]=polyfit(Tn,S(x(1):x(2)),nn);
S_Fit=polyval(pp,Tn);
Sd=S_Fit-S(x(1):x(2));
DSep_Fit=gradient(S_Fit,delta);  %%%%%%%%%%%FOR THE VELOcITY OF 1 MICRON/SEC

dF=0.0069*DSep_Fit; %%%%%%for the TM in water 
FS=F(x(1):x(2))+dF/2;
fwdrag=F;
F=[F(1:x(1)-1);FS;F(x(2)+1:end)];