- 追加された行はこのように表示されます。
- 削除された行は
このように表示されます。
緯度・軽度から距離と方位角を求める。
*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 = 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
alfa1 * invrad
alfa2 * invrad
% 終点の方位角
Alfa2 * invrad
λのスペルが間違っている。。(修正はその内)