PROG-vincenty

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

• 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  

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