
function [c iter_num k]=NewtonMethod(bet1,bet2,lon12,ns,c,e)
% NewtonMethod computes  Clairaut constant c  between points on an ellipsoid

PV=(ns>0); %Pass Vertex
Z=bet1*0;
dc=Z*0+1;
iter_num=Z;
dlon=lon12*0+1;

od1=ceil(log2(abs(c)))+1;
od2=ceil(log2(abs(lon12)));
eps1=2.^od1*eps; 
eps2=max(2.^od2*eps,eps*2^4)*2^4; 
lon12c=lon12;
difflon12=Z*0;
cont=~Z;

k=[];
for i=0:20
      
    cont=cont & abs(dc)>eps1; 
    if ~any(cont), break, end 
    iter_num(cont)=i;              
            
       
    difflon12(cont)=dlongitude12(bet1(cont),bet2(cont),c(cont),PV(cont),e);
    lon12c(cont) =longitude12(bet1(cont),bet2(cont),c(cont),PV(cont),e);  

    dlon(cont)=lon12c(cont)-lon12(cont);
    cont=cont & abs(dlon)>eps2;
     if ~any(cont), break, end 
    delta=dlon(cont)./difflon12(cont);            
    slowcon1=abs((dc(cont)+delta))<eps/2^6; 
    slowcon2=abs((dc(cont)-delta))<eps/2^6;
    delta(slowcon1)=0;
    delta(slowcon2)=0;
    
    c(cont)=c(cont)-delta;
    dc(cont)=delta;
    
    criticalc=cosd(bet1(cont));
    ct=c(cont);
    falure=c(cont)>criticalc;
    ct(falure)=criticalc(falure)-(i+2)*eps;               
    c(cont)=ct;
    
%     k=[k;dc c dlon lon12 lon12c 1 eps1 eps2];
      
end


end