
function [c iter_num t1]=NewtonMethodc(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;
sbet1=sind(bet1);
cbet1=cosd(bet1);
%[sbet1,cbet1]=SinCosNorm(sbet1,cbet1);
sbet2=sind(bet2);
cbet2=cosd(bet2);
%[sbet2,cbet2]=SinCosNorm(sbet2,cbet2);     
M1E = seriesM1E(e);
M2X1 = seriesM2X(sbet1,e);
M2X2=seriesM2X(sbet2,e);
I5=ones(5,1);
%k=[];
tic
for i=0:20
      
    cont=cont & abs(dc)>eps1; 
    if ~any(cont), break, end 
    iter_num(cont)=i;              
    cc=c(cont);
    sbet1c=sbet1(cont);
    cbet1c=cbet1(cont);
    sbet2c=sbet2(cont);
    cbet2c=cbet2(cont);
    PVc=PV(cont);
    
    %nrow=size(cc);    
    %C=repmat(cc,1,5).^repmat((1:2:9),nrow,1);
    C=[cc.^1 cc.^3 cc.^5 cc.^7 cc.^9]; 
    %Cdc=repmat(cc,1,5).^repmat((0:2:8),nrow,1).*repmat((1:2:9),nrow,1);
    Cdc=[cc.^0 3*cc.^2 5*cc.^4 7*cc.^6 9*cc.^8];
    
    u1=real(sqrt(cbet1c.^2-cc.^2));
    u2=real(sqrt(cbet2c.^2-cc.^2));
    sig12=atan2(sbet2c,u2)-(-1).^PVc.*atan2(sbet1c,u1);
   % sig12=asin(sbet2c./sqrt(1-cc.^2))-(-1).^PVc.*asin(sbet1c./sqrt(1-cc.^2));
    %[sig12 sig121 iter_num(cont)]
    omg12=atan2(cc.*sbet2c,u2)-(-1).^PVc.*atan2(cc.*sbet1c,u1);
    M1=C*M1E;
    M21=(C.*M2X1(cont,:))*I5;
    M22=(C.*M2X2(cont,:))*I5;
    M2=M22.*u2-M21.*u1.*(-1).^PVc;
    lon12x=PVc.*(1+M1)*pi+omg12+M1.*sig12+M2;
    lon12x(bet1==0 & c==1)=pi*sqrt(1-e^2);
    lon12c(cont)=rad2deg(lon12x);   
    sig1=atan2(sbet1c,u1);
    sig2=atan2(sbet2c,u2);
    %sig12=sig2-(-1).^PVc.*sig1;
    dsig12=cc./(1-cc.^2).*2.*sin(sig12)./(cos(sig2-sig1)+cos(sig2+sig1));
    dw12=dsig12./cc;
    du1=-cc./u1;
    du2=-cc./u2;       
    dM11=Cdc*M1E;
    dM21=Cdc.*M2X1(cont,:)*I5;
    dM22=Cdc.*M2X2(cont,:)*I5;    
    diff=PVc.*dM11*pi+dw12+dM11.*sig12+M1.*dsig12...
    +dM22.*u2-dM21.*u1.*(-1).^PVc ...
    +M22.*du2-M21.*du1.*(-1).^PVc;
    difflon12(cont)=rad2deg(diff);
    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;
    %falure
%     k=[k;dc c dlon lon12 lon12c 1 eps1 eps2];
      
end
t1=toc;
end