function  [Result_String] = forcecalc6(alpha,P,C0,g,plot_flag,fig_hand,normal_sign,zmin,zmax,whichparam)

% get right figure and clear it

figure(fig_hand);
clf;

%g = data_maker2(g,.4,.1);

% plot_flag=3;
% coarsegrain data

    loopsize = max(size(g));
    z  =g(1:loopsize,2);
    r  =g(1:loopsize,1);
    
    %kc = 15;
    kc=1;
    kT = 4.143e-3;
    kT = 1;

            
    LE = 0.1;

    % N is the size of the data array
    N = max(size(z));
    N1 = N-1;
    N2 = N-2;
    N3 = N-3;
    N4 = N-4;
    
    
    
    rp = (r(2:N)+r(1:N1))/2;
    
    
    %ds is differential arc length    
    ds  = sqrt((z(2:N)-z(1:N1)).^2+(r(2:N)-r(1:N1)).^2);
    
    %s is the arclength
    s=1:N;
    s(1)=0;
    for i=2:1:N;
        s(i) = s(i-1)+ds(i-1);
    end;

    area = sum( ds.*rp*2*pi);
    
    %tz and tr are the componentsof the unit tangent vector
    % first we initialize the tangents
    
    tr  = (r(2:N)-r(1:N1))./(ds(1:N1));
    tz  = (z(2:N)-z(1:N1))./(ds(1:N1));

    trp = 0*r;
    tzp = 0*r;
    dsp = 0*r;
    
    trp(2:N1) = (tr(2:N1)+tr(1:N2))/2;   
    tzp(2:N1) = (tz(2:N1)+tz(1:N2))/2; 
    dsp(2:N1) = (ds(2:N1)+ds(1:N2))/2; 
        
    
    trp(1) = tr(1);
    trp(N) = tr(N1);
    
    tzp(1) = tz(1);
    tzp(N) = tz(N1);
    
    dsp(1) = ds(1);
    dsp(N) = ds(N1);

    
        l = dsp;

    nr = tz*normal_sign;
    nz =  -tr*normal_sign;
    nrp = tzp*normal_sign;
    nzp =  -trp*normal_sign;
    
    
    Cr  = 0*r;
    Cz  = 0*r;
    SNN = 0*r;
    STT = 0*r;
    

    Cr(2:N1) = (tr(1:N2)-tr(2:N1))./dsp(2:N1);
    Cz(2:N1) = (tz(1:N2)-tz(2:N1))./dsp(2:N1);

    SNN(2:N1)= -(Cr(2:N1).*nrp(2:N1)+Cz(2:N1).*nzp(2:N1));

    STT(2:N1)= -nrp(2:N1)./r(2:N1);

     SNN(1)=-(z(2)-z(1))./r(2)^2*tz(2);
     STT(1)=SNN(1);

%    SNN(N)=-(z(N)-z(N1))./r(N1)^2*tz(N1);
%    STT(N)=SNN(N);
   
    S = SNN+STT;
    K = SNN.*STT;

    ac = sum(S.*r.*dsp*2*pi/area);
    
    C2r = STT.*nrp;
    C2z = STT.*nzp;

    % DS = dS/ds
    % init
    DS = 0*tr;

    DS(1:N1) = (S(2:N)-S(1:N1))./ds;

    rDS = rp.*DS;
    
    D2S = 0*r;
    D2S(2:N1) = (rDS(2:N1)-rDS(1:N2))./dsp(2:N1);
    D2S(2:N1) = D2S(2:N1)./r(2:N1);


    D2S(1) = -(S(2)-S(1))*2/r(2)^2;
    %D2S(N) = (S(N)-S(N1))*2/r(N1)^2;

    fk = kc*2*K.*S;
    fs = -1/2*S.^3*kc;
    fd = -D2S*kc;
    
    %ffp = 1+0*S;
    ffp = -normal_sign+0*S;
    ffa = S;
    ff0 = fk+fs+fd;
    ffc = -2*K.*kc;
    
    fp = P*ffp;
    fa = alpha*ffa;
    fc = C0*ffc;
    
    f = fp+fa+fk+fs+fd+fc;
    
    loopsize = max(size(ffa));
    z_control = ones(loopsize,1);
    
    %zmin =7;
    %zmax = 20;
    
    for i=1:loopsize;
        if (z(i)<zmin) || (z(i)>zmax)
            z_control(i) = 0;
        end
    end;
    
    
    
    MMM = [[sum(dsp.*r.*ffa.*ffa.*z_control),sum(dsp.*r.*ffa.*ffp.*z_control)];[sum(dsp.*r.*ffa.*ffp.*z_control),sum(dsp.*r.*ffp.*ffp.*z_control)]];
    MM0  = [sum(dsp.*r.*ff0.*ffa.*z_control);sum(dsp.*r.*ff0.*ffp.*z_control)];
    if( det(MMM) ~= 0 )
        guess = -inv(MMM)*MM0;

    else
        guess = [0;0];
    end;
    
    %-----------------------------------------------%
    % Now do the same computation with spont curve  %
    % as well...                                    %
    %-----------------------------------------------%
    
    NNN = [[sum(dsp.*r.*ffa.*ffa.*z_control),sum(dsp.*r.*ffa.*ffp.*z_control),sum(dsp.*r.*ffa.*ffc.*z_control)];[sum(dsp.*r.*ffa.*ffp.*z_control),sum(dsp.*r.*ffp.*ffp.*z_control),sum(dsp.*r.*ffp.*ffc.*z_control)];[sum(dsp.*r.*ffa.*ffc.*z_control),sum(dsp.*r.*ffc.*ffp.*z_control),sum(dsp.*r.*ffc.*ffc.*z_control)]];
    NN0  = [sum(dsp.*r.*ff0.*ffa.*z_control);sum(dsp.*r.*ff0.*ffp.*z_control);sum(dsp.*r.*ff0.*ffc.*z_control)];
    if( det(NNN) ~= 0 );
        guessN = -inv(NNN)*NN0;
    else
        guessN = [0;0;0];
    end; 
    
    
    
    f_error_user = sum((f.*z_control).^2);

    fmin = guess(2)*ffp+guess(1)*ffa+fk+fs+fd+fc;
    f_error_min = sum((fmin.*z_control).^2);
    
    
    
    
    %    feT = sqrt(1/(2*pi)*(kc*pi^6./(6*l.^6)+alpha*pi^4./(4*l.^4)));
 %   feNT = sqrt(kc^2*pi^8*LE^2./(5*l.^8));

    
    
    
            circ = 2*pi*r;

%        Fint = -pi*P*r.^2+2*pi*r.*((kc*SNN.*(S-C0)-(kc*0.5*(S-C0).^2-.5*C0^2+alpha)).*tzp+[nz.*DS*kc;0]);


FFp = -pi*r.^2*normal_sign;
FFa = -2*pi*r.*tzp*normal_sign;
FFc = kc*2*pi*r.*(-SNN+S).*tzp;

Fd  = kc*2*pi*r.*[nz.*DS*kc;0]*normal_sign;
Fs  = kc*2*pi*r.*((SNN.*(S)-(0.5*(S).^2)).*tzp)*normal_sign;
Fds = Fs+Fd;
FF0 = Fds;


%---------------------------------------------------------%
%Manual Force

Fp = P*FFp;
Fa = alpha*FFa;
Fc = C0*FFc;

F = Fp+Fa+Fd+Fs+Fc;


Fcomp = -kT*sum(F.*z_control)/sum(z_control);
Ftemp  = (-kT*F-Fcomp).*z_control;
dFcomp     = sqrt(sum(Ftemp.^2)/sum(z_control));

disp 'Manual:'
disp(sprintf('F = %g +/- %g   a = %g   p = %g   c_0 = %g\n', Fcomp, dFcomp,alpha,P,C0))


%---------------------------------------------------------%
%Min2

Fp = guess(2)*FFp;
Fa = guess(1)*FFa;
Fc = 0*FFc;

F = Fp+Fa+Fd+Fs+Fc;
 
FcompM = kT*sum(F.*z_control)/sum(z_control);
Ftemp  = (kT*F-FcompM).*z_control;
dFcompM     = sqrt(sum(Ftemp.^2)/sum(z_control));

disp 'min2 (f)^2:'
disp(sprintf('F = %g +/- %g   a = %g   p = %g\n', FcompM, dFcompM,guess(1),guess(2)))
  





Fp = guessN(2)*FFp;
Fa = guessN(1)*FFa;
Fc = guessN(3)*FFc;

F = Fp+Fa+Fd+Fs+Fc;
 
FcompN = kT*sum(F.*z_control)/sum(z_control);
Ftemp  = (kT*F-FcompN).*z_control;
dFcompN     = sqrt(sum(Ftemp.^2)/sum(z_control));


                 disp 'min3 (f)^2'
         disp(sprintf('F = %g +/- %g   a = %g   p = %g   c_0 = %g\n',FcompN,dFcompN, guessN(1),guessN(2),guessN(3)))
         




       
        
Result_String = sprintf('%g +/- %g\n\n%g, %g\n%g +/- %g\n\n%g, %g, %g\n%g +/- %g\n\n%g, %g',Fcomp,dFcomp,guess(1),guess(2),FcompM,dFcompM,guessN(1),guessN(2),guessN(3),FcompN,dFcompN,area,ac);

Fint = Fp+Fa+Fd+Fs+Fc;


%-----------------------------------------%
%                                         %
% New fitting strategy for flatness of    %
% force..                                 %
%                                         %
%-----------------------------------------%

    FMMM = [[sum(dsp.*FFa.*FFa.*z_control),sum(dsp.*FFa.*FFp.*z_control),sum(-dsp.*FFa.*z_control)];[sum(dsp.*FFa.*FFp.*z_control),sum(dsp.*FFp.*FFp.*z_control),sum(-dsp.*FFp.*z_control)];[sum(-dsp.*FFa.*z_control),sum(-dsp.*FFp.*z_control),sum(dsp.*z_control)]];
    FMM0  = [sum(dsp.*FF0.*FFa.*z_control);sum(dsp.*FF0.*FFp.*z_control);sum(-dsp.*FF0.*z_control)];
    if( det(FMMM) ~= 0 )
        Fguess = -inv(FMMM)*FMM0;
    else
        Fguess = [0;0;0];
    end;
    
    Fp = Fguess(2)*FFp;
    Fa = Fguess(1)*FFa;
    Fc = 0*FFc;

F = Fp+Fa+Fd+Fs+Fc;
 
FcompF = kT*sum(F.*z_control)/sum(z_control);
Ftemp  = (kT*F-FcompN).*z_control;
dFcompF     = sqrt(sum(Ftemp.^2)/sum(z_control));


    disp 'min2 (dF)^2:'
    disp(sprintf('F = %g +/- %g   a = %g   p = %g\n', FcompF,dFcompF,Fguess(1),Fguess(2)))
    
        
    %-----------------------------------------------%
    % Now do the same computation with spont curve  %
    % as well...                                    %
    %-----------------------------------------------%
    
    FNNN = [[sum(dsp.*FFa.*FFa.*z_control),sum(dsp.*FFa.*FFp.*z_control),sum(dsp.*FFa.*FFc.*z_control),sum(-dsp.*FFa.*z_control)];[sum(dsp.*FFa.*FFp.*z_control),sum(dsp.*FFp.*FFp.*z_control),sum(dsp.*FFp.*FFc.*z_control),sum(-dsp.*FFp.*z_control)];[sum(dsp.*FFa.*FFc.*z_control),sum(dsp.*FFc.*FFp.*z_control),sum(dsp.*FFc.*FFc.*z_control),sum(-dsp.*FFc.*z_control)];[sum(-dsp.*FFa.*z_control),sum(-dsp.*FFp.*z_control),sum(-dsp.*FFc.*z_control),sum(dsp.*z_control)]];
    FNN0  = [sum(dsp.*FF0.*FFa.*z_control);sum(dsp.*FF0.*FFp.*z_control);sum(dsp.*FF0.*FFc.*z_control);sum(-dsp.*FF0.*z_control)];
    if( det(NNN) ~= 0 );
        FguessN = -inv(FNNN)*FNN0;
    else
        FguessN = [0;0;0];
    end; 

    
    Fp = FguessN(2)*FFp;
    Fa = FguessN(1)*FFa;
    Fc = FguessN(3)*FFc;

F = Fp+Fa+Fd+Fs+Fc;
 
FcompFN = kT*sum(F.*z_control)/sum(z_control);
Ftemp  = (kT*F-FcompFN).*z_control;
dFcompFN     = sqrt(sum(Ftemp.^2)/sum(z_control));

    disp 'min3 (dF)^2:'
    disp(sprintf('F = %g +/- %g   a = %g   p = %g   c_0 = %g\n',FcompFN,dFcompFN, FguessN(1),FguessN(2),FguessN(3)))









        
        
%-----------------------------------------%
% whichparam

switch whichparam;
    case 1;
        'manual';
    case 2;
        'min(f^2) fit 2';
        alpha = guess(1);
        P = guess(2);
        C0 = 0;
    case 3;
        'min(f^2) fit 3';
        alpha = guessN(1);
        P = guessN(2);
        C0 = guessN(3);
    case 4;
        'min(dF^2) fit 2';
        alpha = Fguess(1);
        P = Fguess(2);
        C0 = 0;
    case 5;
        'min(dF^2) fit 3';
        alpha = FguessN(1);
        P = FguessN(2);
        C0 = FguessN(3);        
end;

Fp = P*FFp;
Fa = alpha*FFa;
Fc = C0*FFc;
Fint = Fp+Fa+Fd+Fs+Fc;
fp = P*ffp;
fa = alpha*ffa;
fc = C0*ffc;
f = fp+fa+fk+fs+fd+fc;



                
    % pN = 1.;
%         size(r_)
%         size(DS_)
%         size(tz_)
%         size(SNN_)

         Fint_exp = Fint*0;
     
         for i=2:1:N2
             Fint_exp(i) = Fint_exp(i-1)+f(i)*nz(i)*2*pi*r(i)*ds(i);
         end;

        
% Decide which plot to show
switch plot_flag
    case 4
       
        plot(z,r,'.k');
        hold on;
        plot(z,r,'k');
        plot(z,-r,'.k');
        plot(z,-r,'k');
        ll = mean(ds);
        
        for i=1:N
            plot([z(i),z(i)+ll*nzp(i)],[r(i),r(i)+ll*nrp(i)],'r');
            plot([z(i),z(i)+ll*tzp(i)],[r(i),r(i)+ll*trp(i)],'b');
            plot([z(i),z(i)+ll*nzp(i)],-[r(i),r(i)+ll*nrp(i)],'r');
            plot([z(i),z(i)+ll*tzp(i)],-[r(i),r(i)+ll*trp(i)],'b');
        end;
         axis equal;
         
%     plot(z(1:N1),nr);
%     hold on;
%     plot(z(1:N1),nz+1,'--');
%     plot(z(1:N1),tr,'g');
%     plot(z(1:N1),tz,'r--');
    %add_background_conf(g,zmin,zmax)
    case 3
        %plots curvatures
        plot(z(1:N1),S(1:N1),'k');
        hold on;
        plot(z(1:N1),K(1:N1),'g');
        plot(z(1:N1),SNN(1:N1),'c');
        plot(z(1:N1),STT(1:N1),'b');
        %plot(z(1:N2),D2S(1:N2),'r');
        %plot(z(1:N),D2S(1:N),'r--');
  
        legend('S','K','S_{NN}','S_{TT}');
        xlabel('x ({\mu}m)');
        ylabel('S_{ij}');
        
       % plot(z(1:N1),S(1:N1),'k.');
        hold on;
       %plot(z(1:N1),K(1:N1),'g.');
       %plot(z(1:N1),SNN(1:N1),'c.');
       %plot(z(1:N1),STT(1:N1),'b.');
        %plot(z(1:N2),D2S(1:N2),'r');
        %plot(z(1:N),D2S(1:N),'r--');
        grid;
        add_background_conf(g,zmin,zmax)
    case 2
      

    plot(z(1:N),f(1:N),'LineWidth',3);
    hold;
    plot(z(1:N),fd(1:N),'y');
    plot(z(1:N),fs(1:N)+fk(1:N),'c');
    
    %plot(s,fs,'g');
    %plot(s,fk,'y');
    plot(z(1:N),fa(1:N)+fp(1:N),'r');
    plot(z(1:N),fp(1:N),'m');
    plot(z(1:N),fa(1:N),'g');

    
    %plot(z(1:N),fc(1:N),'y');
    
    plot(z(1:N),fc(1:N),'k');
    
    %plot(z(1:N),f(1:N),'b.');
    
    
    legend('f','f_{d}','f_{elast}','f_{\alpha+P}','f_P','f_\alpha','f_C');
    
    grid;
    
    xlabel('x ({\mu}m)');
    ylabel('f (kT/{\mu}m^3)');
    
    add_background_conf(g,zmin,zmax)

    case 1  

        pN = 1;
        
        plot(z(2:N1),kT*Fint(3:N),'LineWidth',3);

                hold on;

        plot(z(2:N1),kT*(Fa(3:N)+Fp(3:N)),'r');
        plot(z(2:N1),pN*(Fp(3:N)),'m');
        plot(z(2:N1),pN*(Fa(3:N)),'g');
        plot(z(2:N1),kT*Fds(3:N),'c');
        plot(z(2:N1),kT*Fs(3:N),'c--');
        plot(z(2:N1),kT*Fd(3:N),'y');
        plot(z(2:N1),kT*Fc(3:N),'k');
        


        grid;
        legend('F','F_{\alpha+P}','F_{P}','F_{\alpha}','F_{elast}','F_s','F_d','F_C');
        
        xlabel('x ({\mu}m)');
        ylabel('FkT/k_C');
        
        
        add_background_conf(g,zmin,zmax)
    case 1
    case 2


end;


function add_background_conf(g,zmin,zmax);

r = g(:,1);
z = g(:,2);
v        = axis;
xlim_val = xlim;
ylim_val = ylim;


xrange = xlim_val(2)-xlim_val(1);
yrange = ylim_val(2)-ylim_val(1);

xcenter = mean(xlim_val);
ycenter = mean(ylim_val);


Rmax = max(abs(r));
rplotter = r/Rmax*yrange/2;

plot(z,ycenter+rplotter,'k:');
plot(z,ycenter-rplotter,'k:');


plot(zmin+[0 0],ylim_val,'k');
plot(zmax+[0 0],ylim_val,'k');

