function [p sol Lat1X Lat2X wpts wptsg vertex dis1 dis az1 az2 bv ns yyy]= GeoidesicInverseV1(Lat11,Long1,Lat22,Long2,e )
%GeoidesicInverseV1: Distance, initial and azimuths between points on a shperoid
%   Detailed explanation goes here
%dest1=[-32-59.998/60 -71-36.675/60]
%dep1=[-33-46.21/60 151+31.964/60]
%dep1=[34+26.178/60 139+51.39/60]
%[Lat1X Lat2X wpts wptsg vertex dis1 dis az1 az2 bv ns yyy]= MyGeoidesicInverse1(30,179.8,-29.9,0,0.0818191908426215 )
%[Lat1X Lat2X wpts wptsg vertex dis1 dis az1 az2 bv ns yyy]= MyGeoidesicInverse1(-30,0,30,179+2/3,0.0818191908426215 )
%[wpts dis az1 az2 bv ns yyy]= MyGeoidesicInverse1(-32-59.998/60,-71-36.675/60,34+26.178/60 ,139+51.39/60,0.0818191908426215 )
% (34.4363,139.8565,-32.9999666666667,-71.61125)


ddlon12=wrapnpi2pi(Long2-Long1);
swapLon0=sign(ddlon12);%Lon sign
dlon12=abs(ddlon12);
%convert geodetic to reduced latitude
Lat1X=atand(tand(Lat11)*sqrt(1-e^2));
Lat2X=atand(tand(Lat22)*sqrt(1-e^2));

swapLon=1
if abs(Lat1X)<abs(Lat2X)
Lat1=Lat2X
Lat2=Lat1X
swapLon=-1
else
Lat1=Lat1X
Lat2=Lat2X  
end

swapLat=1
if Lat1>0
    swapLat=-1  
end
Lat2=sign(-Lat1)*Lat2
Lat1=-abs(Lat1)


bv1=abs(Lat1);

%ns=dlon12-rad2deg((GeodesicLat2Lon(Lat2,bv1,e)-GeodesicLat2Lon(-bv1,bv1,e)));
%if Lat2>0
ew=rad2deg(GeodesicLat2Lon(Lat2,bv1,e))
if bv1>0
ew=ew-rad2deg(GeodesicLat2Lon(-bv1,bv1,e))  
end 

%else
%    ew=180*sqrt(1-e^2)
%end

ns=dlon12-ew;
if abs(ns)<10^(-10);
    bv=abs(Lat1);
    return
end 

%tanbv= c=sqrt(tand(Lat1)^2+tand(Lat2)^2-2*tand(Lat1)*tand(Lat2)*cosd(dlon12))/sind(dlon12)
%w21=dlon12/(1-e^2*((cosd(Lat1)+cosd(Lat2))/2)^2)^(1/2)

if ns<0
   dlon12w=dlon12/sqrt((1-e^2))
   c=sind(dlon12w)/sqrt((tand(Lat1)^2+tand(Lat2)^2-2*tand(Lat1)*tand(Lat2)*cosd(dlon12w)+sind(dlon12w)^2))
   p=[]
   sol=[]
else
    %c=2*(180-dlon12)/e^2/180
    %c=sind(161.89052473633)*cos(Lat1)
    %c=sind(161.914)*cos(Lat1) 
    p=[]
   sol=[]
    b1=Lat1
    db=abs(Lat2+Lat1)
    AR=180-dlon12
   if db>eps    
    p=[cosd(b1)^2*e^4*180^2,...
      -4*cosd(b1)^2*AR*e^2*180,...
      4*db^2-cosd(b1)^4*e^4*180^2+4*cosd(b1)^2*AR^2,...
      4*cosd(b1)^4*AR*e^2*180,...
      -4*cosd(b1)^4*AR^2]
    sol=roots(p)
    c=sol((sol> 0) & (sol<(AR*2/e^2/180)))
    c=(AR*2/(180*e^2))/2
    del0=inf
    i=1
    while  i<10 && abs(del0)>1e-9
        del0=p*(c.^(4:-1:0))'/(p(1:4).*(4:-1:1)*(c.^(3:-1:0))')
        c=c-del0
        i=i+1   
    end
   
    %dlon12w=dlon12
   %c=sind(dlon12w)/sqrt((tand(Lat1)^2+tand(Lat2)^2-2*tand(Lat1)*tand(Lat2)*cosd(dlon12w)+sind(dlon12w)^2))
   %c=cosd(56.4509636593606 )
   %c=cosd(abs(Lat1)+0.1)   
   else
        c=2*AR/e^2/180
        %c=cosd(56.4509636593606 )
    end
end



%c=sind(w21)/sqrt((tand(Lat1)^2+tand(Lat2)^2-2*tand(Lat1)*tand(Lat2)*cosd(w21)+sind(w21)^2))
bv=acosd(c)
%bv2=acosd(c)
%bv=74.3705928163273   74.3705203342286  74.3705203342312 74.3705203405886 74.3705203342279
%bv=74.3705203342312  74.3705203340971 74.3705203342221 74.3705203338976
%74.3705203342261  74.3705203325112   74.370514657566  
%if ns>0   
%bv=bv1+(dlon12-ew)/(180-ew)*(90-bv1)
%else
%bv=bv1+(ew-dlon12)/ew*(90-bv1)
%end
%bv=74.3705203342278 
%bv=74.39060253586
%bv=abs(Lat1)
yyy=[];
%eps1=eps+1
%eps1=(180-dlon12)/1e7
delta=inf
%delta1=inf
%delta2=inf
i=1
%bv2=bv+eps1

while abs(delta)>1e-16 && (i<20)
      if ns > 0 
         t1=-2*bv-Lat1; 
         t2=Lat2;    
%         t11=-2*acosd(cosd(bv)+eps1)-Lat1;
      else
         t1=Lat1; 
         t2=Lat2;
 %        t11=t1;
      end
       %acosd(cosd(bv)-eps1)
      %Lon1=LongGeoidesic1(t1,bv,e)
      %Lon2=LongGeoidesic1(t2,bv,e)
      %DLon1=(Lon2-Lon1)
      %DLon1=(LongGeoidesic1(t2,bv,e)-LongGeoidesic1(t1,bv,e))/eps1
      %DLon2=(LongGeoidesic1(t2,acosd(cosd(bv)+eps1),e)-LongGeoidesic1(t11,acosd(cosd(bv)+eps1),e))
      
    %diffDLonX=diffDLon(t2,bv,e)-diffDLon(t1,bv,e)
    diffDLonX= DerivativeDLon(t2,bv,e)-DerivativeDLon(t1,bv,e)
      
      %diffDLonX=(DLon2-DLon1)/eps1
      %bv=bv-(Lon2-Lon1-dlon12)/(DLon2-DLon1)
      %delta=(Lon2-Lon1-dlon12)
      %eps1=10^-2*power(10,log(delta))
      DLonc=LongGeoidesic1(t2,bv,e)-LongGeoidesic1(t1,bv,e)
      delta=DLonc-dlon12
     % delta2=abs(delta1-delta)
     % delta1=delta         
      %c=cosd(bv)
      bv1=bv
      c=c-delta/diffDLonX
      bv=acosd(c)      
      %return
      %bv=bv+rad2deg((Lon2-Lon1-dlon12)/(DLon2-DLon1))
      %bv=bv-(Lon2-Lon1-dlon12)/(DLon2-DLon1)
      i=i+1
      az1=asind(c/cosd(Lat1))
      az2=asind(c/cosd(Lat2))
      if ns>0 
         az1=180-az1
         az2=az2
      end
       if swapLat==-1 
         az1=180-az1
         az2=180-az2
       end 
        if swapLon==-1 || swapLon0==-1
          temp=az1
          az1=-az2
          az2=-temp          
       end
     if swapLon0==-1
        %az1=-az1 
        %az2=-az2
     end
       
     
     % if ns > 0 
     %    t1=-2*bv-Lat1; 
     %    t2=Lat2;    
     % else
    %     t1=Lat1; 
     %    t2=Lat2;    
    %  end
      [ kk dis2]= LongGeoidesicDist(t2,bv,e )
       [kk dis1]=LongGeoidesicDist(t1,bv,e )
       dis=dis2-dis1
       yyy=[yyy;[swapLat swapLon swapLon0 dis az1 az2 bv bv1 delta i diffDLonX DLonc]]       
end


[ Lat Lon az] = LongGeoidesic(linspace(t1,t2,150)',bv,e )
axesm ('mercator','MapLatLimit',[-bv-5 bv+5],'MapLonLimit',[min([Long1,Long1+ddlon12])-20-60 max([Long1,Long1+ddlon12])+20]+60);

%axesm ('mercator','MapLatLimit',[Lat1 Lat2],'MapLonLimit',[Long1 Long1+ddlon12]);

gridm on; plabel on; mlabel on;
if Lat1X>0
    yyy(:,2)=180-yyy(:,2);
   yyy(:,3)=180-yyy(:,3);
end
Lat=atand(tand(Lat)/sqrt(1-e^2))
if swapLon==-1
    Lon=flipud(Lon)
    Lat=flipud(Lat)
end
wptsg=[Lat*swapLat (Lon-Lon(1))*swapLon*swapLon0+Long1 (swapLat==-1)*180+sign(swapLat)*az]
plotm(Lat*swapLat,(Lon-Lon(1))*swapLon*swapLon0+Long1,'r')
wpts=1
[wpts vertex]=wptsgreatellipse_long2([Lat11 Long1],[Lat22 Long2],[6378137 0.0818191908426215])  
 plotm(wpts(:,1:2),'k')
 dis1=wpts(end,3)
 %textm(wpts(1,1),wpts(1,2),'Hokohama-Japan')
%textm(wpts(end,1),wpts(end,2),'Valparaiso-Chile')

load coast
plotm(lat,long)
vertex(2)=atand(tand(vertex(2)*sqrt(1-e^2))) 

%if overhalf==0
%%Lon=p*dlonv+rad2deg(GeodesicLat2Lon(dt,c,e))
%%else
%Lon=(p+sign(t))*dlonv-rad2deg(GeodesicLat2Lon(sign(t)*bv-dt,c,e))  
%end
%Lat=asind(sind(t./bv*90))*bv./90
end

function  [diff]=diffDLon(t,bv,e)
eps1=-1e-12
Lonx1=LongGeoidesic1(t,bv,e)
if abs(t)<abs(bv)
   Lonx2=LongGeoidesic1(t,acosd(cosd(bv)+eps1),e)
else
   t=t+2*bv-2*acosd(cosd(bv)+eps1)
   Lonx2=LongGeoidesic1(t,acosd(cosd(bv)+eps1),e) 
end
diff=(Lonx2-Lonx1)/eps1
end
function  [diffx]=DerivativeDLon(t,bv,e)

if abs(t)<=abs(bv)
   diffx=DGeodesicLat2Lon( t,bv,e )
else
    cycle=-2*DGeodesicLat2Lon(bv,bv,e )
    diffx=cycle+DGeodesicLat2Lon(t+2*bv,bv,e )
end
diffx=rad2deg(diffx)
end
