
function [c iter_num k]=vincenty(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;
od1=ceil(log2(abs(c)))+1;
eps1=2.^od1*eps


cont=~Z;

k=[];
sbet1=sind(bet1)
cbet1=cosd(bet1)

%[sbet1,cbet1]=SinCosNorm(sbet1,cbet1);
sbet2=sind(bet2)
cbet2=cosd(bet2)

%[sbet2,cbet2]=SinCosNorm(sbet2,cbet2);          

for i=0:20
      
    cont=cont & abs(dc)>eps1; 
    if ~any(cont), break, end 
    
    iter_num(cont)=i;                                  
    u1=real(sqrt(cbet1(cont).^2-c(cont).^2));
    u2=real(sqrt(cbet2(cont).^2-c(cont).^2));
    sig12=atan2(sbet2(cont),u2)-(-1).^PV(cont).*atan2(sbet1(cont),u1);
    omg12=atan2(c(cont).*sbet2(cont),u2)-(-1).^PV(cont).*atan2(c(cont).*sbet1(cont),u1);
    M1=seriesM1(sbet1(cont),c(cont));
    M21=seriesM2(sbet1(cont),c(cont));
    M22=seriesM2(sbet2(cont),c(cont));
    n=5;
    E=e^2.^(1:n); 
    M2=M22*E'.*u2-M21*E'.*u1.*(-1).^PV(cont);
    
    omg12c=PV(cont).*(M1*E')*pi+omg12./c(cont)+(M1*E'.*sig12+M2)
    PV
    omg12c=rad2deg(omg12c);
    c1=(lon12(cont)-PV(cont)*180)./omg12c;
    dc(cont)=c(cont)-c1;
    c(cont)=c1;
   
     
    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 omg12 lon12 eps1];
      
end


end