function [Em ctp err]=afmcntp(z,d,tipf,cntpr,plt)
%% analyze force indenation curve. with specified contact point and
%% indentation value.
% z and d are required inputs for force curves
% maxind is the maximum indentation amount for elasticity
%tipf is tip paraments
%cntpr is a two element vector cntpr(1) is the contact
%point, cntpr(2) is indentation range.
switch nargin
    case {0,1}
        disp('Not enough inputs')
        return
    case 2
        tipf.type=input('Tip type? cone/bead [cone]: ','s');
        if isempty(tipf.type)
            tipf.type='cone';
        end
        tipf.k=input('Force constant k of tip (N/m): ');
        if tipf.type=='cone'|tipf.type=='Cone';
            tipf.geom=input('Angle of the cone (degree): ');
             pbeta=2;
        else
            tipf.geom=input('Diameter of the bead (nm): ');
             pbeta=1.5;
        end
        cntp=input('Contact point: ');
        endp=input('End point: ');
        plt=1;
    case 3
        figure(1),plot(d);
        cntp=input('Contact point: ');
        endp=input('Maximum Indentation: ');
        if tipf.type=='cone'|tipf.type=='Cone';
             pbeta=2;
        else
             pbeta=1.5;
        end
        plt=1;
    case 4
        cntp=cntpr(1);
        endp=cntpr(2);
        if tipf.type=='cone'|tipf.type=='Cone';
             pbeta=2;
        else
             pbeta=1.5;
        end
        plt=1;
    otherwise
        cntp=cntpr(1);
        endp=cntpr(2);
        if tipf.type=='cone'|tipf.type=='Cone';
             pbeta=2;
        else
             pbeta=1.5;
        end
        
end

% Flatten the curve before the contact point use linear fit.
 pf=polyfit(z(1:cntp),d(1:cntp),1);
 d=d-polyval(pf,z);
if max(d)<1e-3
    d=d*1e9;
    z=z*1e9;
end
w=z-d;

%nhertz=round(maxind/mean(diff(w)));
%mbnd=min([cntp+nhertz length(w)]);

wg=w(cntp:end)-w(cntp);
ids=find(wg>=endp);
if isempty(ids)==1
w_fit=w(cntp:end)-w(cntp);
d_fit=d(cntp:end)-d(cntp);
else
w_fit=w(cntp:cntp+ids(1)-1)-w(cntp);
d_fit=d(cntp:cntp+ids(1)-1)-d(cntp);
end
a=1e-4;
emu=fminsearch(@(g)errsum(d_fit,w_fit,pbeta,g),a);

x = linspace(min(w_fit),max(w_fit),1000);
d_hz=hertm(w_fit,pbeta,emu);
d_hz_plot = hertm(x,pbeta,emu);

%calculate Elastic modulus at the rough estimation
 if pbeta==2
 E=emu*(1-0.5^2)*pi/2/tan(tipf.geom/180*pi)*tipf.k*1e9;
 else
 E=emu*(1-0.5^2)*3/4/(tipf.geom/2)^0.5*tipf.k*1e9;
 end

%calculate Elastic modulus with maximum indentation
if plt==1
h=figure(2);
set(h,'Position',[1 100 450 800])
figure(2)

subplot(2,1,1),plot(z,d);hold on,plot(z(cntp),d(cntp),'ro')
xlabel('z (nm)'); ylabel('d (nm)');
grid on
hold off
title(['E=',num2str(E),' Pa']);
subplot(2,1,2),plot(w_fit,d_fit,'--ro');
hold on,plot(x,d_hz_plot);
xlabel('\delta (nm)'); ylabel('d (nm)');
grid on
hold off

end
    Em=E;
    ctp=cntp;
    %if E<10
    %err.rel=1e20;
    %else
    %d_fit
    err.rel=mean((d_hz-d_fit).^2)/abs(range(d_fit)^2);
    %end
    err.tot=sum((d_hz-d_fit).^2);
    err.rms=mean((d_hz-d_fit).^2);
    err.data=[d_fit d_hz w_fit];
    err.n = length(d_fit);


function hres=hertm(din,gbeta,falph)
hres=falph*din.^(gbeta);

function ersum=errsum(a,b,c,d)
ersum=sum(((a-hertm(b,c,d))).^2);

% unbounded search

        
        