緯度・軽度から距離と方位角を求める。

*https://ja.wikipedia.org/wiki/Vincenty%E6%B3%95

を元に作成（Octave で実行)

 %% 
 format long
 invrad = 180.0/pi ;
 a = 6378137.06 ;      
 f = 1/298.257223563 ; 
 b = ( 1 - f ) * a ;
 
 % 始点の緯度・経度
 phi1 =  34.69139  / invrad ;
 L1  =  135.18306  / invrad ;
 % 終点の緯度・経度
 phi2 =  35.7766  / invrad ; 
 L2  = 140.3187  / invrad ;
 
 U1 = atan((1-f)*tan(phi1)) ;
 U2 = atan((1-f)*tan(phi2)) ;
 L=L2 - L1
 
 ramda1 = L1 ;
 ramda2 = L2 ;
 
 ramda = L ;
 FF = abs( L ) 
 sinU1    = sin(U1)    ;
 sinU2    = sin(U2)    ;
 cosU1    = cos(U1)    ;
 cosU2    = cos(U2)    ;
 
 while FF > 1e-12
   ramdaOO = ramda ;
   sinRamda = sin(ramda) ; 
   cosRamda = cos(ramda) ;
   sinSigma   = sqrt(((cosU2*sinRamda)**2)+((cosU1*sinU2-sinU1*cosU2*cosRamda)**2)) ;
   cosSigma    = sinU1*sinU2+cosU1*cosU2*cosRamda ;
   sigma = atan2(sinSigma,cosSigma) ;
   sinAlfa   = ( cosU1*cosU2*sinRamda ) / sinSigma ;
   cos2Alfa = 1 - sinAlfa ** 2 ;
   cos2sigmam = cosSigma - (2*sinU1*sinU2)/cos2Alfa ;
   C    = (f/16) * cos2Alfa * (4 + f * (4 - 3 * cos2Alfa)) ;
   ramda   = L + (1-C)*f*sinAlfa*(sigma+C*sinSigma*(cos2sigmam+C*cosSigma*(-1+2*(cos2sigmam**2)))) ;
   
   FF = abs ( ramdaOO - ramda ) ;
 endwhile
 u2 = cos2Alfa*(a**2 - b**2)/(b**2);
 k1 = (sqrt(1+u2)-1)/(sqrt(1+u2)+1) ;
 A = ( 1+(1/4)*(k1**2) ) / ( 1-k1 ) ;
 B = k1*(1-(3/8)*(k1**2)) ;
 DeltaSigma = B * sinSigma*( cos2sigmam +(1/4)*B* (cosSigma*(-1+2*(cos2sigmam**2))-(1/6)*B*cos2sigmam*(-3+4*(sinSigma**2))*(-3+4*(cos2sigmam**2)))) ;
 
 s = b * A * ( sigma - DeltaSigma );
 % alfa1 = atan( (cosU2*sinRamda)/((cosU1*sinU2)-(sinU1*cosU2*cosRamda)) );
 % alfa2 = atan( (cosU1*sinRamda)/((cosU1*sinU2*cosRamda)-(sinU1*cosU2)) ) + pi ;
 Alfa1 = atan2((cosU2 * sinRamda),((cosU1 * sinU2)-(sinU1 * cosU2 * cosRamda)) );
 Alfa2 = atan2((cosU1*sinRamda),((cosU1*sinU2*cosRamda)-(sinU1*cosU2)) ) + pi ;
 
 % 距離
 s
 % 始点の方位角
 Alfa1 * invrad
 % 終点の方位角
 Alfa2 * invrad  

λのスペルが間違っている。。（修正はその内)