
function [alp1 iter_num k1]=NewtonMethod1(bet1,bet2,lon12,ns,alp1,e)
%bv=acosd(c);
Z=bet1*0;
%c=cosd(bv);%Revert Clairnaut constant to avoid the error in iteration
PV=(ns>0); %Pass Vertex
%t1=Z;
%t2=Z;
dalp1=Z+1;
iter_num=Z;
k1=[]
dlon=lon12*0+1;
od2=1
od2=2
eps1=2.^(ceil(log2(abs(alp1))))*eps; 
eps2=2.^(ceil(log2(abs(lon12))))*eps; 
for i=0:20
      
    cont=abs(dalp1)>eps1 & abs(dlon)>eps2;
    if ~any(cont), break, end 
    iter_num(cont)=i;
    
    alp1t=alp1(cont);
    lon12c=lambda12(bet1(cont),bet2(cont),alp1t,PV(cont),e);
   
    difflon12=dlambda12(bet1(cont),bet2(cont),alp1t,PV(cont),e);   
    
    %[difflon12 lon12c lon12 lon12c-lon12 alp1 dalp1 dlon iter_num]
    dlon(cont)=lon12c-lon12(cont);
    delta=dlon(cont)./difflon12;               
    slowcon1=abs((dalp1(cont)+delta))<eps/100; 
    slowcon2=abs((dalp1(cont)-delta))<eps/100;
     
    delta(slowcon1)=0;
    delta(slowcon2)=0;
     dalp1(cont)=delta;
    %delta(slowcon2)=(1+rand)*delta(slowcon2);
    %delta(slowcon1)=(0.5+rand*2/10)*delta(slowcon1);        
    %overbet1=ct>cosd(bet1(cont))
    %ct(overbet1)=cosd(bet1(overbet1))-10*eps;
    alp1t=alp1t-delta;
    alp1(cont)=alp1t;
      %bv(cont)=acosd(ct);         
    %k1=[k1;difflon12 alp1' dalp1' lon12' dlon' ];%dlona' asind(c./cosd(bet1))' ];%dlon'
%    k1=[k1; difflon12 lon12c lon12 lon12c-lon12 alp1 dalp1];
    %dlon'
      
end

alp1=(alp1);
end